2014-02-05 1 views
0

나는 파이썬에 대해 매우 익숙하며 다음 코드를 생성했지만 작동하지 않는다.화성에 접근중인 위성의 움직임을 시뮬레이트

이 코드는 화성에 접근하는 위성의 움직임을 변수 yx 거리 (즉 2 차원 공간)에 매핑하려고 시도하고 있습니다. 현재 화성은 정지 된 것으로 간주됩니다.

import scipy as sp 
import numpy as np 
import pylab as pl 
import scipy.integrate as spi 

G=6.67*(10**-11) 
mm=6.4*(10**23) 
#^mm is the mass of mars and G is the gravitational constant  

def f(b,t): 


    xx=b[0] 
    vx=b[1] 
    yy=b[2] 
    vy=b[3] 

    ax=-(G*mm*b[0])/((b[0]**2) + (b[2]**2))**1.5 
    ay=-(G*mm*b[2])/((b[0]**2) + (b[2]**2))**1.5 


    return [vx,ax,vy,ay] 

t=sp.linspace(0.,10000000.,1000) 

xx0=[800000., 0., 10000, 0] 

soln=spi.odeint(f,xx0,t) 

print soln 

x=soln[:,0] 
v1=soln[:,1] 
y=soln[:,2] 
v2=soln[:,3] 

pl.figure(1) 
pl.plot(x,y) 
pl.xlabel("x displacement") 
pl.ylabel("y displacement") 

pl.show() 

이것은 직선 그래프를 그립니다 ... 아마도 물리학을 오해하지만 코드에 문제가 있습니까?

FYI - 뉴턴의 중력 법칙과 일부 기하학 법칙을 사용하여 가속 방정식을 계산했습니다. 저는 주로 코드의 문제에 대해 질문하고 있지만, 물리학 자의 물리학 문제가 있으면 물리 칠 문제가 있다는 것을 듣고 싶습니다.

+3

나는 scipy이 없어? 즉, 궤도를 그릴 때 직선이어야합니다. 초기 속도를주지 않아 직선으로 떨어지게해야합니다. 'plot (t, y)'를 시도해 가속하고 있는지 확인하고 초기 'x'속도를 시도해보십시오. – tom10

+2

함수 안에 t가 누락되지 않았습니까? – Martin

+0

환호성 그렇습니다. 저는 이전에 너무 많은 숫자가 있다는 것을 깨달았습니다. 올바른 감사입니다. – tcrules

답변

3

방정식을 확인하고 올바른지 (자신을 의심할만한 이유가있는 것은 아닙니다). :)

일이 잘못되었다고 생각하는 이유는 처음에는 G * mm~10**10과 같은 많은 양의 초기 조건에서 매우주의해야합니다. 초기 조건의 크기 차이는 위성 궤도 화성과 태양계 바깥에있는 하나. 둘째 방정식이 (x=0, y=0에 가까워지면 합리적인 해결책이 발견 될지라도 (즉 초기 조건은 맞을 것입니다.) 사물은 단수로 나타나고 위성으로 이어지는 큰 가속도로 결국 지역 규모가 사라집니다.

이 문제를 해결하려면 x, y가 0에 가깝지 않고 적절한 경우 ode 해석기를 끝내는 적절한 검사를 할 수 있습니다.

두 번째로 이런 종류의 문제를 해결할 때는 nondimensionalisation을 사용하는 것이 현명합니다. 이것은 종속 변수 (x->x', y->y')의 크기를 다시 줄여서 준비된 버전이 순서가 일치하도록하는 방법입니다. 귀하의 경우에는

단순히 좌표 시간을 무시합니다 및 X = X '* L Y = y'를 보자 * L은 특성 길이 규모이다 L

. 방정식에서 이것을 대입하면 G * mm/L**3이라는 용어로 끝납니다. 이 수량이 1이되도록 특성 길이를 설정하십시오 : L = (G * mm)**(1/3.0).

이렇게하면 속도가 증가하는 숫자를 쉽게 처리 할 수 ​​있다는 장점이 있습니다. 이제는 단지 xx = 10km 및 xx '= xx/L 등의 적절한 값을 취하여 적절한 초기 조건을 설정하는 것만이 어려워집니다.

이 방법을 사용하면 멋진 궤도를 수행 할 수 있습니다. 내가 코드를 실행할 수 없습니다하지만 ... 당신이 작동하지 않습니다 확신 있도록

enter image description here

+0

L을 설정하면 무슨 뜻인지 자세히 설명해 주실 수 있습니까? 예를 들어 파워 세 번째와 같이 꽤 잘 할 수 없습니다. – tcrules

+0

죄송합니다.이 글을 쓰고 나 죄송합니다. [http] : //imgur.com/P0NJuAe) 나는 이것을 공유하는 가장 쉬운 방법이라고 생각했다. 나는 이전에 게으르다. 당신은 L 단위를 계산하고 (t = t '기억) 이 문제에 대한 일반적인 기술입니다. – Greg

+0

감사합니다. 그러나 코드에 나타났습니다. 필자가 작성한 것이 아니라 [ax, vx, ay, vy]를 반환해야합니다. 조금 바뀌는 – tcrules

관련 문제