本文介绍了用Python创建三体问题的越级算法存在的问题的处理方法，对大家解决问题具有一定的参考价值，需要的朋友们下面随着小编来一起学习吧！
问题描述
我正在尝试用LeapFrog算法编写三体问题的代码。我正在使用Piet Hut和Jun Makino的《Moving Stars Anover》作为指南。
本指南中的代码是用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
我想知道我在翻译代码时是否犯了错误。如有任何帮助，我们将不胜感激！
推荐答案
欢迎使用堆栈溢出！
首先，该错误是典型的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
让我给您一个建议：如果您正在学习书中的内容，请尝试编写一个能够完成其工作的版本(就像我们在这里讨论的那个版本)，然后努力使其更加地道。通过使用NumPy，您可以删除空间维度上的大部分(如果不是全部)循环(至少！)。
这篇关于用Python创建三体问题的越级算法存在的问题的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持编程学习网！
本站部分内容来源互联网,如果有图片或者内容侵犯您的权益请联系我们删除！