ODE 45 with matrices
1 Ansicht (letzte 30 Tage)
Ältere Kommentare anzeigen
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 Kommentare
Antworten (2)
James Tursa
am 10 Sep. 2015
Bearbeitet: James Tursa
am 10 Sep. 2015
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 Kommentare
Star Strider
am 10 Sep. 2015
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 Kommentare
Siehe auch
Kategorien
Mehr zu Ordinary Differential Equations finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!