Błędy podczas rozwiązywania python ODE

Question

Mam projekt uniwersytecki, w którym jesteśmy proszeni o symulację satelitarnego podejścia do Marsa za pomocą ODE i funkcji odint SciPy.

Udaje mi się symulować to w 2D, tworząc ODE drugiej rzędu w dwóch ODE pierwszego rzędu. Jednak utknąłem w ograniczeniu czasowym, ponieważ mój kod używa jednostek SI, dlatego działa w ciągu kilku sekund, a granice obszaru Pythona nie symulują nawet jednej pełnej orbity.

Próbowałem konwertować zmienne i stałe na godziny i kilometry, ale teraz kod nadal podaje błędy.

Śledziłem tę metodę:

http://bulldog2.redlands.edu/facultyfolder/deweerd/tutorials/Tutorial-ODEs.pdf

I kod to:

import numpy 

import scipy 

from scipy.integrate import odeint 

def deriv_x(x,t): 
    return array([ x[1], -55.3E10/(x[0])**2 ]) #55.3E10 is the value for G*M in km and hours 

xinit = array([0,5251]) # this is the velocity for an orbit of period 24 hours 

t=linspace(0,24.0,100) 

x=odeint(deriv_x, xinit, t) 

def deriv_y(y,t): 
    return array([ y[1], -55.3E10/(y[0])**2 ]) 

yinit = array([20056,0]) # this is the radius for an orbit of period 24 hours 

t=linspace(0,24.0,100) 

y=odeint(deriv_y, yinit, t) 

nie wiem jak skopiować/wkleić kod błędu z PyLab więc wziąłem PrintScreen błędu:

Drugi błąd przy t = linspace (0.01,24.0,100) i xinit = array ([0.001,5251]):

Jeśli ktoś ma jakieś sugestie, jak poprawić kod będę bardzo wdzięczny.

Dziękuję bardzo!

Answer 1

odeint(deriv_x, xinit, t) 

wykorzystuje xinit jako wstępnego odgadnięcia x. Ta wartość dla x jest używana przy ocenie deriv_x.

deriv_x(xinit, t) 

zgłasza błąd dzielenia przez zero, ponieważ x[0] = xinit[0] równa 0, a deriv_x podziały według x[0].

---

Wygląda na to, że starają się rozwiązać Ody drugiego rzędu

r'' = - C rhat 
     --------- 
     |r|**2 

gdzie rhat jest wektorem jednostkowym w kierunku promieniowym.

Wygląda do oddzielania współrzędnych x i y w oddzielnych równań różniczkowych drugiego rzędu:

x'' = - C    y'' = - C 
     ----- and   ----- 
     x**2     y**2 

początkowych warunkach X0 = 0 i = y0 20056.

ten jest bardzo problematyczne. Wśród problemów jest to, że gdy x0 = 0, x'' wysadzi w powietrze. Pierwotna ODE dla drugiego rzędu dla r'' nie ma tego problemu - mianownik nie wysadza się, gdy x0 = 0 ponieważ y0 = 20056, a więc jest daleki od zera.

Wniosek: Twoja metoda rozdzielania ODE r'' na dwie jednostki ODE dla i y'' jest niepoprawna.

Spróbuj wyszukać inny sposób rozwiązania r'' ODE.

Podpowiedź:

• Co jeśli "stan" wektor jest z = [x, y, x', y']?
• można spisać ODE pierwszego rzędu dla z' pod względem x, y, x' i y'?
• Czy możesz go rozwiązać, dzwoniąc pod numer do integrate.odeint?