Using for loop with ODE solver

2 ビュー (過去 30 日間)
Aga
Aga 2020 年 12 月 9 日
コメント済み: Aga 2020 年 12 月 11 日
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.

採用された回答

Alan Stevens
Alan Stevens 2020 年 12 月 9 日
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?
  1 件のコメント
Aga
Aga 2020 年 12 月 11 日
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.

サインインしてコメントする。

その他の回答 (0 件)

カテゴリ

Help Center および File ExchangeProgramming についてさらに検索

Community Treasure Hunt

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

Start Hunting!

Translated by