matlab - MATLAB-自适应步长Runge-Kutta

Question

我在MATLAB中编程了一个自适应步长RK4来求解ODE系统。代码运行时没有错误，但是当我尝试将x与y作图时，它不会产生所需的曲线。我不再是环形形状，而是获得了一条扁平线。从r输出恒定值这一事实可以明显看出这一点。在检查每行的输出之后，它们没有输出常量或错误或inf或NaN，而是输出了实部和虚部(复数)。我不知道为什么会这样，我相信这是我麻烦的根源。

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')

Answer 1

您将获得复杂的值，因为有时是npts - time < 0。您可能需要打印出deltat的值以检查错误。

另外，当错误估计值大于您的容忍度时，您的代码似乎没有考虑到这种情况。当您的误差估计值大于公差时，您必须:

移回时间和解决方案

根据公式计算新的步长，然后

重新计算您的解决方案和错误估计。

您不知道必须进行多少次迭代的事实使for循环适用于自适应梯级Kutta有点尴尬。我建议改用while循环。