Problém dvou těles

numerické řešení pohybu dvou těles, které na sebe působí gravitační silou

Zvolíme takové jednotky aby $G = 1$

První těleso má hmotnost $M_1$, druhé hmotnost $M_1$

import numpy as np
import matplotlib.pyplot as plt

G = 1 #gravitační konstanta
M1 = 1 #hmotnost prvního tělesa
M2 = 2 #hmotnost druhého tělesa
mu=1/(1/M1+1/M2) #redukovaná hmotnost

T=6 #doba, po kterou pohyb těles chceme sledovat
dt=0.00001 #časový interval pro numerické řešení
t=np.arange(0,T,dt) #pole časů, ve ketrých budeme počítat polohy těles

souřadnice prvního a druhého tělesa

#souřadnice prvního tělesa
x1=np.empty(np.size(t))
y1=np.empty(np.size(t))
#souřadnice druhého tělesa
x2=np.empty(np.size(t))
y2=np.empty(np.size(t))

relativní poloha obou těles $\vec{r}_{rel} = \vec{r}_2-\vec{r}_1$

xr=np.empty(np.size(t))
yr=np.empty(np.size(t))

relativní rychlost obou těles $\vec{v}_{rel} = \vec{v}_2-\vec{v}_1$

vxr=np.empty(np.size(t))
vyr=np.empty(np.size(t))

Úlohu budeme řešit v těžišťové soustavě, tj. počátek soustavy souřadnic dáme do hmotného středu soustavy těchto dvou hmotných bodů.
V těžišťové soustavě pro polohové vektory hmotných bodů 1 a 2 platí $$M_1\vec{r}_1+M_2\vec{r}_2=0$$ Budeme řešit jednočásticový problám pro fiktivní hmotný bod o hmotnosti $\mu$, polohovém vektoru $\vec{r}_{rel}$ a rychlosti $\vec{v}_{rel}$
Pohybová rovnice, pro takový bod je $$\mu\ddot{\vec{r}}_{rel}=\vec{F}_{1,2},$$ kde $\vec{F}_{1,2}$ je síla, kterou ne sebe obě tělesa působí, v našem případě gravitační síla $$\vec{F}_{1,2} = -G\frac{M_1 M_2}{r_{rel}^2}\frac{\vec{r}_{rel}}{r_{rel}}.$$ Protože pohyb probíhá v rovině, máme po rozepsání do složek následující soustavu diferenciálních rovnic \begin{align} \ddot{x} &= -G\frac{M_1 M_2}{\mu\, r_{rel}^3}x_{rel}\\ \ddot{y} &= -G\frac{M_1 M_2}{\mu\, r_{rel}^3}y_{rel}\\ \end{align} Tuto soustavu diferenciálních rovnic vyřešíme numericky.

počáteční podmínky relativní vzdálenosti a rychlosti obou těles

xr[0]=2
yr[0]=0
vxr[0]=0
vyr[0]=0.5

numerické řešení pohybové rovnice

for i in range(np.size(t)-1):
    xr[i+1]=xr[i]+vxr[i]*dt
    yr[i+1]=yr[i]+vyr[i]*dt
    r=np.sqrt(xr[i+1]**2+yr[i+1]**2)
    vxr[i+1]=vxr[i]-G*M1*M2/mu/r**3*xr[i+1]*dt
    vyr[i+1]=vyr[i]-G*M1*M2/mu/r**3*yr[i+1]*dt

vykreslení trajektorie fiktivního bodu o hmotnosti $\mu$ a polohovém vektoru $\vec{r}_{rel}$

fig,ax=plt.subplots(figsize=(5,5))
ax.plot(xr,yr)
plt.xlim(-1.1,2.1)
plt.ylim(-1.6,1.6)
ax.set_xlabel('$x_{rel}$',fontsize=12)
ax.set_ylabel('$y_{rel}$',fontsize=12)
plt.show()

výpočet vlastních souřadnic obou těles \begin{align} \vec{r}_1 &= -\frac{M_2}{M_1+M_2}\vec{r}_{rel}\\ \vec{r}_2 &= -\frac{M_1}{M_1+M_2}\vec{r}_{rel}\\ \end{align}

x1=-xr*M2/(M1+M2)
x2= xr*M1/(M1+M2)   
y1=-yr*M2/(M1+M2)
y2= yr*M1/(M1+M2)  

Vykreslení trajektorií obou těles

fig,ax=plt.subplots(figsize=(5,5))
ax.plot(x1,y1,c='blue',label="těleso 1")
ax.plot(x2,y2,c='green',label="těleso 2")
plt.xlim(-1.5,1.1)
plt.ylim(-1.3,1.3)
ax.set_xlabel('x',fontsize=12)
ax.set_ylabel('y',fontsize=12)
plt.legend()
plt.show()