1 view (last 30 days)

Show older comments

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

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
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.

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

Start Hunting!