Keplerova úloha - pohyb v centrálním poli

Vložení knihoven

import numpy as np # efektivně implementovaná numerická pole
from scipy.optimize import newton # Newtonova metoda hledání kořene transcendentní rovnice
import matplotlib.pyplot as plt # grafy
import matplotlib.animation as animation # animace

Uvažujeme pohyb v rovině $xy$, centrum umístěme do počátku soustavy souřadnic. Označme $\ell$ velikost momentu hybnosti (nenulová je pouze $z$-tová složka). Dále označme $E$ celkovou energii částice, $$E = \frac{1}{2} m \left(\dot{r}^2+r^2\dot{\phi}^2\right)-\frac{\alpha}{r}.$$ Potom s označením $$p = \frac{\ell^2}{m\alpha}, \quad e = \sqrt{1+\frac{2E \ell^2}{m\alpha^2}}$$ dostaneme rovnici dráhy $$\frac{p}{r} = 1+e\cos \phi$$ a dále parametrickou rovnici trajektorie $$\frac{r}{a} = 1-e\cos\xi, \quad \frac{t}{\tau} = \xi-e\sin\xi,$$ kde $$a = \frac{p}{e^2-1}, \quad \tau = \sqrt{\frac{ma^3}{\alpha}}.$$ V kartézských souřadnicích dostaneme $$\frac{x}{a} = \cos\xi - e, \quad \frac{y}{a} = \sqrt{1-e^2} \sin\xi, \quad \frac{t}{\tau} = \xi-e\sin\xi.$$ Vzdálenosti budeme měřit v jednotkách $a$ a časy v jednotkách $\tau$. Zadat stačí pouze jeden parametr a to excentricitu $e$ elipsy.

e = 0.6

Rozdělíme časový interval $[0,\pi]$ na $N$ podintervalů.

N = 100
t = np.linspace(0,np.pi,N)

Dále nadefinujeme funkci jejíž kořen budeme hledat, tj. třetí z rovnic trajektorie, její derivaci a Newtonovou metodou provedeme její inverzi. Dále provedeme vektorizaci funkce, tj. aplikujeme na pole t a získáme pole xi.

def koren_kepler(t):
    def kepler_t(x):
        return x - e*np.sin(x)-t
    
    def d_kepler_t(x):
        return 1 -e*np.cos(x)
    
    return newton(kepler_t,t,fprime=d_kepler_t)

vek_koren_kepler = np.vectorize(koren_kepler)

xi = vek_koren_kepler(t)

Dále nadefinujeme pomocí parametru $\xi$ $x$-ovou a $y$-ovou souřadnici, rovněž rovnou pomocí polí.

f = (1-e**2)**(1/2)
xp = np.cos(xi) + e
yp = f*np.sin(xi)
x = np.concatenate((xp,np.flipud(xp)))
y = np.concatenate((yp,-np.flipud(yp)))

dt = np.pi/N # časový krok

ax = plt.gca()
%matplotlib inline 
from IPython.display import HTML # pro zobrazení animace v HTML

fig = plt.figure() 

ax.grid(False)
ax = fig.add_subplot(111, autoscale_on=False,aspect='equal', xlim=(-1+e-0.1, 1+e+0.1), ylim=(-f-0.1, f+0.1))

cara, = ax.plot([], [], 'o-', lw=2)
cas_template = 'Čas = %.1f'
cas_text = ax.text(0.05, 0.9, '', transform=ax.transAxes)


def init():
    cara.set_data([], [])
    cas_text.set_text('')
    return cara, cas_text


def animuj(i):
    thisx = [0, x[i]]
    thisy = [0, y[i]]

    cara.set_data(thisx, thisy)
    cas_text.set_text(cas_template % (i*dt))
    return cara, cas_text


ani = animation.FuncAnimation(fig, animuj, np.arange(1, len(y)),
                              interval=25, blit=True, init_func=init)
plt.close()

ani.save('Keplerova_uloha.mp4', fps=15) # animaci uložíme jako video

HTML(ani.to_html5_video())