Startseite | Programme: Mehrkörperproblem | Algorithmen, Lektionen, Analysen | Unsere Vorläufer | Kontakt, Datenschutz

VdS-Journal 73, 74: Das Mehrkörperproblem, Beispiel Sonnensystem

Das Mehrkörperproblem ist eines der zentralen Untersuchungsgegenstände der Astronomie: Von der Positionsberechnung über Stabilitätsuntersuchungen bis hin zur dynamischen Entwicklung von Galaxienhaufen. In dieser ersten Aufsatz-Serie erläutere ich die theoretischen Grundlagen und gebe ein wirklich gut funktionierendes Programm.
Python ist eine Interpretersprache und recht langsam. Wer große Probleme rechnen möchte und sich ein wenig mit der Programmierung auskennt, möge sich an mich wenden. Ich habe eine C++-Version.

Schnell-Anwahl:
VdS-Journal 73 – Das Euler-Cauchy-Verfahren
VdS-Journal 74 – Das Heun-Verfahren
VdS-Journal 74 – Das Runge-Kutta-Fehlberg-Verfahren

VdS-Journal 73: Das Euler-Cauchy-Verfahren

Dieses Verfahren ist sehr leicht verständlich, es benötigt  aber sehr viel Rechenzeit für genaue Resultate.

# Programm EulerCauchy.py
from turtle import *
from random import *
from math import *

# globale Variable
G = 6.67408e-11         #Gravitationskonstante
mProAE=149597870700    #Astronomische Einheit
sekProTag=86400
# sekProJahr=365.2425*86400
# mProLj=9460730472580800

N=10 # Anzahl der Körper
r0 = [[0 for i in range(3)] for j in range(10)]
v0 = [[0 for i in range(3)] for j in range(10)]
a0 = [[0 for i in range(3)] for j in range(10)]
m = [0 for i in range (10)]

def startwerte():
    # Orte und Geschwindigkeiten
    r0[0][0]=   -1063436181;    v0[0][0]=      9.29323348
    r0[0][1]=    -403178315;    v0[0][1]=    -12.81658147
    r0[0][2]=      27924311;    v0[0][2]=     -0.15918117
    r0[1][0]=  -20525703382;    v0[1][0]=  37004.00242937
    r0[1][1]=  -67316663267;    v0[1][1]= -11177.58979901
    r0[1][2]=   -3651721547;    v0[1][2]=  -4307.67352332
    r0[2][0]= -108522074771;    v0[2][0]=   1392.28564030
    r0[2][1]=   -5295518213;    v0[2][1]= -35152.38582183
    r0[2][2]=    6163650722;    v0[2][2]=   -560.19899547
    r0[3][0]=  -27570384049;    v0[3][0]= -29777.08241127
    r0[3][1]=  144289541085;    v0[3][1]=  -5492.11604817
    r0[3][2]=      27924311;    v0[3][2]=     -0.15918117
    r0[4][0]=  206972399720;    v0[4][0]=   1173.70663100
    r0[4][1]=   -2404931177;    v0[4][1]=  26284.74809179
    r0[4][2]=   -5126803613;    v0[4][2]=    522.10695479
    r0[5][0]=  597504078214;    v0[5][0]=  -7900.54555615
    r0[5][1]=  439201609603;    v0[5][1]=  11143.31683419
    r0[5][2]=  -15198533263;    v0[5][2]=    130.70502469
    r0[6][0]=  957322522683;    v0[6][0]=  -7422.72126397
    r0[6][1]=  982453470791;    v0[6][1]=   6723.09678597
    r0[6][2]=  -55183935225;    v0[6][2]=    178.08847902
    r0[7][0]= 2157914958773;    v0[7][0]=   4646.31018375
    r0[7][1]=-2055025724501;    v0[7][1]=   4614.84545075
    r0[7][2]=  -35599134212;    v0[7][2]=    -43.05247532
    r0[8][0]= 2513984445571;    v0[8][0]=   4475.19788359
    r0[8][1]=-3739120513758;    v0[8][1]=   3063.80412436
    r0[8][2]=   19061170936;    v0[8][2]=   -166.22362888
    r0[9][0]=-1478397309067;    v0[9][0]=   5253.44576775
    r0[9][1]=-4182980594505;    v0[9][1]=  -2675.42458857
    r0[9][2]=  875240628161;    v0[9][2]=  -1233.31598335
    # Massen
    m[0]=1.9884754159665367161e+30; m[1]=3.3011412045397049632e+23
    m[2]=4.8674688627864414779e+24; m[3]=6.0458266739227718677e+24
    m[4]=6.4171071702951773445e+23; m[5]=1.8985686953960564018e+27
    m[6]=5.6837942431514554871e+26; m[7]=8.6950693776139604302e+25
    m[8]=1.0295513181974406248e+26; m[9]=6.6282513865551224955e+24

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 eulerCauchy(tsteps, dt):
    for t in range(tsteps):
        beschleunigung(N, m, r0, a0);						
        for i in range(N):
            for k in range(3):
                r0[i][k] = r0[i][k] + dt * v0[i][k]	
                v0[i][k] = v0[i][k] + dt * a0[i][k]

def ausgabe(ii):
    for i in range(4,5):
        x=r0[i][0]-r0[0][0]
        y=r0[i][1]-r0[0][1]
        z=r0[i][2]-r0[0][2]
        r=sqrt(x*x+y*y+z*z)
        l=atan2(y,x)
        if (l<0): l+=2*pi;
        b=atan(z/r)
        print(ii, r/mProAE, 180*l/pi, 180*b/pi)

    
# Hauptprogramm
startwerte()
for i in range(400):
    eulerCauchy(10, 0.1*sekProTag)
ausgabe(i+1)

name = input("Fertig?")

VdS-Journal 74: Das Heun-Verfahren

Das Heun-Verfahren ist immer noch leicht verständlich, konvergiert aber wesentlich schneller.

# Programm Heun.py
from turtle import *
from random import *
from math import *

# globale Variable
G = 6.67408e-11         #Gravitationskonstante
mProAE=149597870700    #Astronomische Einheit
sekProTag=86400
# sekProJahr=365.2425*86400
# mProLj=9460730472580800

N=10 # Anzahl der Körper
r0 = [[0 for i in range(3)] for j in range(10)]
v0 = [[0 for i in range(3)] for j in range(10)]
a0 = [[0 for i in range(3)] for j in range(10)]
r1 = [[0 for i in range(3)] for j in range(10)]
v1 = [[0 for i in range(3)] for j in range(10)]
a1 = [[0 for i in range(3)] for j in range(10)]

m = [0 for i in range (10)]

def startwerte():
    # Orte und Geschwindigkeiten
    r0[0][0]=   -1063436181;    v0[0][0]=      9.29323348
    r0[0][1]=    -403178315;    v0[0][1]=    -12.81658147
    r0[0][2]=      27924311;    v0[0][2]=     -0.15918117
    r0[1][0]=  -20525703382;    v0[1][0]=  37004.00242937
    r0[1][1]=  -67316663267;    v0[1][1]= -11177.58979901
    r0[1][2]=   -3651721547;    v0[1][2]=  -4307.67352332
    r0[2][0]= -108522074771;    v0[2][0]=   1392.28564030
    r0[2][1]=   -5295518213;    v0[2][1]= -35152.38582183
    r0[2][2]=    6163650722;    v0[2][2]=   -560.19899547
    r0[3][0]=  -27570384049;    v0[3][0]= -29777.08241127
    r0[3][1]=  144289541085;    v0[3][1]=  -5492.11604817
    r0[3][2]=      27924311;    v0[3][2]=     -0.15918117
    r0[4][0]=  206972399720;    v0[4][0]=   1173.70663100
    r0[4][1]=   -2404931177;    v0[4][1]=  26284.74809179
    r0[4][2]=   -5126803613;    v0[4][2]=    522.10695479
    r0[5][0]=  597504078214;    v0[5][0]=  -7900.54555615
    r0[5][1]=  439201609603;    v0[5][1]=  11143.31683419
    r0[5][2]=  -15198533263;    v0[5][2]=    130.70502469
    r0[6][0]=  957322522683;    v0[6][0]=  -7422.72126397
    r0[6][1]=  982453470791;    v0[6][1]=   6723.09678597
    r0[6][2]=  -55183935225;    v0[6][2]=    178.08847902
    r0[7][0]= 2157914958773;    v0[7][0]=   4646.31018375
    r0[7][1]=-2055025724501;    v0[7][1]=   4614.84545075
    r0[7][2]=  -35599134212;    v0[7][2]=    -43.05247532
    r0[8][0]= 2513984445571;    v0[8][0]=   4475.19788359
    r0[8][1]=-3739120513758;    v0[8][1]=   3063.80412436
    r0[8][2]=   19061170936;    v0[8][2]=   -166.22362888
    r0[9][0]=-1478397309067;    v0[9][0]=   5253.44576775
    r0[9][1]=-4182980594505;    v0[9][1]=  -2675.42458857
    r0[9][2]=  875240628161;    v0[9][2]=  -1233.31598335
    # Massen
    m[0]=1.9884754159665367161e+30; m[1]=3.3011412045397049632e+23
    m[2]=4.8674688627864414779e+24; m[3]=6.0458266739227718677e+24
    m[4]=6.4171071702951773445e+23; m[5]=1.8985686953960564018e+27
    m[6]=5.6837942431514554871e+26; m[7]=8.6950693776139604302e+25
    m[8]=1.0295513181974406248e+26; m[9]=6.6282513865551224955e+24

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 heun(tsteps, dt):
    for t in range(tsteps):
        beschleunigung(N, m, r0, a0);						
        for i in range(N):
            for k in range(3):
                r1[i][k] = r0[i][k] + dt * v0[i][k]	
                v1[i][k] = v0[i][k] + dt * a0[i][k]
        beschleunigung(N, m, r1, a1);						
        for i in range(N):
            for k in range(3):
                r0[i][k] = r0[i][k] + dt * ( 0.5*v0[i][k] + 0.5*v1[i][k] )
                v0[i][k] = v0[i][k] + dt * ( 0.5*a0[i][k] + 0.5*a1[i][k] )
                
                

def ausgabe(ii):
    for i in range(4,5):
        x=r0[i][0]-r0[0][0]
        y=r0[i][1]-r0[0][1]
        z=r0[i][2]-r0[0][2]
        r=sqrt(x*x+y*y+z*z)
        l=atan2(y,x)
        if (l<0): l+=2*pi;
        b=atan(z/r)
        print(ii, r/mProAE, 180*l/pi, 180*b/pi)

    
# Hauptprogramm
startwerte()
for i in range(400):
    heun(1, 1*sekProTag)
ausgabe(i+1)

name = input("Fertig?")

VdS-Journal 74: Das Runge-Kutta-Fehlberg-Verfahren

Runge-Kutta-Verfahren bilden eine ganze Klasse sehr leistungsfähiger Methoden. Die Modifikation nach Erwin Fehlberg erlaubt eine automatische Schrittweitensteuerung. Diese wird hier nicht benutzt, ich werde sie aber für künftige Fragestellungen benötigen.

# Programm RKF5.py
from turtle import *
from random import *
from math import *

# globale Variable
G = 6.67408e-11         #Gravitationskonstante
mProAE=149597870700    #Astronomische Einheit
sekProTag=86400
# sekProJahr=365.2425*86400
# mProLj=9460730472580800

N=10 # Anzahl der Körper
r0 = [[0 for i in range(3)] for j in range(10)]
v0 = [[0 for i in range(3)] for j in range(10)]
a0 = [[0 for i in range(3)] for j in range(10)]
r1 = [[0 for i in range(3)] for j in range(10)]
v1 = [[0 for i in range(3)] for j in range(10)]
a1 = [[0 for i in range(3)] for j in range(10)]
r2 = [[0 for i in range(3)] for j in range(10)]
v2 = [[0 for i in range(3)] for j in range(10)]
a2 = [[0 for i in range(3)] for j in range(10)]
r3 = [[0 for i in range(3)] for j in range(10)]
v3 = [[0 for i in range(3)] for j in range(10)]
a3 = [[0 for i in range(3)] for j in range(10)]
r4 = [[0 for i in range(3)] for j in range(10)]
v4 = [[0 for i in range(3)] for j in range(10)]
a4 = [[0 for i in range(3)] for j in range(10)]
r5 = [[0 for i in range(3)] for j in range(10)]
v5 = [[0 for i in range(3)] for j in range(10)]
a5 = [[0 for i in range(3)] for j in range(10)]
r0_4 = [[0 for i in range(3)] for j in range(10)] # Ergebnis für RKF 4. Ordnung
m = [0 for i in range (10)]

def startwerte():
    # Orte und Geschwindigkeiten
    r0[0][0]=   -1063436181;    v0[0][0]=      9.29323348
    r0[0][1]=    -403178315;    v0[0][1]=    -12.81658147
    r0[0][2]=      27924311;    v0[0][2]=     -0.15918117
    r0[1][0]=  -20525703382;    v0[1][0]=  37004.00242937
    r0[1][1]=  -67316663267;    v0[1][1]= -11177.58979901
    r0[1][2]=   -3651721547;    v0[1][2]=  -4307.67352332
    r0[2][0]= -108522074771;    v0[2][0]=   1392.28564030
    r0[2][1]=   -5295518213;    v0[2][1]= -35152.38582183
    r0[2][2]=    6163650722;    v0[2][2]=   -560.19899547
    r0[3][0]=  -27570384049;    v0[3][0]= -29777.08241127
    r0[3][1]=  144289541085;    v0[3][1]=  -5492.11604817
    r0[3][2]=      27924311;    v0[3][2]=     -0.15918117
    r0[4][0]=  206972399720;    v0[4][0]=   1173.70663100
    r0[4][1]=   -2404931177;    v0[4][1]=  26284.74809179
    r0[4][2]=   -5126803613;    v0[4][2]=    522.10695479
    r0[5][0]=  597504078214;    v0[5][0]=  -7900.54555615
    r0[5][1]=  439201609603;    v0[5][1]=  11143.31683419
    r0[5][2]=  -15198533263;    v0[5][2]=    130.70502469
    r0[6][0]=  957322522683;    v0[6][0]=  -7422.72126397
    r0[6][1]=  982453470791;    v0[6][1]=   6723.09678597
    r0[6][2]=  -55183935225;    v0[6][2]=    178.08847902
    r0[7][0]= 2157914958773;    v0[7][0]=   4646.31018375
    r0[7][1]=-2055025724501;    v0[7][1]=   4614.84545075
    r0[7][2]=  -35599134212;    v0[7][2]=    -43.05247532
    r0[8][0]= 2513984445571;    v0[8][0]=   4475.19788359
    r0[8][1]=-3739120513758;    v0[8][1]=   3063.80412436
    r0[8][2]=   19061170936;    v0[8][2]=   -166.22362888
    r0[9][0]=-1478397309067;    v0[9][0]=   5253.44576775
    r0[9][1]=-4182980594505;    v0[9][1]=  -2675.42458857
    r0[9][2]=  875240628161;    v0[9][2]=  -1233.31598335
    # Massen
    m[0]=1.9884754159665367161e+30; m[1]=3.3011412045397049632e+23
    m[2]=4.8674688627864414779e+24; m[3]=6.0458266739227718677e+24
    m[4]=6.4171071702951773445e+23; m[5]=1.8985686953960564018e+27
    m[6]=5.6837942431514554871e+26; m[7]=8.6950693776139604302e+25
    m[8]=1.0295513181974406248e+26; m[9]=6.6282513865551224955e+24

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(tsteps, dt):
    for t in range(tsteps):
        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]
                # Ergebnisse 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] \
                )
                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. Damit kann man nur etwas anfangen,
        # wenn man die for-Schleife durch eine while-Schleife ersetzt
        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;



rMin = [0 for i in range(N)]
rMax = [0 for i in range(N)]

# Hauptprogramm
startwerte()
for i in range(2000):
    # minimale und maximalen Sonnenabstand löschen
    for j in range(N):
        rMin[j]=7e77
        rMax[j]=0
    # je 100 Jahre rechnen
    for j in range(36524):     # 36524 Schleifen ...
        rkf5(2, 0.5*sekProTag) # .. zu 1 Tagen = 36524 Tage 
    # minimalen und maximalen Sonnenabstand ermitteln
        for j in range(N):
            r=sqrt( sq(r0[j][0]) + sq(r0[j][1]) + sq(r0[j][2]))
            if (rMin[j]>r):
                rMin[j]=r;
            if (rMax[j]<r):
                rMax[j]=r;
    # minimalen und maximalen Sonnenabstand ausgeben
    print("Minimum, t=",i/10," : ",end='')
    for j in range(N):
        print (rMin[j]/mProAE," ",end='')
    print("")
    print("Maximum, t=",i/10," : ",end='')
    for j in range(N):
        print (rMax[j]/mProAE," ",end='')
    print("")

name = input("Fertig?")