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의 값을 (deltat이 일정) 반환 u_1의 계산이다.

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 

는 상대적으로 작은 차이를 나타낸다. 두 가지 첫 번째 값의 차이는 입니다. 매우입니다. 따라서이 계산이 잘못되었다고 가정합니다. 당신이 계산이 정확하다는 것을 발견하면

, 다음 오류는이 라인 중 하나에 :

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)