Using for loop with ODE solver

Hello MATLAB central!

I have written a function that reads Arrhenius equation and uses different parameters in different temperature ranges, like so:

 function dxdT_MSK =MSK_pyrolysis(T_MSK,x_MSK)
    for k = 1:421
    if T_MSK(k,1) < 500
        dxdT_MSK(k,1) = A_MSK(1).*exp(-E_MSK(1)./(R.*T_MSK(k,1))).*(1-x_MSK(k,1))^n_MSK(1);
    elseif (T_MSK(k,1) >= 500) && (T_MSK(k,1) <620)
               dxdT_MSK(k,1) = A_MSK(2).*exp(-E_MSK(2)./(R.*T_MSK(k,1))).*(1-x_MSK(k,1))^n_MSK(2);
    elseif (T_MSK(k,1) >= 620) && (T_MSK(k,1) <690)
        dxdT_MSK(k,1)=A_MSK(3).*exp(-E_MSK(3)./(R.*T_MSK(k,1))).*(1-x_MSK(k,1))^n_MSK(3); 
    elseif T_MSK(k,1) >=690
        dxdT_MSK(k,1)=0;
    end
    end
 end 
 

I would like to solve it with an ODE solver. I have attempted to do it by using this command:

 T_span_MSK = linspace(460,860)';
 E_MSK = []; %matrix of acitvation energies for pure MSK (primary, secondary, and teriary degradations respectively)
n_MSK = []; %matrix of reaction orders for pure MSK
A_MSK =[]; %matrix of preexponential factors for pure MSK
E_MSK(1) = 8.309;
n_MSK(1) =3;
A_MSK(1) = 9.76*10^(-2);
E_MSK(2) = 37.832;
n_MSK(2) = 0.9;
A_MSK(2) = 3.64*10^2;
E_MSK(3)=14.99;
n_MSK(3) = 1;
A_MSK(3) = 1.61;

[T_MSK,x_MSK] = ode45(@MSK_pyrolysis, T_span_MSK,0);

When I attempt to do so, I get this error message:

Index in position 1 exceeds array bounds (must not exceed 1).

Error in MSK_pyrolysis (line 6)
    if T_MSK(k,1) < 500

Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ode45 (line 115)
  odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in Untitled2 (line 13)
[T_MSK,x_MSK] = ode45(@MSK_pyrolysis, T_span_MSK,0);

I would really appreciate some pointers as to what I could improve or change in my code. I know it would be possible to just do it seperately for each temperature range and then use "brute force" to combine them together, but obviously there must be some more sleek way of doing it.

KInd regards.

Answer 1

[Alan Stevens]

Structurally, your code needs to look more like the following:

 T_span_MSK = linspace(460,860)';
%  E_MSK = []; %matrix of acitvation energies for pure MSK (primary, secondary, and teriary degradations respectively)
% n_MSK = []; %matrix of reaction orders for pure MSK
% A_MSK =[]; %matrix of preexponential factors for pure MSK
E_MSK(1) = 8.309;
n_MSK(1) =3;
A_MSK(1) = 9.76*10^(-2);
E_MSK(2) = 37.832;
n_MSK(2) = 0.9;
A_MSK(2) = 3.64*10^2;
E_MSK(3)=14.99;
n_MSK(3) = 1;
A_MSK(3) = 1.61;

[T_MSK,x_MSK] = ode45(@(T_MSK,x_MSK) MSK_pyrolysis(T_MSK,x_MSK,A_MSK,E_MSK,n_MSK), T_span_MSK,0);

 function dxdT_MSK =MSK_pyrolysis(T_MSK,x_MSK, A_MSK,E_MSK,n_MSK)
    R = 8.3145; % J mol^-1 K^-1
    if T_MSK < 500
        dxdT_MSK = A_MSK(1).*exp(-E_MSK(1)./(R.*T_MSK)).*(1-x_MSK)^n_MSK(1);

    elseif (T_MSK >= 500) && (T_MSK <620)
               dxdT_MSK = A_MSK(2).*exp(-E_MSK(2)./(R.*T_MSK)).*(1-x_MSK)^n_MSK(2);

    elseif (T_MSK >= 620) && (T_MSK <690)
        
        dxdT_MSK=A_MSK(3).*exp(-E_MSK(3)./(R.*T_MSK)).*(1-x_MSK)^n_MSK(3); 
    elseif T_MSK >=690
        dxdT_MSK=0;
    end
   
 end 

However,

(i) when x_MSK gets to 1 you start getting complex values because of the value of n_MSK(2).

(ii) you don't give units so I might have used the wrong value for R (which you didn't include).

(iii) you don't have a beta in your rate equations (dx.dT = A/beta*exp(-E/(RT))*(1-x)^n)- is that because that constant is wrapped up in your A values?

[Aga]

Alan,

thank you for your answer!

The values that I have used are based on the results of a study I have read and yes, I noticed that the results obtained with the values listed give out complex numbers for the mass decomposition. Also, I was using kJ, so R is 8.314*10^(-3), and the heating rare is indeed included in the A value. (I realize I should be more careful with littl details like that, but at this stage the code is only for my own use).

Again, thanks for your help.