"James" wrote in message <knit3q$a1d$1@newscl01ah.mathworks.com>...
> "Torsten" wrote in message <knhpdu$na5$1@newscl01ah.mathworks.com>...
> > "James" wrote in message <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>...
> > > > > Dear Matlab Community,
> > > > >
> > > > > I am trying to fit the last vector of a coupled ODE (5 ode) to a set of simulated data by optimizing the set of parameters in the coupled ODE. I could be getting lsqnonlin syntax wrongly, so please advise.
> > > > >
> > > > > This is the file (runode1.m) that runs the simulated data to be fitted:
> > > > >
> > > > > function runode1
> > > > > tspan = 0:0.1:100;
> > > > > init_cond1 = [.1];
> > > > >
> > > > > [t, x] = ode45(@ode_fun1, tspan, init_cond1);
> > > > >
> > > > > save('ode1data');
> > > > >
> > > > > figure; plot(t, x(:, 1),'b');;
> > > > >
> > > > > end
> > > > >
> > > > > function dxdt = ode_fun1(t, x)
> > > > > dxdt = zeros(1, 1);
> > > > > n=2;
> > > > > k=.18; %
> > > > > beta=0.36;
> > > > >
> > > > > dx1dt = beta * (x(1)^n / (k^n + x(1)^n)) 0.18*x(1);
> > > > > dxdt(1) = dx1dt;
> > > > > end
> > > > >
> > > > >
> > > > > The data to be fitted can be loaded from ode1data as (t,x).
> > > > > Next this is the coupled ODE system (fit_simp.m). I would like to fit the last vector, in this case sol(:,5) to the simulated data above (t,x) optimizing the param.
> > > > >
> > > > >
> > > > > function diff = fit_simp(param)
> > > > >
> > > > > init_cond2 = [2 2 0.0015 0.0005 0.1];
> > > > > tspan = 0:0.1:100;
> > > > >
> > > > > global kf1 kr2 kfd1 krd1 kf2
> > > > >
> > > > > kf1=param(1);
> > > > > kr2=param(2);
> > > > > kfd1=param(3);
> > > > > krd1=param(4);
> > > > > kf2=param(5);
> > > > >
> > > > > sol=ode45(@myeqn, tspan, init_cond2);
> > > > > diff = sol(:,5);
> > > > >
> > > > > function dydt=myeqn(y)
> > > > > 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
> > > > >
> > > > >
> > > > > To call lsgnonlin i use the following command:
> > > > > t=tspan;
> > > > > param0=[1 1 1 1 1]';
> > > > > param=lsqnonlin(@fit_simp,param0,[],[],[],t,x);
> > > > >
> > > > > The error seems to point that my number of input arguments is not just sol(:,5). I would appreciate if someone could point out my mistakes?
> > > > > Thanks in advance.
> > > > >
> > > > > Yours sincerely,
> > > > > James
> > > >
> > > > If t and x are vectors of the same length, call lsqnonlin as
> > > >
> > > > param=lsqnonlin(@(param)fit_simp(param,t,x),param0);
> > > >
> > > > and organize fit_simp as
> > > >
> > > > function diff = fit_simp(param,t,x)
> > > >
> > > > global kf1 kr2 kfd1 krd1 kf2
> > > >
> > > > kf1=param(1);
> > > > kr2=param(2);
> > > > kfd1=param(3);
> > > > krd1=param(4);
> > > > kf2=param(5);
> > > >
> > > > init_cond2 = [2 2 0.0015 0.0005 0.1];
> > > > sol=ode45(@myeqn, t, init_cond2);
> > > > diff = sol(:,5)x;
> > > >
> > > > function dydt=myeqn(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
> > > >
> > > > Best wishes
> > > > Torsten.
> > >
> > >
> > >
> > > 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]);
> > > [t,y]=ode45(@myeqn, tspan, init,options);
> > > save('test');
> > > diff = y(:,5);
> > >
> > > 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
> >
> > Didn't you notice in my partial code that you have to return
> > diff=y(:,5)x
> > to lsqnonlin ?
> >
> > Best wishes
> > Torsten.
>
>
> Hi Torsten,
> Yes I did, and i did try it. I know it makes sense as lsqnonlin is minimizing the difference between the ouput of ODE45 and X. I still get 6 itereation without making the difference any smaller. Sorry about my previous fit_simp.m being incomplete.
>
> Here's what I have tried and still work in progress:
>
> %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);
>
> tspan=0:0.1:100;
>
> init = [2 2 0.0015 0.0005 0.1];
> [t,y]=ode45(@myeqn, tspan, init);
> save('test');
> diff = y(:,5)x;
> 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
Dear all,
Thank you all for reading and providing me suggestions/corrections.
I would like to report that I have successfully utilize the lsqnonlin for my problem.
Here I want to recap so that the solution may be used for others in the future.
First you need your data to be fitted to. [t,x] where t and x are 1001x1 matrix respectively. Then you need to set your initial parameter guesstimate =>
param0 = [0.150000000000000
0.0100000000000000
1
0.0100000000000000
2]
Then call the lsqnonlin using:
newparam=lsqnonlin(@(param)fit_simp(param,t,x),param0,[0 0 0 0 0],[]);
Here is the code that works:
%fit_simp.m
function diff = fit_simp(param,t,x)
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 .1];
[t,y]=ode23s(@myeqn, tspan, init);
save('test');
diff = y(:,5)x;
function dydt=myeqn(~,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 hope this solution can help others. Thank you!
