La simulación de N-cuerpos sigue perdiendo precisión a pesar de utilizar un integrador simpléctico y aritmético de precisión arbitraria

Hay varias causas posibles para que un simulador de órbita se comporte mal, incluso cuando se utiliza un excelente paquete de precisión arbitraria como mpmath . Al proporcionar valores numéricos, debe pasarlos a mpmath como cadenas o números enteros, no flotantes, para que pueda convertirlos correctamente. Consulte los documentos de mpmath para obtener más detalles.

Otra posible causa de problemas es la cancelación catastrófica que puede ocurrir al combinar números de tamaños muy diferentes. A menudo, puede minimizar ese problema organizando cuidadosamente el orden de los cálculos y trabajando en una escala en la que el rango de magnitudes de sus números no sea demasiado grande. Por supuesto, cuando usa una biblioteca como mpmath, simplemente puede aumentar la precisión, aunque eso ralentizará las cosas.

Lo estándar que se debe hacer cuando un simulador de órbita codificado correctamente está dando errores inaceptablemente altos es hacer que el paso de tiempo sea más pequeño. Sin embargo, el uso de una gran cantidad de pasos de tiempo aumenta los errores de redondeo acumulativos, por lo que si hace que el paso de tiempo sea demasiado pequeño, el resultado será peor, como mencionó anteriormente David Hammen. Por supuesto, al usar mpmath puede solucionar ese problema aumentando la precisión.

Como referencia, aquí hay un código de Python que realiza una simulación de órbita de 1 cuerpo simple utilizando la forma sincronizada de integración de Leapfrog , que en realidad es lo mismo que el algoritmo Verlet de velocidad . Utiliza la aritmética de punto flotante de Python simple, no mpmath. El floattipo Python usa el formato IEEE-754 binary64 , que tiene 53 bits de precisión (casi 16 dígitos decimales).

Este programa utiliza unidades naturales, es decir, el parámetro gravitatorio del cuerpo central es$GM=\mu=4\pi^2$, por lo que una órbita circular de radio 1 unidad espacial tiene un período de 1 unidad de tiempo.

Para una órbita circular simple, utilizando 10 000 pasos de tiempo por órbita, el error máximo en el radio es de menos de 2 partes en 10 000 000. Se requiere un paso de tiempo más pequeño para órbitas excéntricas.

''' Simple orbit sim using the
    synchronised Leapfrog algorithm
    Written by PM 2Ring 2017.05.20
'''

from math import pi

# Gravitational parameter of central body in "natural" units,
# i.e. where R**3 == T**2
mu = 4.0 * pi * pi

# Speed of a body in a circular orbit of radius r
def speed(r):
    return (mu / r) ** 0.5

# Acceleration vector due to central gravity at (x, y)
def acc(x, y):
    a = -mu / (x * x + y * y) ** 1.5
    return a * x, a * y

class Body(object):
    ''' An orbiting body '''
    def __init__(self, x, y, vx, vy, delta_time):
        # Current position
        self.x, self.y = x, y

        # Current velocity
        self.vx, self.vy = vx, vy

        # Time step
        self.delta_time = delta_time

        # Current acceleration due to central gravity
        self.ax, self.ay = acc(x, y)

        self.points = [(x, y)]

    def update(self):
        ''' Update the body's orbit parameters using
            the synchronised leapfrog algorithm
        '''
        dt = self.delta_time
        dt2 = dt * dt

        # Update position
        x, y = self.x, self.y
        x += self.vx * dt + 0.5 * self.ax * dt2
        y += self.vy * dt + 0.5 * self.ay * dt2

        # Update velocity using mean acceleration
        ax, ay = acc(x, y)
        self.vx += 0.5 * (self.ax + ax) * dt
        self.vy += 0.5 * (self.ay + ay) * dt

        self.ax, self.ay = ax, ay
        self.x, self.y = x, y
        lx, ly = self.points[-1]
        if abs(x - lx) > 1 or abs(y - ly) > 1:
            self.points.append((x, y))
        return x, y

# Radius and velocity of a circular orbit
radius = 100.0
print('Semi-major axis:', radius)

v = speed(radius)
period = radius ** 1.5
print('Period:', period)

steps = 10000
delta_time = period / steps
istep = steps // 100

body = Body(radius,0, 0,v, delta_time)

print("step radius x y")
for i in range(1, steps + 2):
   x, y = body.update()
   if i % istep == 0: print(i, (x*x + y*y)**0.5, x, y)

# Plot
P = points(body.points, size=1, color="blue")
P.show(aspect_ratio=1)

Eso es en su mayoría código Python simple, aparte de las últimas líneas que usan Sage para trazar (a través de matplotlib).

Aquí está la salida. Como puede ver, el error máximo se produce en el punto medio.

Semi-major axis: 100.0
Period: 1000.0
step radius x y
100 100.00000001947542 99.8026728363567 6.279052365941193
200 100.00000007782481 99.2114701058467 12.533324180026883
300 100.00000017481786 98.22872501631164 18.738132688008204
400 100.00000031007181 96.85831601489033 24.86899034488533
500 100.0000004830529 95.10565148152806 30.901701456144753
600 100.00000069307836 92.97764838452808 36.812457667191076
700 100.0000009393194 90.48270498238338 42.577931924118325
800 100.00000122080417 87.63066767962191 48.17537053499955
900 100.00000153642179 84.43279216747251 53.582682968369106
1000 100.00000188492673 80.90169900271201 58.7785290345032
1100 100.0000022649435 77.05132380000519 63.74240310543257
1200 100.00000267497244 72.89686223430738 68.45471504130965
1300 100.00000311339528 68.4547100703815 72.89686750374857
1400 100.0000035784818 63.74239845610735 77.05132935101909
1500 100.00000406839649 58.77852473495271 80.90170482543533
1600 100.00000458120591 53.582679050660076 84.43279825988925
1700 100.00000511488614 48.17536703380712 87.63067404816206
1800 100.00000566733111 42.577928875365075 90.48271164233866
1900 100.00000623636048 36.812455106635035 92.97765536027522
2000 100.0000068197286 30.90169941794281 95.10565880655513
2100 100.00000741513318 24.86898886015765 96.85832373162384
2200 100.00000802022444 18.73813178343009 98.22873317574528
2300 100.00000863261437 12.53332387647376 99.21147876697717
2400 100.00000924988619 6.279052677211943 99.8026820654341
2500 100.00000986960366 9.316423505156823e-07 100.00000986960366
2600 100.00001048932121 -6.279050817680861 99.80268342431145
2700 100.00001110659296 -12.533322029108925 99.21148149399528
2800 100.00001171898288 -18.738129959342203 98.22873728916358
2900 100.00001232407413 -24.868987074894257 96.8583292581682
3000 100.00001291947869 -30.901697693138058 95.10566578063481
3100 100.00001350284674 -36.81245347153712 92.97766382296042
3200 100.00001407187602 -42.577927368209636 90.48272164011684
3300 100.0000146243209 -48.17536570300814 87.63068563149452
3400 100.00001515800102 -53.582677955791624 84.43281148159379
3500 100.00001567081028 -58.7785239475003 80.9017197389242
3600 100.00001616072485 -63.74239805996919 77.05134600842264
3700 100.00001662581123 -68.45471016210631 72.89688595395924
3800 100.0000170642339 -72.8968629230579 68.45473532798135
3900 100.00001747426262 -77.0513252072431 63.742425264971644
4000 100.00001785427928 -80.90170126160325 58.77855309408762
4100 100.00001820278406 -84.43279542201651 53.58270894403392
4200 100.00001851840162 -87.6306720835148 48.17539842982273
4300 100.00001879988628 -90.48271069762764 42.577961726543805
4400 100.00001904612719 -92.97765557980723 36.81248934952026
4500 100.00001925615258 -95.10566033036771 30.9017349732278
4600 100.00001942913354 -96.85832669363184 24.869025633039456
4700 100.00001956438729 -98.22873770192699 18.738169664187293
4800 100.00001966138026 -99.21148497362778 12.533362741262588
4900 100.00001971972958 -99.80269005751201 6.279092389071799
5000 100.00001973920497 -100.00001973919642 4.13416988284121e-05
5100 100.00001971972962 -99.80269524924384 -6.279009868830847
5200 100.00001966138022 -99.21149533660193 -12.533280709847839
5300 100.0000195643871 -98.2287531952455 -18.738088445339077
5400 100.00001942913315 -96.85834725614961 -24.868945547291286
5500 100.00001925615217 -95.10568588093415 -30.901656336641377
5600 100.00001904612671 -92.97768601758592 -36.812412472438055
5700 100.00001879988581 -90.48274590249471 -42.577886912364384
5800 100.00001851840106 -87.6307119165327 -48.17532597380332
5900 100.00001820278348 -84.43283972598266 -53.582639132125166
6000 100.00001785427867 -80.90174986167018 -58.778486201805165
6100 100.00001747426198 -77.05137791160872 -63.74236155630852
6200 100.00001706423315 -72.89691952372205 -68.45467505436653
6300 100.00001662581037 -68.4547704356921 -72.89682935326505
6400 100.00001616072397 -63.7424617686046 -77.05129330402596
6500 100.00001567080935 -58.77859083975651 -80.9016711388251
6600 100.00001515800002 -53.582747767675556 -84.4327671775946
6700 100.00001462431985 -48.175438159004266 -87.63064579844279
6800 100.0000140718749 -42.578002182367385 -90.4826864352153
6900 100.00001350284549 -36.81253034859931 -92.97763338514665
7000 100.00001291947744 -30.90177632970622 -95.10564023003296
7100 100.00001232407287 -24.86906716062583 -96.8583086956148
7200 100.00001171898157 -18.738211178175494 -98.22872179580955
7300 100.0000111065917 -12.533404060510405 -99.21147113098588
7400 100.00001048931996 -6.279133337910164 -99.80267823254468
7500 100.00000986960247 -8.175174521653228e-05 -100.00000986956906
7600 100.00000924988495 6.278970156979581 -99.8026872571317
7700 100.00000863261317 12.53324184506617 -99.21148912991777
7800 100.0000080202232 18.73805056458767 -98.228748669031
7900 100.00000741513185 24.868908774414006 -96.85834429410973
8000 100.00000681972718 30.90162078135972 -95.10568435709045
8100 100.00000623635901 36.81237822955507 -92.97768579802373
8200 100.00000566732959 42.577854061186905 -90.48274684717644
8300 100.0000051148845 48.17529457778811 -87.63071388115166
8400 100.00000458120421 53.582609238750955 -84.43284256382812
8500 100.00000406839483 58.7784578426692 -80.90175342547614
8600 100.0000035784801 63.74233474744264 -77.05138205535967
8700 100.00000311339362 68.45464979676468 -72.89692410438892
8800 100.00000267497069 72.89680563361073 -68.45477531487283
8900 100.00000226494173 77.05127109560578 -63.74246681404648
9000 100.00000188492496 80.90165040261006 -58.77859592673901
9100 100.00000153642003 84.4327478634704 -53.582752780233704
9200 100.00000122080235 87.63062784656721 -48.17544299097739
9300 100.00000093931756 90.48266977747897 -42.578006738258736
9400 100.00000069307654 92.97761794671167 -36.812534544236804
9500 100.000000483051 95.10562593092368 -30.901780092697216
9600 100.00000031006994 96.85829545233462 -24.869070430601987
9700 100.0000001748159 98.22870952295547 -18.73821390682732
9800 100.00000007782285 99.21145974283534 -12.533406211414873
9900 100.00000001947349 99.80266764458823 -6.279134886157624
10000 99.99999999999805 99.99999999996388 -8.268337527927994e-05

Aquí hay una versión en vivo que se ejecuta en el servidor SageMathCell.

Como menciona Wikipedia, existen versiones de orden superior de Leapfrog, descubiertas de forma independiente en las décadas de 1980 y 1990 por varias personas (Forest & Ruth, Yoshida y Candy & Rozmus). La versión de cuarto orden es ciertamente superior a Leapfrog o Verlet simple, pero las mejoras frente al mayor número de cálculos disminuyen con los órdenes más altos, por lo que no tiene mucho sentido ir más allá del cuarto u octavo orden.

Aquí hay una versión del código anterior que usa coeficientes de Yoshida de cuarto orden.

''' Simple orbit sim
    4th order Leapfrog with Yoshida coefficients
    Written by PM 2Ring 2017.05.20
    Added Yoshida 2022.02.12
'''

from math import pi, radians, cos, sin

# Gravitational parameter of central body in "natural" units,
# i.e. where R**3 == T**2
mu = 4.0 * pi * pi

# Speed of a body in a circular orbit of radius r
def speed(r):
    return (mu / r) ** 0.5

# Acceleration vector due to central gravity at (x, y)
def acc(x, y):
    a = -mu / (x * x + y * y) ** 1.5
    return a * x, a * y

# Yoshida 4th order coefficients
def yosh4():
    q = 2 ** (1 / 3)
    w1 = 1 / (2 - q)
    w0 = -q * w1
    c1 = w1 / 2
    c2 = (w0 + w1) / 2
    return (c1, c2, c2, c1), (w1, w0, w1)

class Body(object):
    ''' An orbiting body '''
    def __init__(self, x, y, vx, vy, delta_time):
        # Current position
        self.x, self.y = x, y

        # Current velocity
        self.vx, self.vy = vx, vy

        # The Yoshida coefficients, multiplied by the time step
        yc, yd = yosh4()
        self.yc = [k * delta_time for k in yc]
        self.yd = [k * delta_time for k in yd]

        self.points = [(x, y)]

    def update(self):
        ''' Update the body's orbit parameters using
            the 4th order synchronised leapfrog algorithm
        '''
        x, y = self.x, self.y
        vx, vy = self.vx, self.vy

        # 4th order Leapfrog integration
        for c, d in zip(self.yc, self.yd):
            x += vx * c
            y += vy * c
            ax, ay = acc(x, y)
            vx += ax * d
            vy += ay * d
        c = self.yc[-1]
        x += vx * c
        y += vy * c

        self.vx, self.vy = vx, vy
        self.x, self.y = x, y

        lx, ly = self.points[-1]
        if abs(x - lx) > 1 or abs(y - ly) > 1:
            self.points.append((x, y))
        return x, y

# Radius and velocity of a circular orbit
radius = 100.0
print('Semi-major axis:', radius)

v = speed(radius)
period = radius ** 1.5
print('Period:', period)

steps = 1000
delta_time = period / steps
print('Step:', delta_time)
istep = steps // 100

# Initial direction. Use 0 for a circular orbit
th = radians(0)
vx, vy = v * sin(th), v * cos(th)

body = Body(radius,0, vx,vy, delta_time)
for i in range(1, steps + 2):
   x, y = body.update()
   if i % istep == 0: print(i, (x*x + y*y)**0.5, x, y)

P = points(body.points, size=1, color="blue")
P.show(aspect_ratio=1)

Con solo 1000 pasos por órbita, el error máximo en el radio se ha reducido a 1,1 partes por billón. Aquí está el enlace de SageMathCell .

Aquí hay una gráfica de 10 períodos de una órbita excéntrica creada usando ese código. Utiliza el mismo paso de tiempo de 1.0, con un ángulo inicial de 60° a la vertical.

El error de radio ha saltado hasta 5 partes en diez millones, por órbita, aunque eso no es visible a esta escala.

Si aumentamos el paso de tiempo a 2,0, el error pasa a ser de 6 partes por millón, y empieza a ser evidente en el gráfico.

Entonces, hay una serie de cosas que pueden suceder aquí bajo el capó proverbial que pueden causar que surja este tipo de error.

En primer lugar, el error de redondeo estará en el centro de todos sus problemas. Digamos que la precisión de su máquina es de aproximadamente$10^{-16}$, que es si está utilizando una representación de punto flotante de 64 bits (esto es bastante estándar en todos los ámbitos, tiene que hacer todo lo posible para cambiar esto). Digamos que estás resolviendo la ODE$\frac{d^2 r}{dt^2}=\frac{GM}{r^2}$y lo divides en componentes individuales, etc. Tu distancia es 10 000 000 o$10^7$, que luego se eleva al cuadrado para que trabaje con$10^{14}$en la parte inferior de esa ecuación,$1.9 \times 10^{27}$kg en la parte superior, y todo eso multiplicado por$6.67 \times 10^{-11}$. Como puede adivinar, estos números enormemente diferentes pueden comenzar a crear problemas con errores de redondeo (si esto no está claro, hágamelo saber y puedo agregar un ejemplo más concreto). Comenzará a perder higos de firma efectivos, y ese error puede conducir a cosas como las que está viendo.

Además, no todos los solucionadores son iguales. Los solucionadores simplécticos conservan energía, lo cual es genial, pero no hay comida gratis en física computacional; tienen un costo. La naturaleza exacta de ese costo no está clara para mí, pero siempre me han aconsejado evitarlos a menos que la conservación de energía a nivel granular sea importante (mi investigación actual es con objetos relativistas de N-cuerpos y mi asesor nos ha desaconsejado integradores simplécticos , nos dio una razón, pero no puedo recordar bien cuál era, lo revisaré)

Más allá de eso, no todos los solucionadores se crean por igual. Si puede encontrar un solucionador implícito, será más lento pero más preciso y, por lo general, cuanto mayor sea el orden, mejor. Algunos solucionadores tienen capacidades para lidiar con la rigidez y otras cosas, mientras que otros no; comienza a convertirse más en un arte oscuro que en una ciencia en cierto nivel porque las razones por las que algunos funcionan más que otros se vuelven más y más turbias.

tl; dr, pruebe con varios solucionadores diferentes y escale sus unidades más cerca de la unidad (por ejemplo, use la masa como 1 para reimprimir 1 masa de Júpiter, use algo cercano a la unidad para la distancia y cambie la escala de G en consecuencia. Métalo hasta que G esté dentro de un par de órdenes de magnitud a la unidad también)

Lo que sigue es demasiado largo para ser un comentario. Además, dado que los comentarios se pueden eliminar o mover al chat, estoy respondiendo esto.

Todo integrador numérico basado en el avance del estado cartesiano de un paso de tiempo al siguiente tiene un tamaño de paso óptimo que depende de

• La naturaleza del problema en cuestión,
• Las unidades utilizadas para representar el problema en cuestión,
• Los mecanismos utilizados para representar los números reales, y
• Los mecanismos utilizados para avanzar de estado de un paso de tiempo al siguiente.

Veré los dos últimos elementos, con un enfoque en las técnicas de integración que avanzan la posición cartesiana y la velocidad en función de la aceleración.

La integración numérica de la posición de un vehículo en el espacio es inherentemente un problema de ecuación diferencial ordinaria (ODE) de segundo orden. Algunas técnicas de integración numérica (p. ej., RK4 canónica) solo pueden resolver EDO de primer orden. Una EDO de segundo orden se puede resolver usando un solucionador de EDO de primer orden agregando ecuaciones auxiliares, pero al hacerlo, se arroja la geometría por la ventana. (Que esta es una ODE de segundo orden es geometría). Además de arrojar geometría, técnicas como RK4 también arrojan las leyes de conservación por la ventana.

Cada técnica de integración numérica que avanza la velocidad cartesiana y la posición en función de la aceleración se enfrenta a dos desafíos básicos. Si el paso de tiempo es extremadamente pequeño, la velocidad y/o la posición se congelarán. El problema con los pasos de tiempo extremadamente pequeños es que$1+10^{-16}$es exactamente igual a uno en aritmética de doble precisión. Del otro lado del tamaño del paso de tiempo, surgen errores inherentes a la propia técnica si el paso de tiempo se hace extremadamente grande. Ningún integrador numérico funcionará bien con un tamaño de paso de diez órbitas.

Para un problema dado, un conjunto de unidades dado, una técnica de representación dada de los reales y una técnica de integración dada, existe un tamaño de paso óptimo que se encuentra entre demasiado pequeño y demasiado grande. Este tamaño de paso óptimo bien podría no ser constante. Por ejemplo, para una órbita elíptica, es mejor dar pasos más pequeños cerca del periapsis y pasos más grandes cerca del apoapsis.

Este tamaño de paso óptimo es algo fácil de encontrar para órbitas circulares, ya que hace que el tamaño de paso óptimo sea constante. Lo que se verá es una curva de error de bañera. Idealmente, los integradores deberían nadar cerca del extremo profundo de esta bañera.