2016-08-25 5 views
0

특정 잠재력을 가진 은하 궤도의 플롯을 생성하려고합니다. 내 코드는MATLAB - NaN 오류 받기

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"를 출력하고 있다는 것을 깨달았습니다. 나는 왜 내가 이것을 얻고 있는지 알 수 없다. sqrt()를()^0.5로 대체하려고 시도했지만 여전히 NaN을 산출합니다.

+3

나는 그것이 있다고 생각합니다.^찾고 계십니까? – GameOfThrows

+0

어디에서'NaN'을 생성합니까? 일반적으로 0으로 나누는 경우에만 발생합니다. 언급 된'. ^'@GameOfThrows는 요소의 현명한 힘을 사용합니다. 단순한'^'는 완전한 매트릭스 파워를 취할 것입니다. – Adriaan

+0

@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을 생성했습니다. 코드에서 볼 수 있듯이, 나는 이것을 포함하지 않았다; 출력 내용을 볼 수 있도록 –

답변

0

다소 큰 숫자로 운영되고 있습니다. 가장 큰 부동 소수점 숫자는 1.7e + 308과 같습니다. 저장 장치 번호는 Inf로 평가됩니다. 마찬가지로, 1.7e-308보다 작은 양수는 0으로 평가되고 나눗셈 중에 NaN이 발생합니까?

+0

나는 이것이 사실이라고 생각하지 않는다. 비슷한 일을했는데 코드가 완벽하게 작동했습니다. –