Runge Kutta integration problem - NaN result

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.

Answer 1

[James Tursa]

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( )?

[Gleb]

Yap, its a stiff system. When using ode45 i got the result, not so nice smooth curve but still a curve not a NaN.

[Gleb]

correction: i got NaN matrix when K1-K4 are transposed, else i recieve "Unable to perform assignment because the size of the left side is 1-by-10 and the size of the right side is 10-by-10.

Error in Gas_system_3/GasPhase1 (line 215)

Yg(i+1,:) =Yg(i,:) + (h/6)*( K1 + 2*K2 + 2*K3 + K4 );

"

[James Tursa]

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]

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]

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.