Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
lsqnonlin curve fitting to the output of ODE45

Subject: lsqnonlin curve fitting to the output of ODE45

From: James

Date: 20 May, 2013 19:22:09

Message: 1 of 11

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

Subject: lsqnonlin curve fitting to the output of ODE45

From: Torsten

Date: 21 May, 2013 06:53:09

Message: 2 of 11

"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.

Subject: lsqnonlin curve fitting to the output of ODE45

From: James

Date: 21 May, 2013 16:32:12

Message: 3 of 11

"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.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]);
[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

Subject: lsqnonlin curve fitting to the output of ODE45

From: Steven_Lord

Date: 21 May, 2013 17:47:15

Message: 4 of 11



"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

Subject: lsqnonlin curve fitting to the output of ODE45

From: James

Date: 21 May, 2013 20:05:16

Message: 5 of 11

"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!

James Ooi

Subject: lsqnonlin curve fitting to the output of ODE45

From: Bjorn Gustavsson

Date: 21 May, 2013 22:30:11

Message: 6 of 11

"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.

Subject: lsqnonlin curve fitting to the output of ODE45

From: James

Date: 21 May, 2013 23:25:10

Message: 7 of 11

"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

Subject: lsqnonlin curve fitting to the output of ODE45

From: Marc

Date: 22 May, 2013 04:57:10

Message: 8 of 11

"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.

Subject: lsqnonlin curve fitting to the output of ODE45

From: Torsten

Date: 22 May, 2013 06:43:10

Message: 9 of 11

"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.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]);
> [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.

Subject: lsqnonlin curve fitting to the output of ODE45

From: James

Date: 22 May, 2013 16:52:10

Message: 10 of 11

"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.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]);
> > [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

Subject: lsqnonlin curve fitting to the output of ODE45

From: James

Date: 29 May, 2013 17:07:09

Message: 11 of 11

"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.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]);
> > > [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!

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us