Python N 体模拟代码给出错误答案
Python N-body simulation code giving incorrect answer
我已经编写了一些 python 代码来使用 Euler 方法解决 N 体问题。代码 运行 没有问题，似乎给出了合理的答案（例如，如果有两个粒子，则开始相互移动）。然而，当我 运行 这个模拟经过大量迭代时，我看到粒子（假设我 运行 它有两个粒子）相互经过（我不考虑碰撞）并继续前进方向无限期。这违反了能量守恒，所以我的代码中肯定有缺陷，但我找不到它。谁能找到它并解释我的错误。
谢谢。
感谢@samgak 指出我更新了两次粒子。我现在已经解决了这个问题，但问题仍然存在。我还复制了我在 运行 此模拟中得到的输出，其中两个静止粒子位于 (0,0) 和 (1,0)，时间步长为 1 秒，迭代次数为 100000 次：
质量为 1 位置为 [234.8268420043934, 0.0] 速度为 [0.011249111128594091, 0.0] 的粒子
质量为 1 位置为 [-233.82684200439311, 0.0] 速度为 [-0.011249111128594091, 0.0]
的粒子
还要感谢@PM2Ring 指出我可以进行的一些优化以及使用欧拉方法的风险。
代码：
import math
class Particle:
"""
Class to represent a single particle
"""
def __init__(self,mass,position,velocity):
"""
Initialize the particle
"""
self.G = 6.67408*10**-11 #fixed throughout the simulation
self.time_interval = 10**0 #fixed throughout the simulation, gives the interval between updates
self.mass = mass
self.position = position #should be a list
self.velocity = velocity #should be a list
self.updated_position = position
self.updated_velocity = velocity
def __str__(self):
"""
String representation of particle
"""
return "Particle with mass: " + str(self.mass) + " and position: " + str(self.position) + " and velocity: " + str(self.velocity)
def get_mass(self):
"""
Returns the mass of the particle
"""
return self.mass
def get_position(self):
"""
returns the position of the particle
"""
return self.position
def get_velocity(self):
"""
returns the velocity of the particle
"""
return self.velocity
def get_updated_position(self):
"""
calculates the future position of the particle
"""
for i in range(len(self.position)):
self.updated_position[i] = self.updated_position[i] + self.time_interval*self.velocity[i]
def update_position(self):
"""
updates the position of the particle
"""
self.position = self.updated_position.copy()
def get_distance(self,other_particle):
"""
returns the distance between the particle and another given particle
"""
tot = 0
other = other_particle.get_position()
for i in range(len(self.position)):
tot += (self.position[i]-other[i])**2
return math.sqrt(tot)
def get_updated_velocity(self,other_particle):
"""
updates the future velocity of the particle due to the acceleration
by another particle
"""
distance_vector = []
other = other_particle.get_position()
for i in range(len(self.position)):
distance_vector.append(self.position[i]-other[i])
distance_squared = 0
for item in distance_vector:
distance_squared += item**2
distance = math.sqrt(distance_squared)
force = -self.G*self.mass*other_particle.get_mass()/(distance_squared)
for i in range(len(self.velocity)):
self.updated_velocity[i] = self.updated_velocity[i]+self.time_interval*force*(distance_vector[i])/(self.mass*(distance))
def update_velocity(self):
"""
updates the velocity of the particle
"""
self.velocity = self.updated_velocity.copy()
def update_particles(particle_list):
"""
updates the position of all the particles
"""
for i in range(len(particle_list)):
for j in range(i+1,len(particle_list)):
particle_list[i].get_updated_velocity(particle_list[j])
particle_list[j].get_updated_velocity(particle_list[i])
for i in range(len(particle_list)):
particle_list[i].update_velocity()
particle_list[i].get_updated_position()
for i in range(len(particle_list)):
particle_list[i].update_position()
#the list of particles
partList = [Particle(1,[0,0],[0,0]),Particle(1,[1,0],[0,0])]
#how many iterations I perform
for i in range(100000):
update_particles(partList)
#prints out the final position of all the particles
for item in partList:
print(item)
-------------------------------------------- ---------------------------------------------- ---------------------------------------------- ---------------------------------------------- ---------------------------------------------- ------------------进一步编辑：
我决定实施 Leapfrog 方法，我已经开发了一些代码，再次 运行s 并且似乎运行良好（至少在命令行中）。但是，当我添加绘图功能并对其进行分析时，似乎出现了另一个问题。系统似乎再次走得太远，能量再次无限制地增加。我附上了我得到的输出图片来展示问题。如果我再次只有两个质量相等的粒子，它们将再次相互通过并继续相互远离而不会停止。因此，我的代码中一定存在我没有发现的错误。
如果有人能提供帮助，将不胜感激。
我的代码：
import math
import matplotlib.pyplot as plt
class Particle:
"""
Represents a single particle
"""
def __init__(self,mass,position,velocity):
"""
Initialize the particle
"""
self.G = 6.67408*10**-11
self.time_step = 10**2
self.mass = mass
self.dimensions = len(position)
self.position = position
self.velocity = velocity
self.acceleration = [0 for i in range(len(position))]
self.next_position = position
self.next_velocity = velocity
self.next_acceleration = [0 for i in range(len(position))]
def __str__(self):
"""
A string representation of the particle
"""
return "A Particle with mass: " + str(self.mass) + " and position: " + str(self.position) + " and velocity:" + str(self.velocity)
def get_mass(self):
return self.mass
def get_position(self):
return self.position
def get_velocity(self):
return self.velocity
def get_acceleration(self):
return self.acceleration
def get_next_position(self):
return self.next_position
def put_next_position(self):
for i in range(self.dimensions):
self.next_position[i] = self.position[i] + self.time_step*self.velocity[i]+0.5*self.time_step**2*self.acceleration[i]
def put_next_velocity(self):
for i in range(self.dimensions):
self.next_velocity[i] = self.velocity[i] + 0.5*self.time_step*(self.acceleration[i]+self.next_acceleration[i])
def update_position(self):
self.position = self.next_position.copy()
def update_velocity(self):
self.velocity = self.next_velocity.copy()
def update_acceleration(self):
self.acceleration = self.next_acceleration.copy()
def reset_acceleration(self):
self.acceleration = [0 for i in range(self.dimensions)]
def reset_future_acceleration(self):
self.next_acceleration = [0 for i in range(self.dimensions)]
def calculate_acceleration(self,other_particle):
"""
Increments the acceleration of the particle due to the force from
a single other particle
"""
distances = []
other = other_particle.get_position()
distance_squared = 0
for i in range(self.dimensions):
distance_squared += (self.position[i]-other[i])**2
distances.append(self.position[i]-other[i])
distance = math.sqrt(distance_squared)
force = -self.G*self.mass*other_particle.get_mass()/distance_squared
acc = []
for i in range(self.dimensions):
acc.append(force*distances[i]/(distance*self.mass))
for i in range(self.dimensions):
self.acceleration[i] += acc[i]
def calculate_future_acceleration(self,other_particle):
"""
Increments the future acceleration of the particle due to the force from
a single other particle
"""
distances = []
other = other_particle.get_next_position()
distance_squared = 0
for i in range(self.dimensions):
distance_squared += (self.next_position[i]-other[i])**2
distances.append(self.next_position[i]-other[i])
distance = math.sqrt(distance_squared)
force = -self.G*self.mass*other_particle.get_mass()/distance_squared
acc = []
for i in range(self.dimensions):
acc.append(force*distances[i]/(distance*self.mass))
for i in range(self.dimensions):
self.next_acceleration[i] += acc[i]
def update_all(particleList):
for i in range(len(particleList)):
particleList[i].reset_acceleration()
for j in range(len(particleList)):
if i != j:
particleList[i].calculate_acceleration(particleList[j])
for i in range(len(particleList)):
particleList[i].put_next_position()
for i in range(len(particleList)):
particleList[i].reset_future_acceleration()
for j in range(len(particleList)):
if i != j:
particleList[i].calculate_future_acceleration(particleList[j])
for i in range(len(particleList)):
particleList[i].put_next_velocity()
for i in range(len(particleList)):
particleList[i].update_position()
particleList[i].update_velocity()
partList = [Particle(1,[0,0],[0,0]),Particle(1,[1,0],[0,0])]
Alist = [[],[]]
Blist = [[],[]]
for i in range(10000):
Alist[0].append(partList[0].get_position()[0])
Alist[1].append(partList[0].get_position()[1])
Blist[0].append(partList[1].get_position()[0])
Blist[1].append(partList[1].get_position()[1])
update_all(partList)
plt.scatter(Alist[0],Alist[1],color="r")
plt.scatter(Blist[0],Blist[1],color="b")
plt.grid()
plt.show()
for item in partList:
print(item)
谁能告诉我我的代码哪里出错了。
代码中的主要问题是它使用了 Euler 方法，随着迭代次数的增加，该方法非常不准确（仅 O(h) 与其他方法相比，可以是 O(h^4) 甚至更好).要解决这个问题，需要对代码进行根本性的重组，所以我会说这段代码对于 N 体模拟来说并不是很准确（它只适用于 2 个粒子，随着我添加的越来越多，错误只会增加）。
感谢@samgak 和@PM2Ring 帮助我消除了错误并优化了我的代码，但总的来说这段代码无法使用...
编辑：我已经从头开始实施了评论中提到的越级方法，并且发现它工作得很好。它非常易于理解和实施，而且效果很好！
进一步编辑：我以为我有 leapfrog 方法工作。结果发现其中还有另一个错误，我是在添加 GUI 功能时才看到的。
我已经编写了一些 python 代码来使用 Euler 方法解决 N 体问题。代码 运行 没有问题，似乎给出了合理的答案（例如，如果有两个粒子，则开始相互移动）。然而，当我 运行 这个模拟经过大量迭代时，我看到粒子（假设我 运行 它有两个粒子）相互经过（我不考虑碰撞）并继续前进方向无限期。这违反了能量守恒，所以我的代码中肯定有缺陷，但我找不到它。谁能找到它并解释我的错误。
谢谢。
感谢@samgak 指出我更新了两次粒子。我现在已经解决了这个问题，但问题仍然存在。我还复制了我在 运行 此模拟中得到的输出，其中两个静止粒子位于 (0,0) 和 (1,0)，时间步长为 1 秒，迭代次数为 100000 次：
质量为 1 位置为 [234.8268420043934, 0.0] 速度为 [0.011249111128594091, 0.0] 的粒子
质量为 1 位置为 [-233.82684200439311, 0.0] 速度为 [-0.011249111128594091, 0.0]
的粒子还要感谢@PM2Ring 指出我可以进行的一些优化以及使用欧拉方法的风险。
代码：
import math
class Particle:
"""
Class to represent a single particle
"""
def __init__(self,mass,position,velocity):
"""
Initialize the particle
"""
self.G = 6.67408*10**-11 #fixed throughout the simulation
self.time_interval = 10**0 #fixed throughout the simulation, gives the interval between updates
self.mass = mass
self.position = position #should be a list
self.velocity = velocity #should be a list
self.updated_position = position
self.updated_velocity = velocity
def __str__(self):
"""
String representation of particle
"""
return "Particle with mass: " + str(self.mass) + " and position: " + str(self.position) + " and velocity: " + str(self.velocity)
def get_mass(self):
"""
Returns the mass of the particle
"""
return self.mass
def get_position(self):
"""
returns the position of the particle
"""
return self.position
def get_velocity(self):
"""
returns the velocity of the particle
"""
return self.velocity
def get_updated_position(self):
"""
calculates the future position of the particle
"""
for i in range(len(self.position)):
self.updated_position[i] = self.updated_position[i] + self.time_interval*self.velocity[i]
def update_position(self):
"""
updates the position of the particle
"""
self.position = self.updated_position.copy()
def get_distance(self,other_particle):
"""
returns the distance between the particle and another given particle
"""
tot = 0
other = other_particle.get_position()
for i in range(len(self.position)):
tot += (self.position[i]-other[i])**2
return math.sqrt(tot)
def get_updated_velocity(self,other_particle):
"""
updates the future velocity of the particle due to the acceleration
by another particle
"""
distance_vector = []
other = other_particle.get_position()
for i in range(len(self.position)):
distance_vector.append(self.position[i]-other[i])
distance_squared = 0
for item in distance_vector:
distance_squared += item**2
distance = math.sqrt(distance_squared)
force = -self.G*self.mass*other_particle.get_mass()/(distance_squared)
for i in range(len(self.velocity)):
self.updated_velocity[i] = self.updated_velocity[i]+self.time_interval*force*(distance_vector[i])/(self.mass*(distance))
def update_velocity(self):
"""
updates the velocity of the particle
"""
self.velocity = self.updated_velocity.copy()
def update_particles(particle_list):
"""
updates the position of all the particles
"""
for i in range(len(particle_list)):
for j in range(i+1,len(particle_list)):
particle_list[i].get_updated_velocity(particle_list[j])
particle_list[j].get_updated_velocity(particle_list[i])
for i in range(len(particle_list)):
particle_list[i].update_velocity()
particle_list[i].get_updated_position()
for i in range(len(particle_list)):
particle_list[i].update_position()
#the list of particles
partList = [Particle(1,[0,0],[0,0]),Particle(1,[1,0],[0,0])]
#how many iterations I perform
for i in range(100000):
update_particles(partList)
#prints out the final position of all the particles
for item in partList:
print(item)
-------------------------------------------- ---------------------------------------------- ---------------------------------------------- ---------------------------------------------- ---------------------------------------------- ------------------进一步编辑：
我决定实施 Leapfrog 方法，我已经开发了一些代码，再次 运行s 并且似乎运行良好（至少在命令行中）。但是，当我添加绘图功能并对其进行分析时，似乎出现了另一个问题。系统似乎再次走得太远，能量再次无限制地增加。我附上了我得到的输出图片来展示问题。如果我再次只有两个质量相等的粒子，它们将再次相互通过并继续相互远离而不会停止。因此，我的代码中一定存在我没有发现的错误。
如果有人能提供帮助，将不胜感激。
我的代码：
import math
import matplotlib.pyplot as plt
class Particle:
"""
Represents a single particle
"""
def __init__(self,mass,position,velocity):
"""
Initialize the particle
"""
self.G = 6.67408*10**-11
self.time_step = 10**2
self.mass = mass
self.dimensions = len(position)
self.position = position
self.velocity = velocity
self.acceleration = [0 for i in range(len(position))]
self.next_position = position
self.next_velocity = velocity
self.next_acceleration = [0 for i in range(len(position))]
def __str__(self):
"""
A string representation of the particle
"""
return "A Particle with mass: " + str(self.mass) + " and position: " + str(self.position) + " and velocity:" + str(self.velocity)
def get_mass(self):
return self.mass
def get_position(self):
return self.position
def get_velocity(self):
return self.velocity
def get_acceleration(self):
return self.acceleration
def get_next_position(self):
return self.next_position
def put_next_position(self):
for i in range(self.dimensions):
self.next_position[i] = self.position[i] + self.time_step*self.velocity[i]+0.5*self.time_step**2*self.acceleration[i]
def put_next_velocity(self):
for i in range(self.dimensions):
self.next_velocity[i] = self.velocity[i] + 0.5*self.time_step*(self.acceleration[i]+self.next_acceleration[i])
def update_position(self):
self.position = self.next_position.copy()
def update_velocity(self):
self.velocity = self.next_velocity.copy()
def update_acceleration(self):
self.acceleration = self.next_acceleration.copy()
def reset_acceleration(self):
self.acceleration = [0 for i in range(self.dimensions)]
def reset_future_acceleration(self):
self.next_acceleration = [0 for i in range(self.dimensions)]
def calculate_acceleration(self,other_particle):
"""
Increments the acceleration of the particle due to the force from
a single other particle
"""
distances = []
other = other_particle.get_position()
distance_squared = 0
for i in range(self.dimensions):
distance_squared += (self.position[i]-other[i])**2
distances.append(self.position[i]-other[i])
distance = math.sqrt(distance_squared)
force = -self.G*self.mass*other_particle.get_mass()/distance_squared
acc = []
for i in range(self.dimensions):
acc.append(force*distances[i]/(distance*self.mass))
for i in range(self.dimensions):
self.acceleration[i] += acc[i]
def calculate_future_acceleration(self,other_particle):
"""
Increments the future acceleration of the particle due to the force from
a single other particle
"""
distances = []
other = other_particle.get_next_position()
distance_squared = 0
for i in range(self.dimensions):
distance_squared += (self.next_position[i]-other[i])**2
distances.append(self.next_position[i]-other[i])
distance = math.sqrt(distance_squared)
force = -self.G*self.mass*other_particle.get_mass()/distance_squared
acc = []
for i in range(self.dimensions):
acc.append(force*distances[i]/(distance*self.mass))
for i in range(self.dimensions):
self.next_acceleration[i] += acc[i]
def update_all(particleList):
for i in range(len(particleList)):
particleList[i].reset_acceleration()
for j in range(len(particleList)):
if i != j:
particleList[i].calculate_acceleration(particleList[j])
for i in range(len(particleList)):
particleList[i].put_next_position()
for i in range(len(particleList)):
particleList[i].reset_future_acceleration()
for j in range(len(particleList)):
if i != j:
particleList[i].calculate_future_acceleration(particleList[j])
for i in range(len(particleList)):
particleList[i].put_next_velocity()
for i in range(len(particleList)):
particleList[i].update_position()
particleList[i].update_velocity()
partList = [Particle(1,[0,0],[0,0]),Particle(1,[1,0],[0,0])]
Alist = [[],[]]
Blist = [[],[]]
for i in range(10000):
Alist[0].append(partList[0].get_position()[0])
Alist[1].append(partList[0].get_position()[1])
Blist[0].append(partList[1].get_position()[0])
Blist[1].append(partList[1].get_position()[1])
update_all(partList)
plt.scatter(Alist[0],Alist[1],color="r")
plt.scatter(Blist[0],Blist[1],color="b")
plt.grid()
plt.show()
for item in partList:
print(item)
谁能告诉我我的代码哪里出错了。
代码中的主要问题是它使用了 Euler 方法，随着迭代次数的增加，该方法非常不准确（仅 O(h) 与其他方法相比，可以是 O(h^4) 甚至更好).要解决这个问题，需要对代码进行根本性的重组，所以我会说这段代码对于 N 体模拟来说并不是很准确（它只适用于 2 个粒子，随着我添加的越来越多，错误只会增加）。
感谢@samgak 和@PM2Ring 帮助我消除了错误并优化了我的代码，但总的来说这段代码无法使用...
编辑：我已经从头开始实施了评论中提到的越级方法，并且发现它工作得很好。它非常易于理解和实施，而且效果很好！
进一步编辑：我以为我有 leapfrog 方法工作。结果发现其中还有另一个错误，我是在添加 GUI 功能时才看到的。