Solving a system of ODEs coupled with an algebraic equation
2 次查看(过去 30 天)
显示 更早的评论
Hi everyone,
I have tried to couple an fsolve function with a system ODEs found in:
while following Star Strider ‘Igor_Moura’ function.
The system of ODEs are solved by:
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
pHCalculation;
tv = linspace(min(t), max(t));
Cfit = kinetics5(theta, tv);
figure(1)
plot(t, c, 'p')
hold on
hlp = plot(tv, Cfit,t,ph.*0.1,'r-');
hold on
plot(t,ph.*0.1,'r-')
hold on
plot (t,pH.*0.1,'kd')
grid
xlabel('Time')
ylabel('Concentration')
legend(hlp, 'C_1(t)', 'C_2(t)', 'C_3(t)', 'C_4(t)', 'C_5(t)','pH-Mod','Location','best')
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
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
And the algebraic equation is given by:
K1 = 6.24*2;
K2 = 5.68e-5*4e1;
K3 = 1.7e-3*1.1e3;
K4 = 6.55e-8;
K5 = 5.3e-8;
x0 = 7.41e-06;
t = [0.1
0.2
0.4
0.6
0.8
1
1.5
2
3
4
5
6];
A = [0.04879
6.650241034
9.984349357
12.80945024
15.8072316
18.85920165
21.81901476
24.15458161
23.99797436
22.9179308
22.08000295
21.3980477
];
B = [0
3.228125432
3.039568727
2.353698976
1.924485852
1.669257876
1.577120575
1.848544338
3.227063
5.416789376
7.795831728
10.27649137
];
D = [0.04879
6.89871101
13.90837796
20.90460441
27.81415915
34.59918462
41.11877195
46.61855348
48.32509945
47.89047761
47.18901391
46.2708157
];
pH = [8.13
7.76
7.43
7.45
7.45
7.00
4.39
3.61
3.21
2.88
2.88
2.64
];
for i = 1:length(t)
% BosnianKingdom;
% A = c(1);
% B = c(2);
% D = c(3);
sol(i) = fsolve(@(x) x + 2.* A(i) - ((B(i).* K1.*x)/(x.^2 + K1.*x + K1*K2))- 2.*((B(i).*K1*K2)/(x.^2 ...
+ K1.*x + K1*K2))-((D(i).*K3.*x)/(x.^2 + K3.*x + K3*K4))- 2.*((D(i).*K3.*K4)/(x.^2 ...
+ K3.*x + K3*K4))- K5./x, x0);
end
ph = 3 - log10(sol);
plot(t,ph.*0.1,'r--')
hold on
plot (t,pH.*0.1,'kd')
xlabel('Time (sec)')
ylabel('pH')
hold off
legend('pH-Mod','pH-Exp')
The code works well when I solve ODEs seperately and manually input the ODEs-calculated values c(1), c(2) and c(3) as columns of A, B and D. I am interested in inputing c(1), c(2) and c(3) step-wise as they are being calculated from the system of ODEs. Is it possible, or should carryon the way I am doing?
Thanks for your help.
0 个评论
回答(0 个)
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!