Hufeisenbahnen und Integration in C.

8

Ich untersuche einen bestimmten Fall des eingeschränkten Drei-Körper-Problems. Es wurde festgestellt, dass einige Objekte einem Hufeisenbahnmuster folgen, und ich versuche, etwas durch einen Integrationscode in C zu sortieren. Ich folge einigen Ratschlägen im Artikel Familien periodischer Hufeisenbahnen im eingeschränkten Drei-Körper-Problem . Das gibt mir ideale Anfangsbedingungen und die Gleichungen im Schwerpunktsystem. (m ist die Masse der Erde und die sich daraus ergebende Position der Sonne im Massenschwerpunkt-Bezugssystem (x, y) sind die Koordinaten des dritten Körpers, die als masselos angenommen werden (wie es das eingeschränkte Problem erfordert).

Ö=(x2+y2)/.2+(1- -m)r1+mr2+(1- -m)m2
r12=(x- -m)2+y2

r22=(x- -m+1)2+y2
ein(x)=dÖdx+2v(y)
ein(y)=dÖdy- -2v(x)

Die Positionen von "Sonne" und "Erde" sind im gleichen Bezugssystem auf (m, 0) und (m-1,0) festgelegt. (rotierendes Bezugssystem, vorausgesetzt, die Erde hat eine kreisförmige Umlaufbahn.)

Aus all dem habe ich die Gleichungen berechnet, um das System zu beschreiben:

ein(x)=x+(m- -1)(x- -m)((x- -m)2+y2)1.5- -2m(x- -m+1)((x- -m+1)2+y2)1.5+2v(y)
ein(y)=y- -y(1- -m)((x- -m)2+y2)1.5- -2ym((x- -m+1)2+y2)1.5- -2v(x)

Ich habe den Algorithmus von Runge-Kutta 4 verwendet, um diese Gleichungen zu integrieren. (Ich weiß, dass der Code ziemlich umwerfend ist, aber ich kann einfach keine Zeiger verwenden und ich verwende überall Strukturen).

#include<stdio.h>
#include<stdlib.h>
#include<math.h>

#define dt 0.0001
#define N 100

typedef struct{
    long double x,y;
}vec;

typedef struct{
    vec k1,k2,k3,k4;
}runge;

typedef struct{
    runge r,v;
}big;

double dS,dE,m;

double accx(double,double,double);
double accy(double,double,double);
void rad(vec);
big rungekutta(vec,vec);
vec moto(vec,runge);
double jacobi(vec);

int main(){
  vec r,v;
    big f;
    double J,t;
    int i,Num;
    FILE* s1;
    s1=fopen("HorseShoe.dat","w");

    Num=(int)N/dt;
    scanf("%Lf",&r.x);
    scanf("%Lf",&r.y);
    scanf("%Lf",&v.x);
    scanf("%Lf",&v.y);
    scanf("%lf",&m);

    for(i=0;i<Num;i++){
        t=(i+1)*dt;
        rad(r);
        f=rungekutta(r,v);
        r=moto(r,f.r);
        v=moto(v,f.v);
        J=jacobi(r);
        fprintf(s1,"%lf\t%Lf\t%Lf\t%Lf\t%Lf\t%lf\n",t,r.x,r.y,v.x,v.y,J);
    }
return 0;
}

void rad(vec r){
    dS=pow(r.x-m,2)+pow(r.y,2);
    dE=pow(r.x-m+1,2)+pow(r.y,2);
}

double jacobi(vec r){
    return pow(r.x,2)+pow(r.y,2)+2*(1-m)/dS+2*m/dE+m*(1-m);
}

double accx(double x,double y,double v){
    return x-(x-m)*(1-m)/pow(pow(x-m,2)+pow(y,2),1.5)-m*(x-m+1)/pow(pow(x-m+1,2)+pow(y,2),1.5)+2*v;
}

double accy(double x,double y,double v){
    return y-(1-m)*y/pow(pow(y,2)+pow(x-m,2),1.5)-m*y/pow(pow(y,2)+pow(x-m+1,2),1.5)-2*v;
}

big rungekutta(vec r,vec v){
    big f;
    f.r.k1.x=v.x;
    f.r.k1.y=v.y;
    f.v.k1.x=accx(r.x,r.y,v.y);
    f.v.k1.y=accy(r.x,r.y,v.x);
    f.r.k2.x=v.x+f.v.k1.x*dt/2;
    f.r.k2.y=v.y+f.v.k1.y*dt/2;
    f.v.k2.x=accx(r.x+f.r.k1.x*dt/2,r.y+f.r.k1.y*dt/2,v.y+f.v.k1.y*dt/2);
    f.v.k2.y=accy(r.x+f.r.k1.x*dt/2,r.y+f.r.k1.y*dt/2,v.x+f.v.k1.x*dt/2);
    f.r.k3.x=v.x+f.v.k2.x*dt/2;
    f.r.k3.y=v.y+f.v.k2.y*dt/2;
    f.v.k3.x=accx(r.x+f.r.k2.x*dt/2,r.y+f.r.k2.y*dt/2,v.y+f.v.k2.y*dt/2);
    f.v.k3.y=accy(r.x+f.r.k2.x*dt/2,r.y+f.r.k2.y*dt/2,v.x+f.v.k2.x*dt/2);
    f.r.k4.x=v.x+f.v.k3.x*dt;
    f.r.k4.y=v.y+f.v.k3.y*dt;
    f.v.k4.x=accx(r.x+f.r.k3.x*dt,r.y+f.r.k3.y*dt,v.y+f.v.k3.y*dt);
    f.v.k4.y=accy(r.x+f.r.k3.x*dt,r.y+f.r.k3.y*dt,v.x+f.v.k3.x*dt);
    return f;
}

vec moto(vec r,runge rk){
    r.x+=(rk.k1.x+2*rk.k2.x+2*rk.k3.x+rk.k4.x)*dt/6;
    r.y+=(rk.k1.y+2*rk.k2.y+2*rk.k3.y+rk.k4.y)*dt/6;
    return r;
}

Beim Zeichnen der Ergebnisse erhalte ich nur eine Spirale, während ich bei Verwendung der angegebenen Eingaben eine Hufeisenbahn erhalten sollte. Ich habe viele verschiedene Eingaben versucht (m = 0,0001 und m = 0,000003, letzteres entspricht den tatsächlichen Werten der Erd- und Sonnenmasse (Sonnenmasse ist 1 m)).

Ich kann einfach nicht sehen, was los ist (wahrscheinlich alles: D). Kann mir bitte jemand helfen?

Elisa
quelle
Menschen können möglicherweise besser helfen, wenn Sie erklären, was Ihre Variablen (in den Gleichungen) sind: O, x, y, m, r1 / 2 usw. Außerdem macht Ihre Frage nicht klar, was das Problem ist, dass Sie sind Begegnung. "Ich kann nicht sehen, was los ist" - was sollst du bekommen und was bekommst du stattdessen?
Alex
Bevor ich fortfahre, möchte ich zwei Dinge sehen: Ein Link zu dem Referenzpapier, das Sie verwendet haben, um zu sehen, ob Sie einen Fehler in Ihrer Mathematik gemacht haben, und Ihre Eingaben. Ich vermute, Ihre Eingaben sind möglicherweise ausgeschaltet. Dies ist ein nicht standardmäßiges Einheitensystem, in dem die Masse der Erde angegeben istm sollte etwa sein 3×10- -6ist die Masse der Sonne 1 weniger als diese kleine Zahl und die Größe des Positionsvektors (x,y)ist ungefähr 1. Die Zeit scheint ebenfalls skaliert zu sein. Dies ist kein Ort für SI-Einheiten.
David Hammen
Ich habe gerade die Frage bearbeitet und die von Ihnen gestellten Details hinzugefügt (Entschuldigung, wenn ich etwas chaotisch bin, erste Frage hier). Ich habe G = 1 verwendet, natürlich konnte ich keine SI-Einheiten verwenden;)
Elisa
5
Sie können gerne eine Antwort auf Ihre eigene Frage veröffentlichen, in der genau angegeben ist, was falsch war und wie Sie es behoben haben. Dies könnte sich in Zukunft für andere als nützlich erweisen, die zu dieser Frage kommen!
Zephyr
2
Nur eine kleine, aber wichtige Bemerkung: Runge-Kutta 4ᵗʰ Ordnung ist nicht symplektisch , was bedeutet, dass es die gesamte Orbitalenergie im Laufe der Zeit verändert, zunächst langsam, aber exponentiell an Schwere zunimmt. Abhängig von Ihrem gewünschten Integrationsintervall sollten Sie RK4 zugunsten von etwas bevorzugen, das besser für die Himmelsmechanik geeignet ist. In diesem Dokument finden Sie beispielsweise Hintergrundinformationen und Empfehlungen.
Rody Oldenhuis

Antworten:

1

Ich werde weder Ihre Gleichungen noch Ihren Code sorgfältig durchgehen, aber da Sie die genauen Startpositionen, die Sie verwendet haben, nicht erwähnt haben, ist meine Vermutung, dass Sie nicht zuerst nach einem stabilen Hufeisen-Orbit- Zustandsvektor gesucht haben.

Die Positionen für Sonne und Erde sind im synodischen (rotierenden) Rahmen festgelegt. Aber wo hast du dein Astroid angefangen? Für eine bestimmte Startposition müssen Sie dem Asteroiden auch die richtige Startgeschwindigkeit (Geschwindigkeit und Richtung) geben, da er sonst möglicherweise nicht stabil ist und einfach abwandert oder sich wie erwähnt dreht.

Schauen Sie sich diese Antwort an, in der ich die richtigen Gleichungen zeige und dann nach stabilen Hufeisenbahnen suche, bevor ich sie zeichne.

Ich beginne mit einer Position auf der x-Achse gegenüber der Erde und etwas weiter oder näher von der Sonne entfernt als die Erde. Dann starte ich es mit einer geringen Geschwindigkeit (im rotierenden Rahmen) genau in Richtung der y-Achse. Ich beobachte, wie es zur Erde driftet und dann zurück in die Region des Startgebiets bumerang. Wenn es die x-Achse wieder kreuzt, überprüfe ich, wie nahe die Geschwindigkeit an y liegt. Ich benutze einen Optimierungs-Nulllöser brentq, um die Anfangsgeschwindigkeit anzupassen, bis die Rücklaufgeschwindigkeit genau in die entgegengesetzte Richtung ist.

Wenn dies zutrifft, ist die Umlaufbahn periodisch und wiederholt sich, und wir können sie unter den Bedingungen des CR3BP als veraltete Umlaufbahn bezeichnen.


Aus dieser Antwort (es gibt dort viel mehr zu lesen, einschließlich mehrerer Referenzen!):

x¨=x+2y˙- -(1- -μ)(x+μ)r13- -μ(x- -1+μ)r23
y¨=y- -2x˙- -(1- -μ)yr13- -μyr23
z¨=- -(1- -μ)zr13- -μzr23

Geben Sie hier die Bildbeschreibung ein

oben: Halbzyklen einiger wackeliger Hufeisenbahnen

Geben Sie hier die Bildbeschreibung ein

oben: Zeiten bis zu den ersten x-Achsen-Kreuzungen derselben wackeligen Hufeisenbahnen, die zur Berechnung der Halbzykluszeiten verwendet werden.

Geben Sie hier die Bildbeschreibung ein

oben: Zykluszeiten aus dieser Berechnung (schwarze Punkte) im Vergleich zur Methode zur Schätzung der Synodenperiode (rote Punkte). Gute qualitative Übereinstimmung. Auch die Start-y-Geschwindigkeiten an jedem Startpunkt in x.

unten: Python-Skript für diese Diagramme.

def x_acc(x, ydot):
    r1    = np.abs(x-x1)
    r2    = np.abs(x-x2)
    xddot = x + 2*ydot  -  ((1-mu)/r1**3)*(x+mu) - (mu/r2**3)*(x-(1-mu))
    return xddot

def C_calc(x, y, z, xdot, ydot, zdot):
    r1 = np.sqrt((x-x1)**2 + y**2 + z**2)
    r2 = np.sqrt((x-x2)**2 + y**2 + z**2)
    C = (x**2 + y**2 + 2.*(1-mu)/r1 + 2.*mu/r2 - (xdot**2 + ydot**2 + zdot**2))
    return C

def deriv(X, t): 
    x, y, z, xdot, ydot, zdot = X
    r1 = np.sqrt((x-x1)**2 + y**2 + z**2)
    r2 = np.sqrt((x-x2)**2 + y**2 + z**2)
    xddot = x + 2*ydot  -  ((1-mu)/r1**3)*(x+mu) - (mu/r2**3)*(x-(1-mu))
    yddot = y - 2*xdot  -  ((1-mu)/r1**3)*y      - (mu/r2**3)*y
    zddot =             -  ((1-mu)/r1**3)*z      - (mu/r2**3)*z
    return np.hstack((xdot, ydot, zdot, xddot, yddot, zddot))

# http://cosweb1.fau.edu/~jmirelesjames/hw4Notes.pdf

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

halfpi, pi, twopi = [f*np.pi for f in (0.5, 1, 2)]

mu = 0.001

x1 = -mu
x2 = 1. - mu

x = np.linspace(-1.4, 1.4, 1201)
y = np.linspace(-1.4, 1.4, 1201)

Y, X = np.meshgrid(y, x, indexing='ij')
Z    = np.zeros_like(X)

xdot, ydot, zdot = [np.zeros_like(X) for i in range(3)]

C = C_calc(X, Y, Z, xdot, ydot, zdot)
C[C>8] = np.nan

if True:
    plt.figure()
    plt.imshow(C)
    plt.colorbar()
    levels = np.arange(2.9, 3.2, 0.04) 
    CS = plt.contour(C, levels,
                 origin='lower',
                 linewidths=2) 
    plt.show()

ydot0s   = np.linspace(-0.08, 0.08, 20)
x0ydot0s = []
for ydot0 in ydot0s:
    x0, infob =  brentq(x_acc, -1.5, -0.5, args=(ydot0), xtol=1E-11, rtol=1E-11,
                           maxiter=100, full_output=True, disp=True)
    x0ydot0s.append((x0, ydot0))

states = [np.array([x0, 0, 0, 0, ydot0, 0]) for (x0, ydot0) in x0ydot0s]

times  = np.arange(0, 150, 0.01)

results = []
for X0 in states:
    answer, info = ODEint(deriv, X0, times, atol = 1E-11, full_output=True)
    results.append(answer.T.copy())

resultz = []
for x0ydot0, thing in zip(x0ydot0s, results):
    y     = thing[1]
    check = y[2:]*y[1:-1] < 0
    zc    = np.argmax(y[2:]*y[1:-1] < 0) + 1
    if zc > 10:
        resultz.append((thing, zc, x0ydot0))

if True:
    plt.figure()
    hw = 1.6
    for j, (thing, zc, x0ydot0) in enumerate(resultz):
        x, y = thing[:2,:zc]
        plt.plot(x, y)
    plt.xlim(-hw, hw)
    plt.ylim(-hw, hw)
    plt.plot([x1], [0], 'ok')
    plt.plot([x2], [0], 'ok')
    plt.show()

if True:
    plt.figure()
    for j, (thing, zc, x0ydot0) in enumerate(resultz):
        x, y = thing[:2]
        plt.plot(times[:zc], y[:zc])
    plt.show()

if True:
    plt.figure()
    for j, (thing, zc, x0ydot0) in enumerate(resultz):
        x0, ydot0 = x0ydot0
        cycle_time = 2. * times[zc] / twopi
        ratio = abs(x0/x2)
        T_simple_model = twopi * abs(x0/x2)**1.5
        T_synodic_simple_model = 1. / (1. - twopi/T_simple_model) # https://astronomy.stackexchange.com/a/25002/7982
        plt.subplot(2, 1, 1)
        plt.plot(x0, cycle_time, 'ok')
        plt.plot(x0, abs(T_synodic_simple_model), 'or')
        plt.subplot(2, 1, 2)
        plt.plot(x0, ydot0, 'ok')
    plt.subplot(2, 1, 1)
    plt.xlabel('x0', fontsize=16)
    plt.ylabel('cycle times (periods)', fontsize=16)
    plt.subplot(2, 1, 2)
    plt.xlabel('x0', fontsize=16)
    plt.ylabel('ydot0', fontsize=16)
    plt.show()
uhoh
quelle