MATLAB - 適応ステップサイズルンゲクッタ

私はMATLABでのODEのシステムを解決するための適応ステップサイズRK4をプログラムしました。コードはエラーなく実行されますが、yに対してxをプロットしようとすると、目的のカーブが生成されません。トロイダル型ではなく、単にフラットラインを得るだけです。これは、rが一定値を出力していることから明らかである。各行の出力を確認した後、定数や誤差、またはinfやNaNを出力しないで、実数成分と虚数成分（複素数）の両方を出力します。なぜこのようなことが起こっているのか、私はそれが自分のトラブルの原因であると信じています。あなたが "私" あなたのコード内で使用しているMATLAB - 適応ステップサイズルンゲクッタ

function AdaptRK4() 
parsec = 3.08*10^18; 
r_1  = 8.5*1000.0*parsec; % in cm 
theta_1 = 0.0; 
a  = 0.5*r_1; 
gam  = 1; 

grav = 6.6720*10^-8; 
amsun = 1.989*10^33; 
amg  = 1.5d11*amsun; 
gm  = grav*amg;  

u_1  = 20.0*10^5; 
v  = sqrt(gm/r_1); 

time = 0.0; 
epsilon = 0.00001; 
m1  = 0.5; 
m2  = 0.5; 
m3  = 0.5; 
i  = 1; 
nsteps = 50000; 
deltat = 5.0*10^12; 

angmom = r_1*v; 
angmom2 = angmom^2.0; 
e  = -2*10^5.0*gm/r_1+u_1*u_1/2.0+angmom2/(2.0*r_1*r_1); 

for i=1:nsteps 
deltat = min(deltat,nsteps-time); 

fk3_1 = deltat*u_1; 
fk4_1 = deltat*(-gm*r_1*r_1^(-gam)/(a+r_1)^(3- gam)+angmom2/(r_1^3.0)); 
fk5_1 = deltat*(angmom/(r_1^2.0)); 
r_2  = r_1+fk3_1/4.0; 
u_2  = u_1+fk4_1/4.0; 
theta_2 = theta_1+fk5_1/4.0; 

fk3_2 = deltat*u_2; 
fk4_2 = deltat*(-gm*r_2*r_2^(-gam)/(a+r_2)^(3-gam)+angmom2/(r_2^3.0)); 
fk5_2 = deltat*(angmom/(r_2^2.0));   
r_3  = r_1+(3/32)*fk3_1 + (9/32)*fk3_2; 
u_3  = u_1+(3/32)*fk4_1 + (9/32)*fk4_2; 
theta_3 = theta_1+(3/32)*fk5_1 + (9/32)*fk5_2; 

fk3_3 = deltat*u_3;       
fk4_3 = deltat*(-gm*r_3*r_3^(-gam)/(a+r_3)^(3-gam)+angmom2/(r_3^3.0)); 
fk5_3 = deltat*(angmom/(r_3^2.0));    
r_4  = r_1+(1932/2197)*fk3_1 - (7200/2197)*fk3_2 + (7296/2197)*fk3_3; 
u_4  = u_1+(1932/2197)*fk4_1 - (7200/2197)*fk4_2 + (7296/2197)*fk4_3; 
theta_4 = theta_1+(1932/2197)*fk5_1 - (7200/2197)*fk5_2 + (7296/2197)*fk5_3; 

fk3_4 = deltat*u_4;       
fk4_4 = deltat*(-gm*r_4*r_4^(-gam)/(a+r_4)^(3-gam)+angmom2/(r_4^3.0)); 
fk5_4 = deltat*(angmom/(r_4^2.0));    
r_5  = r_1+(439/216)*fk3_1 - 8*fk3_2 + (3680/513)*fk3_3 -  (845/4104)*fk3_4; 
u_5  = u_1+(439/216)*fk4_1 - 8*fk4_2 + (3680/513)*fk4_3 - (845/4104)*fk4_4; 
theta_5 = theta_1+(439/216)*fk5_1 - 8*fk5_2 + (3680/513)*fk5_3 -  (845/4104)*fk5_4; 

fk3_5 = deltat*u_5;       
fk4_5 = deltat*(-gm*r_5*r_5^(-gam)/(a+r_5)^(3-gam)+angmom2/(r_5^3.0)); 
fk5_5 = deltat*(angmom/(r_5^2.0));    
r_6  = r_1-(8/27)*fk3_1 - 2*fk3_2 - (3544/2565)*fk3_3 + (1859/4104)*fk3_4-(11/40)*fk3_5; 
u_6  = u_1-(8/27)*fk4_1 - 2*fk4_2 - (3544/2565)*fk4_3 + (1859/4104)*fk4_4-(11/40)*fk4_5; 
theta_6 = theta_1-(8/27)*fk5_1 - 2*fk5_2 - (3544/2565)*fk5_3 + (1859/4104)*fk5_4-(11/40)*fk5_5; 

fk3_6 = deltat*u_6;   
fk4_6 = deltat*(-gm*r_6*r_6^(-gam)/(a+r_6)^(3-gam)+angmom2/(r_6^3.0)); 
fk5_6 = deltat*(angmom/(r_6^2.0)); 

fm3_1 = m1 + 25*fk3_1/216+1408*fk3_3/2565+2197*fk3_4/4104-fk3_5/5; 
fm4_1 = m2 + 25*fk4_1/216+1408*fk4_3/2565+2197*fk4_4/4104-fk4_5/5; 
fm5_1 = m3 + 25*fk5_1/216+1408*fk5_3/2565+2197*fk5_4/4104-fk5_5/5; 
fm3_2 = m1 + 16*fk3_1/135+6656*fk3_3/12825+28561*fk3_4/56430-9*fk3_5/50+2*fk3_6/55; 
fm4_2 = m2 + 16*fk4_1/135+6656*fk4_3/12825+28561*fk4_4/56430-9*fk4_5/50+2*fk4_6/55; 
fm5_2 = m3 + 16*fk5_1/135+6656*fk5_3/12825+28561*fk5_4/56430-9*fk5_5/50+2*fk5_6/55; 
R3 = abs(fm3_1-fm3_2)/deltat; 
R4 = abs(fm4_1-fm4_2)/deltat; 
R5 = abs(fm5_1-fm5_2)/deltat; 
err3 = 0.84*(epsilon/R3)^(1/4); 
err4 = 0.84*(epsilon/R4)^(1/4); 
err5 = 0.84*(epsilon/R5)^(1/4); 


if R3<= epsilon 
    time = time+deltat; 
    fm3 = fm3_1; 
    i = i+1; 
    deltat = err3*deltat; 
end 

if R4<= epsilon 
    time = time+deltat; 
    fm4 = fm4_1; 
    i = i+1; 
    deltat = err4*deltat; 
end 

if R5<= epsilon 
    time = time+deltat; 
    fm5 = fm5_1; 
    i = i+1; 
    deltat = err5*deltat; 
end 

e=2*gm^2.0/(2*angmom2); 
ecc=(1.0+(2.0*e*angmom2)/(gm^2.0))^0.5; 
x(i)=r_1*cos(theta_1)/(1000.0*parsec); 
y(i)=r_1*sin(theta_1)/(1000.0*parsec); 
time=time+deltat; 
r(i)=r_1; 
time1(i)=time; 
end 

figure() 
plot(x,y, '-k'); 
TeXString = title('Plot of Orbit in Gamma Potential Obtained Using  RK4') 
xlabel('x') 
ylabel('y')

出典

2016-09-16 A Bit Too Curious