【发布时间】：2016-08-25 09:21:12
【问题描述】：
我正在尝试生成具有特定潜力的银河轨道图。我的代码由
给出 function Eulersystem_MNmodel()
parsec = 3.08*10^18;
r_1 = 8.5*1000.0*parsec; % This converts our value of r into cm.
z_1 = 1.0;
theta_1 = 0.0; %Initial Value for Theta.
U_1 = 100.0*10^5; %Initial value for U in cm/sec
V = 156.972*10^5; %Proposed value of the radial velocity in cm/sec
W_1 = 150*10^5.0;
grav = 6.6720*10^-8; %Universal gravitational constant
amsun = 1.989*10^33; %Proposed Angular momentum of the sun
amg = 1.5d11*amsun;
gm = grav*amg; %Note this is constant
nsteps = 50000; %The number of steps
deltat = 5.0*10^11; %Measured in seconds
angmom = r_1*V; %The angular momentum
angmom2 = angmom^2.0; %The square of the angular momentum
E = -gm/r_1 + U_1*U_1/2 + angmom2/(2.0*r_1*r_1); %Total Energy of the system
time = 0.0;
for i=1:nsteps
r_1 = r_1 + deltat*U_1;
U_1 = U_1 + deltat*((-gm*r_1)/((r_1^2.0 + (1+sqrt(z_1^2.0+1))^2.0)^1.50))
z_1 = z_1 + deltat*W_1;
W_1 = W_1 + deltat*(gm*z_1*(1+sqrt(z_1^2.0+1))/(sqrt(z_1^2.0+1))*(r_1^2.0+(1+sqrt(z_1^2.0+1))^2.0)^1.5);
E = -gm/r_1+U_1/2.0+angmom2/(2.0*r_1*r_1);
ecc = (1.0 + (2.0*E*angmom2)/(gm^2.0))^0.5;
time1(i) = time;
time = time+deltat;
r(i) = r_1;
z(i) = z_1;
end
figure()
plot(r,z)
我不断得到直线，而不是更有趣的曲线，这导致我研究 U_1 函数的输出。经过这次调查，我意识到它一直在输出 NaN “Not a Number”。我不明白为什么我会得到这个。我尝试用 ()^0.5 替换 sqrt()，这仍然产生 NaN。
【问题讨论】：
-
我想是。^你在找吗？
-
它在哪里生成
NaN？通常，这只发生在除以零时。提到的.^@GameOfThrows 正在采取元素明智的力量。简单的^将占用完整的矩阵功率。 -
@Adriaan 行 U_1 = U_1 + deltat*((-gm*r_1)/((r_1^2.0 + (1+sqrt(z_1^2.0+1))^2.0)^1.50) ) 生成 NaN。正如您在代码中看到的，我没有包含 ;这样我就可以看到它在输出什么
-
将该行替换为
U_1 = U_1 + deltat.*((-gm.*r_1)/((r_1.^2.0 + (1+sqrt(z_1.^2.0+1)).^2.0).^1.50))。和其他人一样 -
@AnderBiguri 如果我这样做，我仍然会得到 NaN。
标签： matlab physics numerical-methods astronomy