MATLAB - 自适应步长Runge-Kutta

我在MATLAB中编写了一个自适应步长RK4来解决一个ODE系统。代码运行时没有错误，但是当我尝试对y绘制x时，它不会生成所需的曲线。而不是一个环形的形状，我只是得到一条平坦的线。这从r正在输出一个常数值的事实可以看出。在检查每行的输出后，它们不输出常量或错误或inf或NaN，而是输出实数和虚数分量（复数）。我不知道为什么会发生这种情况，我相信这是我麻烦的根源。MATLAB - 自适应步长Runge-Kutta

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