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 有上限，而是u_1 的计算返回值-3.1445e+24 / deltat（deltat 是常数）。

u_1在行中计算：

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

第一次迭代后，返回：

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

基于等式u_1(n+1) = u_1(n) + du，看起来du 代表了一个相对较小的差异。前两个值之间的差异非常很大，所以我假设这个计算不正确。

如果您发现该计算是正确的，那么您的错误出现在以下几行之一：

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)