ODE 45 with matrices
2 ビュー (過去 30 日間)
古いコメントを表示
Hello,
I am trying to solve a system of differential equations with ode 45. basically, my vectors should be x_dot = Ax + Bu Where A and B are matrices, and x_dot, x, and u are vectors. All of the components of u are constant.
function G = parte(t,x,u)
V = 150; % L
k = 0.02; % L/mol*min
beta = 0.15; % kJ L^.5 / mol^.5 min
DeltaH = -15; % kJ/mol A
rho = 4.2; % kg/L
cp = 1.2; % kJ/kg K
u = zeros(3,1);
u(1) = 1.4; % L/min
u(2) = 300; % K
u(3) = 40; % mol/L
A = [-u(1)/V, u(1)*DeltaH/(rho*V*cp), beta/(2*rho*V*cp)*x(3)^(-0.5);
0, -(u(1)/V), -2*k*x(2);
0, 2*k*x(2), -u(1)/V];
B = [(u(2)-x(1))/V + (x(2)-u(3))*DeltaH/(rho*V*cp), u(1)/V, -u(1)*DeltaH/(rho*V*cp);
(u(3)-x(2))/V, 0, u(1)/V;
-x(3)/V, 0, 0];
G = A*x + B*u;
end
Well, here is my script that I try to run and plot
Ti = 300; % K
CAi = 40; % mol/L
CPi = 0; % mol/L
[T4,Y4] = ode45(@parte,[0 10],[Ti CAi CPi]);
subplot(1,2,1)
plot(T4,[Y4(:,2),Y4(:,3)])
xlabel('time (minutes)')
ylabel('Concentration (lb mol/ft^{3})')
legend('A','P','location','best')
title('Concentration vs. time')
And my output
ans =
0 300.0000 40.0000 0
0.2500 NaN NaN NaN
0.5000 NaN NaN NaN
0.7500 NaN NaN NaN
1.0000 NaN NaN NaN
1.2500 NaN NaN NaN
1.5000 NaN NaN NaN
1.7500 NaN NaN NaN
2.0000 NaN NaN NaN
2.2500 NaN NaN NaN
2.5000 NaN NaN NaN
2.7500 NaN NaN NaN
3.0000 NaN NaN NaN
3.2500 NaN NaN NaN
3.5000 NaN NaN NaN
3.7500 NaN NaN NaN
4.0000 NaN NaN NaN
4.2500 NaN NaN NaN
4.5000 NaN NaN NaN
4.7500 NaN NaN NaN
5.0000 NaN NaN NaN
5.2500 NaN NaN NaN
5.5000 NaN NaN NaN
5.7500 NaN NaN NaN
6.0000 NaN NaN NaN
6.2500 NaN NaN NaN
6.5000 NaN NaN NaN
6.7500 NaN NaN NaN
7.0000 NaN NaN NaN
7.2500 NaN NaN NaN
7.5000 NaN NaN NaN
7.7500 NaN NaN NaN
8.0000 NaN NaN NaN
8.2500 NaN NaN NaN
8.5000 NaN NaN NaN
8.7500 NaN NaN NaN
9.0000 NaN NaN NaN
9.2500 NaN NaN NaN
9.5000 NaN NaN NaN
9.7500 NaN NaN NaN
10.0000 NaN NaN NaN
I'm wondering how do do ode 45 when you have a system that uses a matrix. Thanks
0 件のコメント
回答 (2 件)
James Tursa
2015 年 9 月 10 日
編集済み: James Tursa
2015 年 9 月 10 日
CPi is 0, so x(3) in your parte function is 0, so x(3)^(-0.5) is inf (the A(1,3) element). This leads to all of your NaN results.
If the equations are correct, then maybe just do the A*x manually (turn x(3)*x(3)^(-0.5) into x(3)^(0.5) or sqrt(x(3)) to avoid that inf*0 in the result. E.g.,
Ax = [-u(1)/V*x(1) + u(1)*DeltaH/(rho*V*cp)*x(2) + beta/(2*rho*V*cp)*x(3)^(0.5);
0 - (u(1)/V)*x(2) - 2*k*x(2)*x(3);
0 + 2*k*x(2)^2 - u(1)/V*x(3)];
:
G = Ax + B*u;
0 件のコメント
Star Strider
2015 年 9 月 10 日
Unless I’m reading the ‘A’ and ‘B’ matrices wrong, your system is linear and time-invariant with constant coefficients. Your ‘G’ equation (where G=dx/dt) integrates (using Laplace transforms) to:
x = expm(A*t)*B*u
You can put that in a for loop, since it has to be evaluated for each value of ‘t’. You will likely have to experiment with that to get the result you want, although you will likely have to take the references to ‘u’ out of ‘A’ to do it.
0 件のコメント
参考
カテゴリ
Help Center および File Exchange で Ordinary Differential Equations についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!