MATLAB Answers

ODE 45 with matrices

1 view (last 30 days)
Rick
Rick on 10 Sep 2015
Edited: James Tursa on 10 Sep 2015
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

Answers (2)

James Tursa
James Tursa on 10 Sep 2015
Edited: James Tursa on 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;

Star Strider
Star Strider on 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.

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!