Matlab 'Matrix dimensions must agree' ode23s

0
votes

The following is my code. I try to model PFR in Matlab using ode23s. It works well with one component irreversible reaction. But when extending more dependent variables, 'Matrix dimensions must agree' problem shows. Have no idea how to fix it. Is possible to use other software to solve similar problems? Thank you.

function PFR_MA_length
    clear all; clc; close all;
    function dCdt = df(t,C)  
        dCdt = zeros(N,2);
        dCddt = [0; -vo*diff(C(:,1))./diff(V)-(-kM*C(2:end,1).*C(2:end,2)-kS*C(2:end,1))];
        dCmdt = [0; -vo*diff(C(:,2))./diff(V)-(-kM*C(2:end,1).*C(2:end,2))];
        dCdt(:,1) = dCddt;
        dCdt(:,2) = dCmdt;
    end
    kM = 1;
    kS = 0.5;                           % assumptions of the rate constants
    C0 = [2, 2];                        % assumptions of the entering concentration
    vo = 2;                             % volumetric flow rate
    volume = 20;                        % total volume of reactor, spacetime = 10

    N = 100;                            % number of points to discretize the reactor volume on

    init = zeros(N,2);                  % Concentration in reactor at t = 0
    init(1,:) = C0;                      % concentration at entrance

    V = linspace(0,volume,N)';          % discretized volume elements, in column form

    tspan = [0 20];
    [t,C] = ode23s(@(t,C) df(t,C),tspan,init);
end

'''

2
matlabode

2 Answers

0
votes

You can put a break point on the line that computes dCddt and observe that the size of the matrices C and V are different.

>> size(C)
 ans =
  200     1

>> size(V)
 ans =
  100     1

The element-wise divide operation, ./, between these two variables would then result in the error that you mentioned.

Per ode23s's help, the output of the call to dCdt = df(t,C) needs to be a vector. However, you are returning a matrix of size 100x2. In the next call to the same function, ode32s converts it to a vector when computing the value of C, hence the size 200x1.

0
votes

In the GNU octave interpretation of Matlab behavior, one has to explicitly make sure that the solver only sees flat one-dimensional state vectors. These have to be translated forward and back in the application of the model.

Explicitly reading the object A as flat array A(:) forgets the matrix dimension information, these can be added back with the reshape(A,m,n) command.

    function dCdt = df(t,C)  
        C = reshape(C,N,2);
        ...
        dCdt = dCdt(:);
    end

    ...

    [t,C] = ode45(@(t,C) df(t,C), tspan, init(:));