Drielichamen.ipynb

Kernel: Python 3 (Ubuntu Linux)

Het drielichamenprobleem

We kunnen de positie $\mathbf{x}_i \in \mathbb{R}^3$ en de impuls $\mathbf{p}_i \in \mathbb{R}^3$ van 3 lichamen van gelijke massa onder invloed van elkaars zwaartekracht beschijven met een stelsel differentiaalvergelijkingen

$$\dot{\mathbf{x}}_i(t) = \mathbf{p}_i(t),$$$$\dot{\mathbf{p}}_i(t) = \sum_{i\not=j} \frac{-1}{\|\mathbf{x}_i - \mathbf{x}_j\|_2^3}(\mathbf{x}_i - \mathbf{x}_j),$$

voor $i=1,2,3$.

We kunnen dit schrijven als

$$\dot{\mathbf{x}}(t) = F_1(\mathbf{x}(t),\mathbf{p}(t)),$$$$\dot{\mathbf{p}}(t) = F_2(\mathbf{x}(t),\mathbf{p}(t)),$$

waar $\mathbf{x}(t) \in \mathbb{R}^9$ de posities van alle lichamen bevat, etc. Dit geeft uiteindelijk een stelsel van 18 differentiaalvergelijkingen.

Er zijn een aantal stabiele banen bekend, één daarvan heeft beginposities

$$\mathbf{x}_1(0) = (-1,0,0)^T, \mathbf{x}_2(0) = (1,0,0)^T, \mathbf{x}_3(0) = (0,0,0)^T,$$

en beginimpuls

$$\mathbf{p}_1(0) = (a,b,0)^T, \mathbf{p}_2(0) = (a,b,0)^T, \mathbf{p}_3(0) = -2(a,b,0)^T,$$

met $a = 0.347111$ en $b = 0.532728$. De periode van de baan is $T = 6.324449$.

De totale energy

$$E(t) = \frac{1}{2} \sum_{i=1}^3 \|\mathbf{p}_i(t)\|_2^2 - \sum_{i=1}^3\sum_{j = 1}^{i-1} \frac{1}{\|\mathbf{x}_i - \mathbf{x}_j\|_2},$$

is behouden.

Vragen voor het verslag

Basis

• Probeer of je met Euler Forward de bovenstaande stabiele baan kunt berekenen.

• Onderzoek de stabiliteit van EF en verklaar de resultaten. Hint: bereken de Jacobiaan en kijk naar de eigenwaarden op het gevonden traject.

• In Matlab is ook een adaptieve methode beschikbaar (ode23) die zelf een geschikte tijdstap kiest om de relatieve fout binnen een bepaalde marge te houden. Probeer of het met deze methode lukt om en gesloten baan te produceren.

• Kijk ook eens naar het behoud van energie.

Extra

• Een uitgebreid overzicht van stabiele banen is hier te vinden. Probeer of je ook andere stabiele configuraties kunt simuleren.

• Populair voor dit soort systemen zijn zogenaamde symplectische methoden. De eenvoudigste is symplectische Euler: $$\mathbf{x}^{n+1} = \mathbf{x}^n + \Delta t F_1(\mathbf{x}^n,\mathbf{p}^n),$$ $$\mathbf{p}^{n+1} = \mathbf{p}^n + \Delta t F_2(\mathbf{x}^{n+1},\mathbf{p}^n).$$ Lukt het beter om een gesloten baan de simuleren? En hoe zit het met behoud van energie?

• Pas Euler Backward op de vergelijking toe, werkt deze beter dan de symplectische methode?

Code

Om je opweg te helpen kun je voortbouwen op onderstaand voorbeeld

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# definieer stelsel vergelijkingen
F = lambda u,t : np.array([u[9:12],u[12:15],u[15:18],- ((u[0:3]-u[3:6])/np.linalg.norm(u[0:3]-u[3:6])**3 + (u[0:3]-u[6:9])/np.linalg.norm(u[0:3]-u[6:9])**3),- ((u[3:6]-u[0:3])/np.linalg.norm(u[3:6]-u[0:3])**3 + (u[3:6]-u[6:9])/np.linalg.norm(u[3:6]-u[6:9])**3),- ((u[6:9]-u[0:3])/np.linalg.norm(u[6:9]-u[0:3])**3 + (u[6:9]-u[3:6])/np.linalg.norm(u[6:9]-u[3:6])**3)]).flatten()


# definieer de Jacobiaan (dit is een dummy, dus zelf aanpassen)
DF = lambda u,t : np.eye(18)

# beginwaarden
a = 0.347111
b = 0.532728

x10 = np.array([-1,0,0])
x20 = np.array([1,0,0])
x30 = np.array([0,0,0])
p10 = np.array([a,b,0])
p20 = np.array([a,b,0])
p30 = np.array([-2*a,-2*b,0])

u0 = np.array([x10,x20,x30,p10,p20,p30]).flatten()

# eindtijd en tijdstap
dt = 1e-3
T  = 10
t  = np.linspace(0,T,int(T/dt)+1)

# los op met odeint ('exacte' oplossing)
u1 = odeint(F,u0,t)

# los op met Forward Euler (benadering)
u2 = np.zeros((len(t),18))
u2[0] = u0
for k in range(len(t)-1):
    u2[k+1] = u2[k] + dt*F(u2[k],t)

# reken total energie uit
E1 = 0.5*np.linalg.norm(u1[:,9:18],axis=1) - 1/np.linalg.norm(u1[:,0:3]-u1[:,3:6],axis=1) - 1/np.linalg.norm(u1[:,0:3]-u1[:,6:9],axis=1) - 1/np.linalg.norm(u1[:,3:6]-u1[:,6:9],axis=1)
E2 = 0.5*np.linalg.norm(u2[:,9:18],axis=1) - 1/np.linalg.norm(u2[:,0:3]-u2[:,3:6],axis=1) - 1/np.linalg.norm(u2[:,0:3]-u2[:,6:9],axis=1) - 1/np.linalg.norm(u2[:,3:6]-u2[:,6:9],axis=1)

# berekend de eigenwaarden langs het traject
mu1 = np.zeros((len(t),18))
mu2 = np.zeros((len(t),18))
for k in range(len(t)):
    mu1[k] = np.linalg.eigvals(DF(u1[k],t[k]))
    mu2[k] = np.linalg.eigvals(DF(u2[k],t[k]))
    
# plot oplossing
fig1, axs1 = plt.subplots(3,1)

axs1[0].set_title('Baan van x_1')
axs1[0].plot(u1[:,0],u1[:,1],'k--')
axs1[0].plot(u2[:,0],u2[:,1],'r')
axs1[0].set_aspect('equal','box')

axs1[1].set_title('Baan van x_2')
axs1[1].plot(u1[:,3],u1[:,4],'k--')
axs1[1].plot(u2[:,3],u2[:,4],'r')
axs1[1].set_aspect('equal','box')

axs1[2].set_title('Baan van x_3')
axs1[2].plot(u1[:,6],u1[:,7],'k--')
axs1[2].plot(u2[:,6],u2[:,7],'r')
axs1[2].set_aspect('equal','box')

# plot totale energie
fig2, axs2 = plt.subplots(1,1)

axs2.set_title('totale energie')
axs2.plot(t,E1,'k--',t,E2,'r')
axs2.legend(('odeint','forward euler'))

# plot eigenwaarden langs traject
fig3, axs3 = plt.subplots(1,1)
axs3.set_title('eigenwaarden')
axs3.plot(t,mu1,'k--',t,mu2)

[<matplotlib.lines.Line2D at 0x7f91d0a379b0>,
 <matplotlib.lines.Line2D at 0x7f91d0a37b00>,
 <matplotlib.lines.Line2D at 0x7f91d0a37c50>,
 <matplotlib.lines.Line2D at 0x7f91d0a37da0>,
 <matplotlib.lines.Line2D at 0x7f91d0a37ef0>,
 <matplotlib.lines.Line2D at 0x7f91d0b77080>,
 <matplotlib.lines.Line2D at 0x7f91d0b771d0>,
 <matplotlib.lines.Line2D at 0x7f91d0b77320>,
 <matplotlib.lines.Line2D at 0x7f91d0b77470>,
 <matplotlib.lines.Line2D at 0x7f91d0b775c0>,
 <matplotlib.lines.Line2D at 0x7f91d0b77710>,
 <matplotlib.lines.Line2D at 0x7f91d0b77860>,
 <matplotlib.lines.Line2D at 0x7f91d0b779b0>,
 <matplotlib.lines.Line2D at 0x7f91d0b77b00>,
 <matplotlib.lines.Line2D at 0x7f91d0b77c50>,
 <matplotlib.lines.Line2D at 0x7f91d0b77da0>,
 <matplotlib.lines.Line2D at 0x7f91d0b77ef0>,
 <matplotlib.lines.Line2D at 0x7f91d0b6b080>,
 <matplotlib.lines.Line2D at 0x7f91d0b89a20>,
 <matplotlib.lines.Line2D at 0x7f91d0b89ba8>,
 <matplotlib.lines.Line2D at 0x7f91d0b89cf8>,
 <matplotlib.lines.Line2D at 0x7f91d0b89e48>,
 <matplotlib.lines.Line2D at 0x7f91d0b89f98>,
 <matplotlib.lines.Line2D at 0x7f91d0b6e128>,
 <matplotlib.lines.Line2D at 0x7f91d0b6e278>,
 <matplotlib.lines.Line2D at 0x7f91d0b6e3c8>,
 <matplotlib.lines.Line2D at 0x7f91d0b6e518>,
 <matplotlib.lines.Line2D at 0x7f91d0b6e668>,
 <matplotlib.lines.Line2D at 0x7f91d0b6e7b8>,
 <matplotlib.lines.Line2D at 0x7f91d0b6e908>,
 <matplotlib.lines.Line2D at 0x7f91d0b6ea58>,
 <matplotlib.lines.Line2D at 0x7f91d0b6eba8>,
 <matplotlib.lines.Line2D at 0x7f91d0b6ecf8>,
 <matplotlib.lines.Line2D at 0x7f91d0b6ee48>,
 <matplotlib.lines.Line2D at 0x7f91d0b6ef98>,
 <matplotlib.lines.Line2D at 0x7f91d0b59128>]