Parameter estimation in ODE
7 views (last 30 days)
Show older comments
Bosnian Kingdom
on 23 Apr 2019
Commented: Star Strider
on 24 Apr 2019
I tried to use Star Strider ‘Igor_Moura’ function to solve same problem adding another one ode equation (without experimental data, just need to solve), but i have problem:
function C=kinetics(theta,t)
c0=[1;0;0;0;0];
[T,Cv]=ode45(@DifEq,t,c0);
%
function dC=DifEq(t,c)
dcdt=zeros(5,1);
dcdt(1)=-theta(1).*c(1)-theta(2).*c(1);
dcdt(2)= theta(1).*c(1)+theta(4).*c(3)-theta(3).*c(2)-theta(5).*c(2);
dcdt(3)= theta(2).*c(1)+theta(3).*c(2)-theta(4).*c(3)+theta(6).*c(4);
dcdt(4)= theta(5).*c(2)-theta(6).*c(4);
dcdt(5)=theta(1).*c(2)-theta(5);
dC=dcdt;
end
C=Cv;
end
So, I just added 5th equation, for example:
dcdt(5)=theta(1).*c(2)-theta(5)
I assuemd that there is no experimental data for 5th equation, and i just have to solve this equation.
I change nothing in script.
t=[0.1
0.2
0.4
0.6
0.8
1
1.5
2
3
4
5
6];
c=[0.902 0.06997 0.02463 0.00218
0.8072 0.1353 0.0482 0.008192
0.6757 0.2123 0.0864 0.0289
0.5569 0.2789 0.1063 0.06233
0.4297 0.3292 0.1476 0.09756
0.3774 0.3457 0.1485 0.1255
0.2149 0.3486 0.1821 0.2526
0.141 0.3254 0.194 0.3401
0.04921 0.2445 0.1742 0.5277
0.0178 0.1728 0.1732 0.6323
0.006431 0.1091 0.1137 0.7702
0.002595 0.08301 0.08224 0.835];
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)', 'Location','N')
And i get this:
Error using lsqcurvefit (line 251)
Function value and YDATA sizes are not equal.
Error in script (line 28)
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c);
I would really appreciate the help.
Thanks in advance,
0 Comments
Accepted Answer
Star Strider
on 24 Apr 2019
You are trying to fit 5 columns of your differential equation to 4 columns of data. That will throw the error you got.
Change the ‘C’ assignment in the ‘kinetics’ objective function to:
C=Cv(:,1:4);
and it works.
So ‘klinetics’ is now:
function C=kinetics(theta,t)
c0=[1;0;0;0;0];
[T,Cv]=ode45(@DifEq,t,c0);
%
function dC=DifEq(t,c)
dcdt=zeros(5,1);
dcdt(1)=-theta(1).*c(1)-theta(2).*c(1);
dcdt(2)= theta(1).*c(1)+theta(4).*c(3)-theta(3).*c(2)-theta(5).*c(2);
dcdt(3)= theta(2).*c(1)+theta(3).*c(2)-theta(4).*c(3)+theta(6).*c(4);
dcdt(4)= theta(5).*c(2)-theta(6).*c(4);
dcdt(5)=theta(1).*c(2)-theta(5);
dC=dcdt;
end
C=Cv(:,1:4);
end
If you then want to estimate Compartment 5, run your differential equation function again with all the estimated parameters to include it and plot it.
I created a version of ‘kinetics’ called ‘kinetics5’ to calculate and plot Compartment 5 as well. That entire additional code (with all the previously estimated parameters, so all this goes after the lsqcurvefit call) is:
function C=kinetics5(theta,t)
c0=[1;0;0;0;0];
[T,Cv]=ode45(@DifEq,t,c0);
%
function dC=DifEq(t,c)
dcdt=zeros(5,1);
dcdt(1)=-theta(1).*c(1)-theta(2).*c(1);
dcdt(2)= theta(1).*c(1)+theta(4).*c(3)-theta(3).*c(2)-theta(5).*c(2);
dcdt(3)= theta(2).*c(1)+theta(3).*c(2)-theta(4).*c(3)+theta(6).*c(4);
dcdt(4)= theta(5).*c(2)-theta(6).*c(4);
dcdt(5)=theta(1).*c(2)-theta(5);
dC=dcdt;
end
C=Cv;
end
for k1 = 1:length(theta)
fprintf(1, '\t\tTheta(%d) = %8.5f\n', k1, theta(k1))
end
tv = linspace(min(t), max(t));
Cfit = kinetics5(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','SW')
That should do what you want.
12 Comments
More Answers (0)
See Also
Categories
Find more on Assembly 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!