MATLAB Newsgroup

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

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

"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

"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

"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

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

"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

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

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

"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

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

You can think of your watch list as threads that you have bookmarked.

You can add tags, authors, threads, and even search results to your watch list. This way you can easily keep track of topics that you're interested in. To view your watch list, click on the "My Newsreader" link.

To add items to your watch list, click the "add to watch list" link at the bottom of any page.

To add search criteria to your watch list, search for the desired term in the search box. Click on the "Add this search to my watch list" link on the search results page.

You can also add a tag to your watch list by searching for the tag with the directive "tag:tag_name" where tag_name is the name of the tag you would like to watch.

To add an author to your watch list, go to the author's profile page and click on the "Add this author to my watch list" link at the top of the page. You can also add an author to your watch list by going to a thread that the author has posted to and clicking on the "Add this author to my watch list" link. You will be notified whenever the author makes a post.

To add a thread to your watch list, go to the thread page and click the "Add this thread to my watch list" link at the top of the page.

The newsgroups are a worldwide forum that is open to everyone. Newsgroups are used to discuss a huge range of topics, make announcements, and trade files.

Discussions are threaded, or grouped in a way that allows you to read a posted message and all of its replies in chronological order. This makes it easy to follow the thread of the conversation, and to see what’s already been said before you post your own reply or make a new posting.

Newsgroup content is distributed by servers hosted by various organizations on the Internet. Messages are exchanged and managed using open-standard protocols. No single entity “owns” the newsgroups.

There are thousands of newsgroups, each addressing a single topic or area of interest. The MATLAB Central Newsreader posts and displays messages in the comp.soft-sys.matlab newsgroup.

**MATLAB Central**

You can use the integrated newsreader at the MATLAB Central website to read and post messages in this newsgroup. MATLAB Central is hosted by MathWorks.

Messages posted through the MATLAB Central Newsreader are seen by everyone using the newsgroups, regardless of how they access the newsgroups. There are several advantages to using MATLAB Central.

**One Account**

Your MATLAB Central account is tied to your MathWorks Account for easy access.

**Use the Email Address of Your Choice**

The MATLAB Central Newsreader allows you to define an alternative email address as your posting address, avoiding clutter in your primary mailbox and reducing spam.

**Spam Control**

Most newsgroup spam is filtered out by the MATLAB Central Newsreader.

**Tagging**

Messages can be tagged with a relevant label by any signed-in user. Tags can be used as keywords to find particular files of interest, or as a way to categorize your bookmarked postings. You may choose to allow others to view your tags, and you can view or search others’ tags as well as those of the community at large. Tagging provides a way to see both the big trends and the smaller, more obscure ideas and applications.

**Watch lists**

Setting up watch lists allows you to be notified of updates made to postings selected by author, thread, or any search variable. Your watch list notifications can be sent by email (daily digest or immediate), displayed in My Newsreader, or sent via RSS feed.

- Use a newsreader through your school, employer, or internet service provider
- Pay for newsgroup access from a commercial provider
- Use Google Groups
- Mathforum.org provides a newsreader with access to the comp.soft sys.matlab newsgroup
- Run your own server. For typical instructions, see: http://www.slyck.com/ng.php?page=2