Parameter estimation using ODE

Question

Rashedul Khondakar 2020-4-7

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/516209-parameter-estimation-using-ode

评论： Star Strider 2020-4-8

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.

>>

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Star Strider 2020-4-8

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/516209-parameter-estimation-using-ode#answer_424686

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

You have defined it as having 6.

7 个评论
显示 5更早的评论隐藏 5更早的评论

Rashedul Khondakar 2020-4-8

在 MATLAB Online 中打开

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