Startseite
| Programme:
Dreikörperproblem
| Algorithmen,
Lektionen, Analysen
|
Unsere Vorläufer
| Kontakt, Datenschutz
Startseite
| Programme:
Dreikörperproblem
| Algorithmen,
Lektionen, Analysen
|
Unsere Vorläufer
| Kontakt, Datenschutz
Schnell-Anwahl:
VdS-Journal 81 – Lage der
Lagrange-Punkte 1-3
VdS-Journal 82 – Stabilität der
Lagrangepunkte
VdS-Journal 83 – Periodische Bahnen
# Lagrange1-3 : Berechnung der Orte der Lagrange-Punkte
from math import *
def sq(x): return x*x
G = 6.67408e-11
# Sonne - Jupiter
m1=1.988475415966536e+30 # Masse Sonne
m2=1.898568695e+27 # Masse Jupiter
r=7.7835e11 # Abstand für Kreisbahn
x1=-m2*r/(m1+m2) # Koordinate ...
x2=m1*r/(m1+m2) # ...vom Schwerpunkt
omega=sqrt(G*(m1+m2)/r/r/r)
v1=omega*x1 # Bahngeschwindigkeit
v2=omega*x2
print("x1=",x1,"v1=", v1)
print("x2=",x2,"v2=",v2)
def beschl(xL,x1,x2): # effektive Beschleunigung
a1=G*m1/sq(xL-x1)*(x1-xL)/abs(x1-xL)
a2=G*m2/sq(xL-x2)*(x2-xL)/abs(x2-xL)
aw=sq(omega)*xL
a=a1+a2+aw
return a
for L in range(1,4): # alle drei Lagrangepunkte
if (1==L): xL1=0.001*x2;xL2=0.9999999*x2
if (2==L): xL1=1.0001*x2; xL2=2*x2
if (3==L): xL1=-2*x2; xL2=-1.0001*x1
# Berechnung des Lagrangepunktes mit dem Halbierungsverfahren
a1=beschl(xL1,x1,x2)
a2=beschl(xL2,x1,x2)
for i in range(60):
xL=(xL1+xL2)/2
aL=beschl(xL,x1,x2)
if (aL>0): xL2=xL;a2=aL
else: xL1=xL; a1=aL
print("L",L,": xL=",xL, "vL=", omega*xL)
name=input("Fertig?")
So wie abgedruckt berechnet das Programm das Verhalten von Körpern in der Nähe der Lagrangepunkte. Man kann das umstellen auf exakt die Langrangepunkte, wenn man ungenau = false setzt.
# Programm StabLagrange.py - Analyse der Stabilität
from turtle import *
from random import *
from math import *
from plot import *
# globale Variable
G = 6.67408e-11 #Gravitationskonstante
mProAE = 149597870700 #Astronomische Einheit
sekProTag = 86400
true = 1
false = 0
# statische Variable für rkf5
rErrMin = 0
rErrMax = 0
dt = 0
N = 7 # Anzahl der Körper
r0 = [[0 for i in range(7)] for j in range(7)]
v0 = [[0 for i in range(7)] for j in range(7)]
a0 = [[0 for i in range(7)] for j in range(7)]
r1 = [[0 for i in range(7)] for j in range(7)]
v1 = [[0 for i in range(7)] for j in range(7)]
a1 = [[0 for i in range(7)] for j in range(7)]
r2 = [[0 for i in range(7)] for j in range(7)]
v2 = [[0 for i in range(7)] for j in range(7)]
a2 = [[0 for i in range(7)] for j in range(7)]
r3 = [[0 for i in range(7)] for j in range(7)]
v3 = [[0 for i in range(7)] for j in range(7)]
a3 = [[0 for i in range(7)] for j in range(7)]
r4 = [[0 for i in range(7)] for j in range(7)]
v4 = [[0 for i in range(7)] for j in range(7)]
a4 = [[0 for i in range(7)] for j in range(7)]
r5 = [[0 for i in range(7)] for j in range(7)]
v5 = [[0 for i in range(7)] for j in range(7)]
a5 = [[0 for i in range(7)] for j in range(7)]
# r0_4: Ergebnis für RKF 4. Ordnung zur Fehlerberechnung
r0_4 = [[0 for i in range(7)] for j in range(7)]
m = [0 for i in range(7)]
def startwerte():
ungenau = true # kleine Abweichung, um die Stabilität zu prüfen
Y = 0
f = 1.0009548 # Korrekturfaktor für L4L5 bei Jupitermasse
T = 1 # simulierte Zeit in Jupiterumläufen
P = 100 # Anzahl Zwischenwerte
if (ungenau):
Y = 0.001 # bei L1 bis L3 den Ort und die Bahngeschwindigkeit leicht verschieben
f = 1.0 # bei L4 und L5 den Korrekturfaktor weglassen
T = 120 # Hufeisen von L3 erfordert 120 Umläufe
P = 6000 # ... und ordentlich Zwischenwerte
# Massen
m[0] = 1.988475415966536e+30
m[1] = 1.898568695e+27
m[2] = 0
m[3] = 0
m[4] = 0
m[5] = 0
m[6] = 0
# Orte und Geschwindigkeiten
r0[0][0] = -742448884.05 # Sonne
r0[0][1] = 0
r0[0][2] = 0
v0[0][0] = 0
v0[0][1] = -12.461410754
v0[0][2] = 0
r0[1][0] = 777607551115.94 # Jupiter
r0[1][1] = 0
r0[1][2] = 0
v0[1][0] = 0
v0[1][1] = 13051.520863369
v0[1][2] = 0
rJ = r0[1][0] # der Jupiterradius ist konstant in dieser Rechnung
vJ = v0[1][1] # die Bahngeschwindigkeit auch
r0[2][0] = 0.9332557970715 * rJ # L1
r0[2][1] = Y * r0[1][1]
r0[2][2] = 0
v0[2][0] = -Y * v0[1][1]
v0[2][1] = 0.9332557970715 * vJ
v0[2][2] = 0
r0[3][0] = 1.0698510256405 * rJ # L2
r0[3][1] = Y * r0[1][1]
r0[3][2] = 0
v0[3][0] = -Y * v0[1][1]
v0[3][1] = 1.0698510256405 * vJ
v0[3][2] = 0
r0[4][0] = -1.0013526135998 * rJ # L3
r0[4][1] = Y * r0[1][1]
r0[4][2] = 0
v0[4][0] = -Y * v0[1][1]
v0[4][1] = -1.0013526135998 * vJ
v0[4][2] = 0
# L4 und L5
r0[5][0] = f * (0.5 * (r0[0][0] + r0[1][0]) - r0[0][0])
r0[5][1] = f * ((r0[0][0] + r0[1][0]) * sqrt(3) / 2)
r0[5][2] = 0
v0[5][0] = -r0[5][1] * vJ / rJ
v0[5][1] = r0[5][0] * vJ / rJ
v0[5][2] = 0
# L5
r0[6][0] = r0[5][0]
r0[6][1] = -r0[5][1]
r0[6][2] = 0
v0[6][0] = -r0[6][1] * vJ / rJ
v0[6][1] = r0[6][0] * vJ / rJ
v0[6][2] = 0
return rJ, T , P
def sq(x):
return x * x
def abstand(r, p, q):
return sqrt(sq(r[p][0] - r[q][0]) + sq(r[p][1] - r[q][1]) + sq(r[p][2] - r[q][2]))
def beschleunigung(N, m, r, a):
for p in range(N): # alle Beschleunigungen löschen
for k in range(3):
a[p][k] = 0
for p in range(N): # alle Planeten
for q in range(p + 1, N): # gegen alle anderen Planeten
R = abstand(r, p, q)
R3 = R * R * R
for k in range(3):
A = G * (r[q][k] - r[p][k]) / R3
a[p][k]+=m[q] * A
a[q][k]-=m[p] * A
def rkf5(tEnd, init):
global dt, rErrMin, rErrMax
t = 0
while(t < tEnd):
t = t + dt
beschleunigung(N, m, r0, a0)
for i in range(N):
for k in range(3):
r1[i][k] = r0[i][k] + 1 / 4 * dt * v0[i][k]
v1[i][k] = v0[i][k] + 1 / 4 * dt * a0[i][k]
beschleunigung(N, m, r1, a1)
for i in range(N):
for k in range(3):
r2[i][k] = r0[i][k] + 3 / 32 * dt * v0[i][k] + 9 / 32 * dt * v1[i][k]
v2[i][k] = v0[i][k] + 3 / 32 * dt * a0[i][k] + 9 / 32 * dt * a1[i][k]
beschleunigung(N, m, r2, a2)
for i in range(N):
for k in range(3):
r3[i][k] = r0[i][k] + 1932 / 2197 * dt * v0[i][k] + -7200 / 2197 * dt * v1[i][k] + \
7296 / 2197 * dt * v2[i][k]
v3[i][k] = v0[i][k] + 1932 / 2197 * dt * a0[i][k] + -7200 / 2197 * dt * a1[i][k] + \
7296 / 2197 * dt * a2[i][k]
beschleunigung(N, m, r3, a3)
for i in range(N):
for k in range(3):
r4[i][k] = r0[i][k] + 439 / 216 * dt * v0[i][k] - 8 * dt * v1[i][k] + 3680 / 513 * \
dt * v2[i][k] - 845 / 4104 * dt * v3[i][k]
v4[i][k] = v0[i][k] + 439 / 216 * dt * a0[i][k] - 8 * dt * a1[i][k] + 3680 / 513 * \
dt * a2[i][k] - 845 / 4104 * dt * a3[i][k]
beschleunigung(N, m, r4, a4)
for i in range(N):
for k in range(3):
r5[i][k] = r0[i][k] - 8 / 27 * dt * v0[i][k] + 2 * dt * v1[i][k] - 3544 / 2565 * \
dt * v2[i][k] + 1859 / 4104 * dt * v3[i][k] - 11 / 40 * dt * v4[i][k]
v5[i][k] = v0[i][k] - 8 / 27 * dt * a0[i][k] + 2 * dt * a1[i][k] - 3544 / 2565 * \
dt * a2[i][k] + 1859 / 4104 * dt * a3[i][k] - 11 / 40 * dt * a4[i][k]
beschleunigung(N, m, r5, a5)
for i in range(N):
for k in range(3):
r0_4[i][k] = r0[i][k]
# Ergebniss 4. Ordnung zur Fehlerberechnung
r0_4[i][k]+= dt * (25 / 216 * v0[i][k] + 0 * v1[i][k] + 1408 / 2565 * v2[i][k] + \
2197.0 / 4104 * v3[i][k] - 1 / 5 * v4[i][k] + 0 * v5[i][k])
# eigentliches Ergebnis 5. Ordnung
r0[i][k]+= dt * (16 / 135 * v0[i][k] + 0 * v1[i][k] + 6656 / 12825 * v2[i][k] + \
28561 / 56430 * v3[i][k] - 9 / 50 * v4[i][k] + 2 / 55 * v5[i][k])
v0[i][k]+= dt * (16 / 135 * a0[i][k] + 0 * a1[i][k] + 6656 / 12825 * a2[i][k] + \
28561 / 56430 * a3[i][k] - 9 / 50 * a4[i][k] + 2 / 55 * a5[i][k])
# Fehler ausrechnen. Benutzung erfordert eine while- statt einer
# for-Schleife
err = 0
for i in range(N): # welcher Körper gibt den größten Fehler?
localErr = 0
for k in range(3):
localErr+=sq(r0_4[i][k] - r0[i][k])
localErr = sqrt(localErr)
if (localErr > err):
err = localErr
if (init > 0): # Fehler wird in der ersten Runde so gesetzt,
rErrMin = err / 2 # wie es dem ersten dt entspricht
rErrMax = err
init = false
if (err < rErrMin): # Schrittweite verändern, falls der Fehler
dt*=1.2 # außerhalb der Grenzen
elif (err > rErrMax):
dt/=1.2
if (dt / sekProTag < 0.01):
Err = input("Schrittweitensteuerung versagt: dt erhöhen")
exit()
return t
# Hauptprogramm
rJ, T, P = startwerte()
initKoor(-1.3,1.3,0.2, -1.1,1.1,0.2)
dt = 2.5 # 2,5 Tage
# den anderen Teil der Startwerte in SI umwandeln
dt = dt * sekProTag
UJup = 2 * pi * rJ / v0[1][1]
T = T * UJup
# Punkte für Jupiter und Sonne
Plot(r0[0][0] / rJ,r0[0][1] / rJ,20, "gold")
Plot(r0[1][0] / rJ,r0[1][1] / rJ,8,"coral")
t = 0
init = true
while (t < T):
t = t + rkf5(T / P, init)
init = false # nur in der ersten Runde soll die interne Fehlergrenze gesetzt werden.
xJ = r0[1][0]
yJ = r0[1][1]
for i in range(2,7):
x = r0[i][0]
y = r0[i][1]
# in mitrotierene Koordinaten umrechnen (Skalarprodukt)
rP = sqrt(x * x + y * y)
xP = xJ * x + yJ * y
xP = xP / rJ / rJ
yP = -yJ * x + xJ * y
yP = yP / rJ / rJ
if i==2 and t<5*UJup: # L1
Plot(xP,yP,3,"blue")
if i==3 and t<2*UJup: # L2
Plot(xP,yP,3,"red")
if i==4: # L3
Plot(xP,yP,3,"magenta")
if i>4:
Plot(xP,yP,3,"green")
update() # Zwischenergebnis anzeigen
update()
done()
Das Programm rechnet wie die beiden vorigen in mitrtierenden Koordinaten. Den Bewegungsablauf in natürlichen Koordinaten erhält man, wenn man die erste Plot-Zeile auskommentiert und die beiden nächsten einkommentiert.
Als Kommentar sind alle im Heft angegebenen Startwertesätze gegeben. Diese müssen auf die Koordinaten X, Y, U, V geschrieben werden.
# 3koerperPeriodisch.py - periodische Bahnen
from turtle import *
from random import *
from math import *
from plot import *
# globale Variable
G = 6.67408e-11 #Gravitationskonstante
mProAE = 149597870700 #Astronomische Einheit
sekProTag = 86400
true = 1
false = 0
# statische Variable für rkf5
rErrMin = 0
rErrMax = 0
dt = 0
N = 3 # Anzahl der Körper
r0 = [[0 for i in range(N)] for j in range(N)]
v0 = [[0 for i in range(N)] for j in range(N)]
a0 = [[0 for i in range(N)] for j in range(N)]
r1 = [[0 for i in range(N)] for j in range(N)]
v1 = [[0 for i in range(N)] for j in range(N)]
a1 = [[0 for i in range(N)] for j in range(N)]
r2 = [[0 for i in range(N)] for j in range(N)]
v2 = [[0 for i in range(N)] for j in range(N)]
a2 = [[0 for i in range(N)] for j in range(N)]
r3 = [[0 for i in range(N)] for j in range(N)]
v3 = [[0 for i in range(N)] for j in range(N)]
a3 = [[0 for i in range(N)] for j in range(N)]
r4 = [[0 for i in range(N)] for j in range(N)]
v4 = [[0 for i in range(N)] for j in range(N)]
a4 = [[0 for i in range(N)] for j in range(N)]
r5 = [[0 for i in range(N)] for j in range(N)]
v5 = [[0 for i in range(N)] for j in range(N)]
a5 = [[0 for i in range(N)] for j in range(N)]
# r0_4: Ergebnis für RKF 4. Ordnung zur Fehlerberechnung
r0_4 = [[0 for i in range(N)] for j in range(N)]
m = [0 for i in range(N)]
def startwerte():
# Orte und Geschwindigkeiten
r0[0][0] = -742448884.05
r0[0][1] = 0
r0[0][2] = 0
v0[0][0] = 0
v0[0][1] = -12.461410754
v0[0][2] = 0
r0[1][0] = 777607551115.94
r0[1][1] = 0
r0[1][2] = 0
v0[1][0] = 0
v0[1][1] = 13051.520863369
v0[1][2] = 0
r0[2][0] = 0 # wird vom Hauptprogramm gesetzt
r0[2][1] = 0
r0[2][2] = 0
v0[2][0] = 0
v0[2][1] = 0 # wird im Hauptprogramm gesetzt
v0[2][2] = 0
# Massen
m[0] = 1.988475415966536e+30
m[1] = 1.898568695e+27
m[2] = 0
def sq(x):
return x * x
def abstand(r, p, q):
return sqrt(sq(r[p][0] - r[q][0]) + sq(r[p][1] - r[q][1]) + sq(r[p][2] - r[q][2]))
def beschleunigung(N, m, r, a):
for p in range(N): # alle Beschleunigungen löschen
for k in range(3):
a[p][k] = 0
for p in range(N): # alle Planeten
for q in range(p + 1, N): # gegen alle anderen Planeten
R = abstand(r, p, q)
R3 = R * R * R
for k in range(3):
A = G * (r[q][k] - r[p][k]) / R3
a[p][k]+=m[q] * A
a[q][k]-=m[p] * A
def rkf5(tEnd, init):
global dt, rErrMin, rErrMax
t = 0
while(t < tEnd):
t = t + dt
beschleunigung(N, m, r0, a0)
for i in range(N):
for k in range(3):
r1[i][k] = r0[i][k] + 1 / 4 * dt * v0[i][k]
v1[i][k] = v0[i][k] + 1 / 4 * dt * a0[i][k]
beschleunigung(N, m, r1, a1)
for i in range(N):
for k in range(3):
r2[i][k] = r0[i][k] + 3 / 32 * dt * v0[i][k] + 9 / 32 * dt * v1[i][k]
v2[i][k] = v0[i][k] + 3 / 32 * dt * a0[i][k] + 9 / 32 * dt * a1[i][k]
beschleunigung(N, m, r2, a2)
for i in range(N):
for k in range(3):
r3[i][k] = r0[i][k] + 1932 / 2197 * dt * v0[i][k] + -7200 / 2197 * dt * v1[i][k] + \
7296 / 2197 * dt * v2[i][k]
v3[i][k] = v0[i][k] + 1932 / 2197 * dt * a0[i][k] + -7200 / 2197 * dt * a1[i][k] + \
7296 / 2197 * dt * a2[i][k]
beschleunigung(N, m, r3, a3)
for i in range(N):
for k in range(3):
r4[i][k] = r0[i][k] + 439 / 216 * dt * v0[i][k] - 8 * dt * v1[i][k] + 3680 / 513 * \
dt * v2[i][k] - 845 / 4104 * dt * v3[i][k]
v4[i][k] = v0[i][k] + 439 / 216 * dt * a0[i][k] - 8 * dt * a1[i][k] + 3680 / 513 * \
dt * a2[i][k] - 845 / 4104 * dt * a3[i][k]
beschleunigung(N, m, r4, a4)
for i in range(N):
for k in range(3):
r5[i][k] = r0[i][k] - 8 / 27 * dt * v0[i][k] + 2 * dt * v1[i][k] - 3544 / 2565 * \
dt * v2[i][k] + 1859 / 4104 * dt * v3[i][k] - 11 / 40 * dt * v4[i][k]
v5[i][k] = v0[i][k] - 8 / 27 * dt * a0[i][k] + 2 * dt * a1[i][k] - 3544 / 2565 * \
dt * a2[i][k] + 1859 / 4104 * dt * a3[i][k] - 11 / 40 * dt * a4[i][k]
beschleunigung(N, m, r5, a5)
for i in range(N):
for k in range(3):
r0_4[i][k] = r0[i][k]
# Ergebniss 4. Ordnung zur Fehlerberechnung
r0_4[i][k]+= dt * (25 / 216 * v0[i][k] + 0 * v1[i][k] + 1408 / 2565 * v2[i][k] + \
2197.0 / 4104 * v3[i][k] - 1 / 5 * v4[i][k] + 0 * v5[i][k])
# eigentliches Ergebnis 5. Ordnung
r0[i][k]+= dt * (16 / 135 * v0[i][k] + 0 * v1[i][k] + 6656 / 12825 * v2[i][k] + \
28561 / 56430 * v3[i][k] - 9 / 50 * v4[i][k] + 2 / 55 * v5[i][k])
v0[i][k]+= dt * (16 / 135 * a0[i][k] + 0 * a1[i][k] + 6656 / 12825 * a2[i][k] + \
28561 / 56430 * a3[i][k] - 9 / 50 * a4[i][k] + 2 / 55 * a5[i][k])
# Fehler ausrechnen. Benutzung erfordert eine while- statt einer for-Schleife
err = 0
for i in range(N): # welcher Körper gibt den größten Fehler?
localErr = 0
for k in range(3):
localErr+=sq(r0_4[i][k] - r0[i][k])
localErr = sqrt(localErr)
if (localErr > err):
err = localErr
if (init > 0): # Fehler wird in der ersten Runde so gesetzt,
rErrMin = err / 2 # wie es dem ersten dt entspricht
rErrMax = err
init = false
if (err < rErrMin): # Schrittweite verändern, falls der Fehler
dt*=1.2 # außerhalb der Grenzen
elif (err > rErrMax):
dt/=1.2
if (dt / sekProTag < 0.01):
Err = input("Schrittweitensteuerung versagt: dt erhöhen")
exit()
return t
# Hauptprogramm
startwerte()
rJ = r0[1][0] # der Jupiterradius ist konstant in dieser Rechnung
vJ = v0[1][1] # die Bahngeschwindigkeit auch
# Koordinatensystem einrichten
# x0 x1 xScale y0 y1 yScale
initKoor(-2.4,2.4,0.2, -2,2,0.2)
# alle im Text beschriebenen Beispiele
# X Y U V
# 0.2 0 0 2.89414 2:1
# 0.9 0 0 0.38 3:1
# 0.9 0 0 0.796 2:1
# 0.495 0.874685658 -0.8653265 0.4930764 1:1
# 0.495 0.874685658 -0.8653265 -0.4930764 1:1
# 0.45 0.953 -0.846264 0.428626 1:1
# 0.45 0.953 -0.846264 -0.428626 1:1
# 0.6 1.04 -0.72176 0.38641 1:1
# 0.6 1.04 -0.72176 0.38641 1:1
# die Ausgangswert in Jupitierradien und Jupitergeschwindigkeiten
X=0.45
Y=0.953
U=-0.846264
V=0.428626
T = 1 # simulierte Zeit in Jupiterumläufen (1)
P = 1000 # Anzahl Zwischenwerte (1000)
dt = 2.5 # 2,5 Tage
# Startwerte in die Felder einsetzen
r0[2][0] = X * rJ
v0[2][1] = V * v0[1][1]
r0[2][1] = Y * rJ
v0[2][0] = U * v0[1][1]
# den anderen Teil der Startwerte in SI umwandeln
dt = dt * sekProTag
UJup = 2 * pi * rJ / v0[1][1]
T = T * UJup
# Punkte für Jupiter und Sonne
Plot(r0[0][0] / rJ,r0[0][1] / rJ,20, "gold")
Plot(r0[1][0] / rJ,r0[1][1] / rJ,8,"coral")
Plot(0.5,sqrt(3)/2,5,"blue")
init = true
t = 0
while (t < T):
t = t + rkf5(T / P, init)
init = false # nur in der ersten Runde soll die interne Fehlergrenze gesetzt werden.
# in mitrotierene Koordinaten umrechnen (Skalarprodukt)
x = r0[2][0]
y = r0[2][1]
xJ = r0[1][0]
yJ = r0[1][1]
rP = sqrt(x * x + y * y)
xP = xJ * x + yJ * y
xP = xP / rJ / rJ
yP = -yJ * x + xJ * y
yP = yP / rJ / rJ
Plot(xP,yP,3,"black") # mitrotierende Koordinaten
#Plot(x/rJ,y/rJ,2,"black") # normale Koordinaten Asteroid
#Plot (xJ/rJ,yJ/rJ,2,"coral") # ... und Jupiter
update() # einkommentieren um den Bahnverlauf zu sehen
update()
done()