Parameter estimation using ODE

Hi I am trying to do parameter estimation with a code prepared by Star Strider ‘Igor_Moura'. I have added one more differential equation and edited the code. But it does not run. I would really appreciate it if you would help me to solve my problem. Thank you.

function Igor_Moura
% 2016 12 03
% NOTES: 
% 
%   1. The ‘theta’ (parameter) argument has to be first in your 
%       ‘kinetics’ funciton,
%   2. You need to return ALL the values from ‘DifEq’ since you are fitting
%       all the values

function C=kinetics(theta,t)
c0=[100;0;0;0;0];
[T,Cv]=ode45(@DifEq,t,c0);
%
    function dC=DifEq(t,c)
    dcdt=zeros(5,1);
    dcdt(1)= theta(2).*c(2)+theta(4).*c(3)-(theta(1)+theta(3)+theta(7)).*c(1);
    dcdt(2)= theta(1).*c(1)+theta(5).*c(3)-(theta(6)+theta(8)+theta(9)+theta(2)).*c(2);
    dcdt(3)= theta(3).*c(1)+theta(6).*c(2)-theta(5).*c(3);
    dcdt(4)= 3/5.*theta(7).*c(1)-theta(10).*c(4);
    dcdt(5)= 4/5.*theta(7).*c(1)+4/5.*theta(7).*c(2);
    dC=dcdt;
    end
C=Cv(:,1:4);
end



t=[2.58906
7.00914
14.01829
28.03657
42.05486
56.07314];

c=[94.75 3.277 0.284 0.84 1.14
92.32	4.108	0.361	1.238 1.6763
89.77898 5.24996 0.44828 1.70899 2.32190
85.32557 6.61436 0.50974 2.49297 3.26241
82.73129 7.75630 0.55955 3.49189 4.30973
79.71524 8.54231 0.60406 3.98515 4.76588]; 
theta0=[1;1;1;1;1;1];

[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c);

fprintf(1,'\tRate Constants:\n')
for k1 = 1:length(theta)
    fprintf(1, '\t\tTheta(%d) = %8.5f\n', k1, theta(k1))
end

tv = linspace(min(t), max(t));
Cfit = kinetics(theta, tv);

figure(1)
plot(t, c, 'p')
hold on
hlp = plot(tv, Cfit);
hold off
grid
xlabel('Time')
ylabel('Concentration')
legend(hlp, 'C_1(t)', 'C_2(t)', 'C_3(t)', 'C_4(t)', 'C_5(t)', 'Location','N')

end

And getting the following error

Index exceeds the number of array elements (6).

Error in Aris_code/kinetics/DifEq (line 16)

dcdt(1)= theta(2).*c(2)+theta(4).*c(3)-(theta(1)+theta(3)+theta(7)).*c(1);

Error in odearguments (line 90)

f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode45 (line 115)

odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in Aris_code/kinetics (line 12)

[T,Cv]=ode45(@DifEq,t,c0);

Error in lsqcurvefit (line 214)

initVals.F = feval(funfcn_x_xdata{3},xCurrent,XDATA,varargin{:});

Error in Aris_code (line 43)

[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c);

Caused by:

Failure in initial objective function evaluation. LSQCURVEFIT cannot continue.

>>

Answer 1

[Star Strider]

Your ‘kinetics’ function has 9 parameters, so ‘theta0’ must have 9 elements.

You have defined it as having 6.

[Rashedul Khondakar]

Thank you very much for your help. I will try again. If solves will comment. Much appreciate it.

[Rashedul Khondakar]

Thank you Star Strider. It has run.

[Star Strider]

As always,, my pleasure!

[Rashedul Khondakar]

Hi,

If I want to put a range of input values of all the reaction rate (theta), for example theta (1) = 0.0015-0.0025 and theta (2) = 0.002-0.0025. what would be the input command for theta0?

theta0=[.001984;.002;.000171;.00015;1;1;.000341;1;1;1];

Thanks in advance.

Regards,

function Igor_Moura
% 2016 12 03
% NOTES: 
% 
%   1. The ‘theta’ (parameter) argument has to be first in your 
%       ‘kinetics’ funciton,
%   2. You need to return ALL the values from ‘DifEq’ since you are fitting
%       all the values

function C=kinetics(theta,t)
c0=[100;0;0;0;0];
[T,Cv]=ode45(@DifEq,t,c0);
%
    function dC=DifEq(t,c)
    dcdt=zeros(5,1);
    dcdt(1)= theta(2).*c(2)+theta(4).*c(3)-(theta(1)+theta(3)+theta(7)).*c(1);
    dcdt(2)= theta(1).*c(1)+theta(5).*c(3)-(theta(6)+theta(8)+theta(9)+theta(2)).*c(2);
    dcdt(3)= theta(3).*c(1)+theta(6).*c(2)-theta(5).*c(3);
    dcdt(4)= 3/5.*theta(7).*c(1)-theta(10).*c(4);
    dcdt(5)= 2/5.*theta(7).*c(1)+2/5.*theta(7).*c(2);
    dC=dcdt;
    end
C=Cv(:,1:5);
end



t=[2.58906
7.00914
14.01829
28.03657
42.05486
56.07314];

c=[94.67042929	3.27466397	0.283768327	1.411034146	1.424156884
92.24398939	4.104443076	0.361063431	2.063880736	2.093690133
89.70079166	5.245389346	0.447888069	2.848294106	2.899949115
85.25125707	6.608597877	0.509300617	4.154932927	4.074595714
82.65923334	7.749544147	0.559065958	5.819791255	5.382649314
79.64582001	8.534870801	0.603537114	6.641879119	5.952365412];

theta0=[.001984;.002;.000171;.00015;1;1;.000341;1;1;1];

[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c);

fprintf(1,'\tRate Constants:\n')
for k1 = 1:length(theta)
    fprintf(1, '\t\tTheta(%d) = %8.5f\n', k1, theta(k1))
end

tv = linspace(min(t), max(t));
Cfit = kinetics(theta, tv);

figure(1)
plot(t, c, 'p')
hold on
hlp = plot(tv, Cfit);
hold off
grid
xlabel('Time')
ylabel('Concentration')
legend(hlp, 'C_1(t)', 'C_2(t)', 'C_3(t)', 'C_4(t)', 'C_5(t)', 'Location','N')

end

[Star Strider]

The lsqcurvefit function allows setting of parameter bounds. The estimated parameters will be within those bounds.

See: Best Fit with Bound Constraints for an example.

[Rashedul Khondakar]

Thank you.

[Star Strider]

As always, my pleasure!