Marc
Posts:
108
Registered:
7/2/10


Re: lsqnonlin curve fitting to the output of ODE45
Posted:
May 22, 2013 12:57 AM


"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.22045e014. > > > > >> 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',1e200,'AbsTol',[1e200 1e200 1e200 1e200 > > > > > 1e200]); > > > > > > > > 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 1e198% accuracy, or > > > > > > > > 0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000001% > > > > accuracy. > > > > > > > > Leaving RelTol at its default value of 1e3, or at most reducing it down to > > > > 1e6 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 > > > > 1e200 and using 100*eps instead. Similarly I would not reduce AbsTol below > > > > 1e10, 1e12 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 1e3 till 1e20, 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 searchspace  or use one of the other odeintegration 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',1e3,'AbsTol',[1e6 1e6 1e6 1e6 1e6]); > [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.

