sortie de k3_1 est plafonnée à + 24 -3.1445e

Question

Je résoudre un système d'équations différentielles à l'aide RK4. Je générer un tracé de ligne droite qui semble être dû au fait que k3_1 est plafonnée à + 24 -3.1445e. Je ne comprends pas pourquoi il est plafonné.

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 

Notez que l'anomale se produit sur la ligne indiquée.

Answer 1

Il n'y a pas k3_1 qui est couvert, il est le calcul de u_1 qui retourne une valeur de -3.1445e+24/deltat (deltat est constant).

u_1 est calculé sur la ligne:

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

Après la première itération, retourne:

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 

Sur la base de l'équation u_1(n+1) = u_1(n) + du il ressemble du représente une différence relativement faible. La différence entre les deux premières valeurs est très grande, donc je suppose qu'il est ce calcul est incorrect.

Si vous trouvez que ce calcul est correct, votre erreur est dans une de ces lignes:

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)