我已经编写了一些
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指出我可以进行的一些优化以及使用Euler方法的危险.
码：
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方法,并且我已经开发了一些代码,它们再次运行并且似乎运行良好(至少在命令行中).然而,当我添加绘图功能并进行分析时,似乎还有另一个问题.系统似乎再次走得太远,能量再次无限制地增加.我附上了一张我输出的图片来展示问题.如果我再次只有两个质量相等的粒子,它们会再次相互通过而不停地相互离开.因此,我的代码中肯定存在一个我找不到的错误.
如果有人可以帮助它将非常感激.
我的代码：
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 ^ 4)甚至更好的方法相比,只有O(h).要解决这个问题,需要对代码进行基本的重组,所以我会说这个代码对于N体仿真来说并不是很准确(它可以播放2个粒子,因为我越来越多地添加错误只会增加).
感谢@samgak和@ PM2Ring帮助我删除错误并优化我的代码,但整体而言这段代码无法使用…
编辑：我已经实现了从头开始的评论中提到的跳过方法,并发现它完美地工作.它理解和实现起来非常简单,也很有效！
进一步编辑：我以为我有越级方法工作.事实证明,当我添加GUI功能时,我才会看到另一个错误.
