Runge Kutta integration problem - NaN result
7 views (last 30 days)
Show older comments
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 Comments
Answers (1)
James Tursa
on 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 Comments
James Tursa
on 14 Jan 2021
Edited: James Tursa
on 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.
See Also
Categories
Find more on Function Creation in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!