MATLAB Answers

Runge Kutta integration problem - NaN result

조회 수: 35(최근 30일)
Gleb
Gleb 14 Jan 2021
편집: James Tursa 14 Jan 2021
Hi! Here is the code
[x, Yg]=GasPhase1(J_g)
VarNames = { 'T', 'q', 'Y_gHMX', 'Y_cTf', 'Y_gTf', 'Y_CF2','Y_C', 'Y_B', 'Y_BF3', 'fg' };
Table1 = table( Yg(:,1),Yg(:,2),Yg(:,3), Yg(:,4), Yg(:,5), Yg(:,6), Yg(:,7), Yg(:,8), ...
Yg(:,9), Yg(:,10), 'VariableNames',VarNames);
disp(Table1)
function [x, Yg]=GasPhase1(J_g)
xspan = [0 L_g/l];
switch solverselect
case 1
AbsTolerance= [5, 20, ones(1, 8)*10^-1 ];
init_cond =[T_fout, J_g, Y_gHMXout, f1Y_cTfout, 0, 0, 0, f1Y_B0, 0, f_gout ];
options = odeset('RelTol', 10^-1 ,'AbsTol', AbsTolerance, 'NonNegative', 3); %
[x,Yg] = ode15s(@(x,Y) odefcn_g(x,Y), xspan, init_cond, options);
case 2
f = @(x,Yg) odefcn_g(x,Yg);
h=10^-1;
N_st = ceil((L_g/l)/h)+1;
x=(0:h:L_g/l)';
Yg = zeros(N_st, 10);
Yg(1, :)= [T_fout, J_g, Y_gHMXout, f1Y_cTfout, 0, 0, 0, f1Y_B0, 0, f_gout ];
for i = 1:N_st
K1 = f( x(i) , Yg(i,:) );
K2 = f( x(i) + h/2, Yg(i,:) + K1*h/2 );
K3 = f( x(i) + h/2, Yg(i,:) + K2*h/2 );
K4 = f( x(i) + h , Yg(i,:) + K3*h );
Yg(i+1,:) =Yg(i,:) + (h/6)*( K1 + 2*K2 + 2*K3 + K4 );
end
end
function dYdx_g = odefcn_g(x,Yg)
% some vector of RHS's
end
when i use the standard solver 15s everithing is ok, but Runge Kutta method returns NaN elements. I think here is a syntax error.

  댓글 수: 0

댓글을 달려면 로그인하십시오.

답변(1개)

James Tursa
James Tursa 14 Jan 2021
You may have a stiff DE and you can't use RK4 with a fixed step size. What happens when you call ode45( ) instead of ode15s( )?

  댓글 수: 5

표시 이전 댓글 수: 2
James Tursa
James Tursa 14 Jan 2021
I am assuming your odefcn_g( ) function returns a column vector to comply with ode15s( ) and ode45( ) requirements, so yes you would need to transpose the K1-K4 into row vectors when adding to the row vector Yg(i,:). But it looks like ode45( ) is returning NaN, so you probably have no hope of making your manually coded RK4 with a fixed step size work.
Gleb
Gleb 14 Jan 2021
with transposed K1-K4 i ve got
T q Y_gHMX Y_cTf Y_gTf Y_CF2 Y_C Y_B Y_BF3 fg
________________ _______________ ________________ ________ ________ ________ ________ _________ ________ ________
725+0i 100+0i 0.9+0i 0+0i 0+0i 0+0i 0+0i 0.03+0i 0+0i 0.3+0i
-1.327e+05+0i 3.1992e+27+0i -2.3527e+21+0i 0+0i 0+0i 0+0i 0+0i 0.03+0i 0+0i 0.3+0i
-1.0392e+06+0i NaN+NaNi Inf+0i 0+0i 0+0i 0+0i NaN+NaNi 0.03+0i 0+0i 0.3+0i
so i think variables grow incredibly fast after the first iteration, and on 3rd become "infinite". Is it an effects of stiffness? when i changed RHS to more linear, then NaN began to appear later
James Tursa
James Tursa 14 Jan 2021
Your initial condition is complex?
If your first step produces numbers in the e27 range, then my first thought would be that something is definitely wrong with the setup or stepsize etc.

댓글을 달려면 로그인하십시오.

Community Treasure Hunt

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

Start Hunting!

Translated by