# Parameter estimation using ODE

12 views (last 30 days)
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.
>>

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.

Star Strider on 8 Apr 2020
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 on 8 Apr 2020
Thank you.
Star Strider on 8 Apr 2020
As always, my pleasure!