ODE 45 with matrices

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

Answer 1

[James Tursa]

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;

Answer 2

[Star Strider]

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.