k3_1の出力は-3.1445e + 24でキャップされます

Question

私はRK4を使ってODEのシステムを解いています。私は、k3_1が-3.1445e + 24でキャップされているという事実によると思われる直線プロットを生成しています。私はなぜそれがキャップされているのか分からない。

function RK4system_MNModel()  
parsec = 3.08*10^18; 
r_1  = 8.5*1000.0*parsec; % in cm 
z_1  = 0.0;    % in cm also 
theta_1 = 0.0; 

grav = 6.6720*10^-8; 
amsun = 1.989*10^33;  % in grams 
amg  = 1.5d11*amsun;  % in grams 
gm  = grav*amg;   % constant 
q  = 0.9;    % axial ratio 

u_1  = 130.0;  % in cm/sec 
w_1  = 95*10^4.0;   % in cm/sec 
v  = 180*10^4.0;  % in cm/sec 
vcirc = sqrt(gm/r_1);  % circular speed (constant) 

nsteps = 50000; 
deltat = 5.0*10^11;   % in seconds 

angmom = r_1*v;    % these are the same 
angmom2 = angmom^2.0; 
e  = -gm/r_1+u_1*u_1/2.0+angmom2/(2.0*r_1*r_1); 

time=0.0; 
for i=1:nsteps 

k3_1 = deltat*u_1 %%%%% THIS LINE 
k4_1 = deltat*(-gm*r_1/((r_1^2.0+(1+sqrt(1+z_1^2.0))^2.0)^1.5) +     angmom2/(r_1^3.0)); % u'=-dphi_dr+lz^2/(r^3.0) with lz=vi*ri this gives deltau 
k5_1 = deltat*(angmom/(r_1^2.0));   % theta'=lz/r^2 this gives deltatheta   
k6_1 = deltat*w_1; 
k7_1 = deltat*(-gm*z_1*(1+sqrt(1+z_1^2.0))/(sqrt(1+z_1^2.0)*(r_1^2.0+(1+sqrt(1+z_1^2.0))^2.0)^1.5)); 
r_2  = r_1+k3_1/2.0; 
u_2  = u_1+k4_1/2.0; 
theta_2 = theta_1+k5_1/2.0; 
z_2  = z_1 + k6_1/2.0; 
w_2  = w_1 + k7_1/2.0; 

k3_2 = deltat*u_2; 
k4_2 = deltat*(-gm*r_2/((r_2^2.0+(1+sqrt(1+z_2^2.0))^2.0)^1.5)+angmom2/(r_2^3.0)); 
k5_2 = deltat*(angmom/(r_2^2.0));     % theta'=lz/r^2 =====> deltatheta  
k6_2 = deltat*w_2; 
k7_2 = deltat*(-gm*z_2*(1+sqrt(1+z_2^2.0))/(sqrt(1+z_2^2.0)*(r_2^2.0+(1+sqrt(1+z_2^2.0))^2.0)^1.5)); 
r_3  = r_1+k3_2/2.0; 
u_3  = u_1+k4_2/2.0; 
theta_3 = theta_1+k5_2/2.0; 
z_3  = z_1 + k6_2/2.0; 
w_3  = w_1 + k7_2/2.0; 

k3_3 = deltat*u_3;       % r'=u      
k4_3 = deltat*(-gm*r_3/((r_3^2.0+(1+sqrt(1+z_3^2.0))^2.0)^1.5)+angmom2/(r_3^3.0));% u'=-dphi_dr+lz^2/(r^3.0) 
k5_3 = deltat*(angmom/(r_3^2.0));   % theta'=lz/r^2 
k6_3 = deltat*w_3; 
k7_3 = deltat*(-gm*z_3*(1+sqrt(1+z_3^2.0))/(sqrt(1+z_3^2.0)*(r_3^2.0+(1+sqrt(1+z_3^2.0))^2.0)^1.5)); 
r_4  = r_1+k3_2; 
u_4  = u_1+k4_2; 
theta_4 = theta_1+k5_2; 
z_4  = z_1 + k6_2; 
w_4  = w_1 + k7_2; 

k3_4 = deltat*u_4;       % r'=u      
k4_4 = deltat*(-gm*r_4/((r_4^2.0+(1+sqrt(1+z_4^2.0))^2.0)^1.5)+angmom2/(r_4^3.0)); % u'=-dphi_dr+lz^2/(r^3.0) 
k5_4 = deltat*(angmom/(r_4^2.0));   % theta'=lz/r^2  
k6_4 = deltat*w_4; 
k7_4 = deltat*(-gm*z_4*(1+sqrt(1+z_4^2.0))/(sqrt(1+z_4^2.0)*(r_4^2.0+(1+sqrt(1+z_4^2.0))^2.0)^1.5)); 

r_1  = r_1+(k3_1+2.0*k3_2+2.0*k3_3+k3_4)/6.0; % New value of R for next step 
u_1  = u_1+(k4_1+2.0*k4_2+2.0*k4_3+k4_4)/6.0; % New value of U for next step 
theta_1 = theta_1+(k5_1+2.0*k5_2+2.0*k5_3+k5_4)/6.0; % New value of  theta 
z_1  = z_1+(k6_1+2.0*k6_2+2.0*k6_3+k6_4)/6.0; 
w_1  = w_1+(k7_1+2.0*k7_2+2.0*k7_3+k7_4)/6.0; 

e  = -gm/r_1+u_1*u_1/2.0+angmom2/(2.0*r_1*r_1); % energy 
ecc  = (1.0+(2.0*e*angmom2)/(gm^2.0))^0.5;  % eccentricity 
x(i) = r_1*cos(theta_1)/(1000.0*parsec);  % X for plotting orbit 
y(i) = r_1*sin(theta_1)/(1000.0*parsec);  % Y for plotting orbit 
time = time+deltat; 
r(i) = r_1; 
z(i) = z_1; 
time1(i)= time; 
end 

指定された行に異常が発生することに注意してください。

Answer 1

それに蓋だとk3_1ではない、それは-3.1445e+24/deltatの値を返すu_1の計算（deltatは定数）です。

u_1が一列に計算される：最初の反復後

u_1  = u_1+(k4_1+2.0*k4_2+2.0*k4_3+k4_4)/6.0; 

、これが返さ：

duが比較的小さい差を表すように見える式 u_1(n+1) = u_1(n) + duに基づい

u_1(1) = 6.500e+13    % Hard coded before the loop 
u_1(2) = -1.432966614767040e+04 % Calculated using the equation above 
u_1(3) = -2.878934017859105e+04 % Calculated using the equation above 
u_1(4) = -4.324903004768405e+04 

を。 2つの最初の値の差は非常にです。この計算が間違っていると仮定しています。

あなたはその計算が正確であることが判明した場合、その後、あなたのエラーがこれらの行のいずれかである：彼らはあなたの質問を解決した場合

k4_1 = deltat*(-gm*r_1/((r_1^2.0+(1+sqrt(1+z_1^2.0))^2.0)^1.5)+angmom2/(r_1^3.0)); % u'=-dphi_dr+lz^2/(r^3.0) with lz=vi*ri this gives delta 
k4_2 = deltat*(-gm*r_2/((r_2^2.0+(1+sqrt(1+z_2^2.0))^2.0)^1.5)+angmom2/(r_2^3.0)); 
k4_3 = deltat*(-gm*r_3/((r_3^2.0+(1+sqrt(1+z_3^2.0))^2.0)^1.5)+angmom2/(r_3^3.0));% u'=-dphi_dr+lz^2/(r^3.0) 
k4_4 = deltat*(-gm*r_4/((r_4^2.0+(1+sqrt(1+z_4^2.0))^2.0)^1.5)+angmom2/(r_4^3.0)); % u'=-dphi_dr+lz^2/(r^3.0)