Parameter estimation using ODE
5 views (last 30 days)
Show older comments
Rashedul Khondakar
on 7 Apr 2020
Commented: Star Strider
on 8 Apr 2020
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 Comments
Accepted Answer
Star Strider
on 8 Apr 2020
Your ‘kinetics’ function has 9 parameters, so ‘theta0’ must have 9 elements.
You have defined it as having 6.
7 Comments
More Answers (0)
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!