Drexel dragonThe Math ForumDonate to the Math Forum



Search All of the Math Forum:

Views expressed in these public forums are not endorsed by Drexel University or The Math Forum.


Math Forum » Discussions » Software » comp.soft-sys.matlab

Topic: lsqnonlin curve fitting to the output of ODE45
Replies: 10   Last Post: May 29, 2013 1:07 PM

Advanced Search

Back to Topic List Back to Topic List Jump to Tree View Jump to Tree View   Messages: [ Previous | Next ]
Marc

Posts: 108
Registered: 7/2/10
Re: lsqnonlin curve fitting to the output of ODE45
Posted: May 22, 2013 12:57 AM
  Click to see the message monospaced in plain text Plain Text   Click to reply to this topic Reply

"James" wrote in message <kngvom$h6p$1@newscl01ah.mathworks.com>...
> "Bjorn Gustavsson" <bjonr@irf.se> wrote in message <kngshj$8ta$1@newscl01ah.mathworks.com>...
> > "James" wrote in message <kngk1s$d41$1@newscl01ah.mathworks.com>...
> > > "Steven_Lord" <slord@mathworks.com> wrote in message <kngbv2$et8$1@newscl01ah.mathworks.com>...
> > > >
> > > >
> > > > "James " <hsukiang@gmail.com> wrote in message
> > > > news:kng7ic$cm$1@newscl01ah.mathworks.com...

> > > > > "Torsten" wrote in message <knf5kl$iia$1@newscl01ah.mathworks.com>...
> > > > >> "James" wrote in message <kndt51$rin$1@newscl01ah.mathworks.com>...
> > > >
> > > > *snip*
> > > >

> > > > > Thanks Torsten for the reply. Vectors t and x are of the same length. I
> > > > > tried your recommendation, with a little tweak, it almost work! However
> > > > > an error event took place due to tolerance issue with the integration
> > > > > process due to tolerance. I could not understand why the coupled ODE
> > > > > integrates smoothly when ran by itself but faces tolerance issue when ran
> > > > > with lsqnonlin.
> > > > > Appreciate your inputs on this.
> > > > >
> > > > > Warning: RelTol has been increased to 2.22045e-014.

> > > > >> In funfun\private\odearguments at 128
> > > > > In ode45 at 114
> > > > > In fit_simp at 15
> > > > > In @(param)fit_simp(param,t,x)
> > > > > In lsqnonlin at 199 .
> > > > >
> > > > > Here's the fit_simp.m that I have tweaked from your recommendation:
> > > > >
> > > > > function diff = fit_simp(param,t,x)
> > > > >
> > > > > global kf1 kr1 kfd1 krd1 kf2
> > > > >
> > > > > kf1=param(1);
> > > > > kr1=param(2);
> > > > > kfd1=param(3);
> > > > > krd1=param(4);
> > > > > kf2=param(5);
> > > > >
> > > > > tspan=0:0.1:100;
> > > > >
> > > > > init = [2 2 0.0015 0.0005 0.1];
> > > > > options = odeset('RelTol',1e-200,'AbsTol',[1e-200 1e-200 1e-200 1e-200
> > > > > 1e-200]);

> > > >
> > > > Don't use tolerances this small. By the definitions in the documentation:
> > > >
> > > > http://www.mathworks.com/help/matlab/ref/odeset.html
> > > >
> > > > your RelTol value is asking for 1e-198% accuracy, or
> > > >
> > > > 0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000001%
> > > > accuracy.
> > > >
> > > > Leaving RelTol at its default value of 1e-3, or at most reducing it down to
> > > > 1e-6 or so, should be sufficient for most problems. The warning you received
> > > > was ODE45 telling you it was ignoring your request for it to use a RelTol of
> > > > 1e-200 and using 100*eps instead. Similarly I would not reduce AbsTol below
> > > > 1e-10, 1e-12 without a strong reason and a solution of very small magnitude
> > > > (and even then, I'd probably rescale the problem instead, like using 100
> > > > meters instead of 0.1 kilometers as the quantities in my problem.)
> > > >
> > > > *snip rest of code*
> > > >
> > > > --
> > > > Steve Lord
> > > > slord@mathworks.com
> > > > To contact Technical Support use the Contact Us link on
> > > > http://www.mathworks.com

> > >
> > >
> > > Thanks Steve, I have tried from 1e-3 till 1e-20, the latter I have suspected of "overkill" tolerance set. Having said that, thank for the advice. I will stick to a reasonable tolerance. I think there is some underlying issue or mistake that I have made when I call ode45 from lsqnonlin. I m still trying to fix the problem as we speak. Thank you all for reading and posting useful hint. I m much closer to solving this problem!
> > >

> > Might one problem be (become) that lsqnonlin runs into some corner of parameter space where your ODE becomes stiff (or worse)? If so you might have to constrain the search-space - or use one of the other ode-integration functions.
>
>
> I did try ode23s thinking along the same line as you but it did really help. What actually works (partially) is that I made "param" a global variable. Then lsqnonlin was able to call fit_simp.m.
> However lsqnonlin did not do what it is intend to, ie. to optimize/ or fit to the data t,x supplied while varying the param. I see 6 iterations that lsqnonlin call fit_simp.m and every iteration did nothing but just print the supplied parameters. At the 6th iteration lsqnonlin basically give up, I am thinking of changing the lsqnonlin stopping criteria. Any suggestions?
>
> "Initial point is a local minimum.
>
> Optimization completed because the size of the gradient at the initial point
> is less than the default value of the function tolerance."
>
>
> here's my latest fit_simp.m
> function diff = fit_simp(param,t,x)
>
> global kf1 kr1 kfd1 krd1 kf2 param
>
> kf1=param(1);
> kr1=param(2);
> kfd1=param(3);
> krd1=param(4);
> kf2=param(5);
>
> % param(1)=kf1;
> % param(2)=kr1;
> % param(3)=kfd1;
> % param(4)=krd1;
> % param(5)=kf2;
>
> tspan=0:0.1:100;
>
> init = [2 2 0.0015 0.0005 0.1];
> % options = odeset('RelTol',1e-3,'AbsTol',[1e-6 1e-6 1e-6 1e-6 1e-6]);
> [t,y]=ode45(@myeqn, tspan, init);
> save('test');
> diff = y(:,5);
> figure;plot(t,y(:,5),'g',t,x,'r');
>
> function dydt=myeqn(t,y)
> dydt = zeros(5, 1);
> dydt(1) = -kf1*y(1)*y(2) + kr1*y(3);
> dydt(2) = -kf1*y(1)*y(2) + kr1*y(3);
> dydt(3) = kf1*y(1)*y(2) - kr1*y(3) - 2*kfd1*y(4) + 2*krd1*y(4);
> dydt(4) = 2*kfd1*y(3)^2 - krd1*y(4) - kf2*y(4);
> dydt(5) = kf2*y(4);
> end
> end


I have not had good results with nlinfit, lsqnonlin when fitting systems of odes. I would set this up like an optimization problem, finding the least squares yourself. This way you can also play with the weights. Then try fminsearch. It's slower but the nelder mead simplex does not require taking differentials, so I have found it to converge with much better results.

Also, I would try ode15s. If I have some time tomorrow i will give your system a shot.

I am not sure why you have the params as global. Since you are nesting the functions, this should not be needed.



Point your RSS reader here for a feed of the latest messages in this topic.

[Privacy Policy] [Terms of Use]

© Drexel University 1994-2014. All Rights Reserved.
The Math Forum is a research and educational enterprise of the Drexel University School of Education.