Parameter estimation using ODE

2020 4 7
1 个回答
回答已采纳
更新时间：2020 4 8
10 次查看（30 天）

0 个投票

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

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

 采纳的回答

Star Strider
Star Strider 2020-4-8

0 个投票

Your ‘kinetics’ function has 9 parameters, so ‘theta0’ must have 9 elements.
You have defined it as having 6.

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

Thank you very much for your help. I will try again. If solves will comment. Much appreciate it.
Thank you Star Strider. It has run.
Star Strider
Star Strider 2020-4-8
As always,, my pleasure!
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
Star Strider 2020-4-8
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.
Star Strider
Star Strider 2020-4-8
As always, my pleasure!

请先登录，再进行评论。

更多回答（0 个）

类别

在 帮助中心 和 File Exchange 中查找有关 Econometrics Toolbox 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by