使用 Python 为 3 体问题创建蛙跳算法的问题

Question

我正在尝试使用 leapfrog 算法为三体问题编写代码。我使用 Piet Hut 和 Jun Makino 的 "Moving Stars Around" 作为指南。

指南中的代码是用 C 语言编写的，但在进行实验之前，我尝试使用 Python 作为起点来遵循确切的工作流程。

以下是我尝试遵循 section 5.1 中的代码。

import numpy as np

N = 3       #number of bodies
m = 1       #mass
dt = 0.01   #timestep
t_end = 10  #duration

r = [] 
v = []
a = [[0, 0, 0] for i in range(N)]

for i in range(N):
    phi = i * 2 * np.pi/3
    r.append([np.cos(phi), np.sin(phi), 0])

for i in range(N):
    for j in range(i+1, N):
        rji = []
        for k in range(3):
            rji.append(r[j][k] - r[i][k])
        r2 = 0
        for k in range(3):
            r2 += rji[k]**2
        r3 = r2 * np.sqrt(r2)
        for k in range(3):
            a[i][k] += m * rji[k] / r3
            a[j][k] -= m * rji[k] / r3


v_abs = np.sqrt(-a[0][0])
for i in range(N):
    phi = i * 2 * np.pi/3
    v.append([-v_abs * np.sin(phi),
              v_abs * np.cos(phi), 0])


ekin = 0
epot = 0
for i in range(N):

    for j in range(i+1, N):
        rji = [0, 0, 0]
        for k in range(3):
            rji[k] = r[j][k] - r[i][k]
        r2 = 0
        for k in range(3):
            r2 += rji[k]**2
        d = np.sqrt(r2)
        epot -= m**2 / d

    for k in range(3):
        ekin += 0.5 * m * v[i][k]**2

e_in = ekin + epot
print('Initial total energy E_in = ', e_in)


dt_out = 0.01
t_out = dt_out

for t in np.arange(0, t_end, dt):

    for i in range(N):
        for k in range(3):
            v[i][k] += a[i][k] * dt/2
        for k in range(3):
            r[i][k] += v[i][k] * dt

    for i in range(N):
        for k in range(3):
            a[i][k] = 0

    for i in range(N):
        for j in range(i+1, N):
            rji = []
        for k in range(3):
            rji.append(r[j][k] - r[i][k])
        r2 = 0
        for k in range(3):
            r2 += rji[k]**2
        r3 = r2 * np.sqrt(r2)
        for k in range(3):
            a[i][k] += m * rji[k] / r3
            a[j][k] -= m * rji[k] / r3

    for i in range(N):
        for k in range(3):
            v[i][k] += a[i][k] * dt/2
        '''
        if t >= t_out:
            for i in range(N):
                print(r[i][k], ' ')
            for k in range(N):
                print(v[i][k], ' ')
        '''
        t_out += dt_out


epot = 0
ekin = 0
for i in range(N):

    for j in range(i+1, N):
        rji = [0, 0, 0]
        for k in range(3):
            rji[k] = r[j][k] - r[i][k]
        r2 = 0
        for k in range(3):
            r2 += rji[k]**2
        d = np.sqrt(r2)
        epot -= m**2 / d

    for k in range(3):
        ekin += 0.5 * m * v[i][k]**2

e_out = ekin + epot
print('Final total energy E_out = ', e_out)
print('absolute energy error: E_out - E_in = ', e_out - e_in)
print('relative energy error: (E_out - E_in)/E_in = ', (e_out - e_in)/e_in)

我定义了时间步长 dt = 0.01 和持续时间 t_end = 10，而不是提示输入。在 section 5.4 中，结果应该是：

|gravity> g++ -o leapfrog2 leapfrog2.C
|gravity> leapfrog2 > leapfrog2_0.01_10.out
Please provide a value for the time step
0.01
and for the duration of the run
10
Initial total energy E_in = -0.866025
Final total energy E_out = -0.866025
absolute energy error: E_out - E_in = 2.72254e-10
relative energy error: (E_out - E_in) / E_in = -3.14372e-10

连同圆形图。但是，我的代码的结果有所不同：

Initial total energy E_in =  -0.8660254037844386
Final total energy E_out =  -0.39922101519288833
absolute energy error: E_out - E_in =  0.46680438859155027
relative energy error: (E_out - E_in)/E_in =  -0.5390192788244604

当然，在我绘制结果后，它们不会绕圈子。

我想知道我在翻译代码时是否犯了错误。任何帮助，将不胜感激！

Answer 1

欢迎来到堆栈溢出！

首先，该错误是一个典型的 python 问题：您的部分代码缩进不正确。具体来说：

for i in range(N):
    for j in range(i+1, N):
        rji = []
    for k in range(3):
        rji.append(r[j][k] - r[i][k])
    r2 = 0
    for k in range(3):
        r2 += rji[k]**2
    r3 = r2 * np.sqrt(r2)
    for k in range(3):
        a[i][k] += m * rji[k] / r3
        a[j][k] -= m * rji[k] / r3

应该是：

for i in range(N):
    for j in range(i+1, N):
        rji = []
        for k in range(3):
            rji.append(r[j][k] - r[i][k])
        r2 = 0
        for k in range(3):
            r2 += rji[k]**2
        r3 = r2 * np.sqrt(r2)
        for k in range(3):
            a[i][k] += m * rji[k] / r3
            a[j][k] -= m * rji[k] / r3

我给你一个建议：如果你想在阅读本书的同时学习 python，请尝试编写一个可以完成其工作的版本（就像我们在这里讨论的那个），然后制作努力使其更加地道。通过使用 numpy，您可以删除空间维度上的大部分（如果不是全部）循环（至少！）。

使用 Python 为 3 体问题创建蛙跳算法的问题

Problem in creating leapfrog algorithm for 3-body problem using Python

python

simulation

orbital-mechanics