N body simulation 自己也有点兴趣所以试着看了一下
import numpy as np
#Numpy是科学计算常用的套件,运算量大时会比原生list快
particle = np.random.standard_normal((nParticles, 3))
#三维阵列, 存放 nParticles个球的x,y,z座标
#ex. nParticles = 100
#这边是用随机产生的座标
particlev = np.zeros_like(particle)
#存放每个球的速度的阵列,形状跟particle一样
#初始值为零
def nbody_np(particle, particlev): # NumPy arrays as input
'''
从初始条件开始,计算每个球受的重力
时间轴每次增加0.01,总共跑五次
'''
t0 = time.time(); nSteps = 5; dt = 0.01
particle_result = np.zeros(nParticles,3,nSteps+1)
#原本我贴的例子只告诉你花多少时间算
#这里用particle_result记录每个时间点的球的位置
particle_log = np.copy(particle) #从起始点开始
for step in range(1, nSteps + 1, 1):
Fp = np.zeros((nParticles, 3)) #每个球受的重力
for i in range(nParticles):
dp = particle - particle[i] #计算第i个球到每个球之间的距离
#也就是dx, dy, dz
drSquared = np.sum(dp ** 2, axis=1)
#(dx)^2 + (dy)^2 +(dz)^2
h = drSquared * np.sqrt(drSquared) #原程式错了
#不是相加而是相乘
#数值运算有分辨率问题
#这边设定重力不能小于10^-10
#至于重力为什么是1/(r^(3/2)),
#是假设重力常数跟所有球的质量都是1
#3/2次方的原因是来自向量方程式
drPowerN32 = 1. / np.maximum(h, 1E-10)
Fp += -(dp.T * drPowerN32).T
#最后每个球会受到所有球的重力所以要相加
particlev += dt * Fp
#球的速度在t+dt后就是原速度+dt*重力
particle_log += particlev * dt
particle_result[:,:,step] += particle_log
#球的新位置是 原位置+球的速度*dt
return particle_result
python科学运算原则上能向量化就向量化
除非是用c之类的编译语言写否则尽量避免for loop会比较好
向量化之后应该没有帮每个球取名的必要吧?
至于Nbody simulation如果球不是一个质点的情况
你的运动方程式会更复杂喔,不知道高中物理现在是有没有这么难 XD
※ 引述《lefan (红气球雯雯)》之铭言：
: 谢谢Neisseria大介绍globals函数让我解决了帮球自动取名的问题
: 但又碰上新的问题，
: 我希望在每一个循环中，自动把每个球的位置塞入新的list中，
: 好让我可以每个循环重新计算球与球间的距离。
: 若不用循环我会这样写:
: b_new_pos_list = []
: b_new_pos_list.append(ball_0.pos)
: b_new_pos_list.append(ball_1.pos)
: b_new_pos_list.append(ball_2.pos)
: b_new_pos_list.append(ball_3.pos)
: 相同的，我想利用for loop自动把每个球的位置放入b_new_pos_list中
: 因此我尝试这样写。
: b_new_pos_list=[]
: for N in range(0,4,1):
: b_new_pos_list.append(ball_N.pos)
: 但当然还是不行，因为系统没办法自动判断出ball_N.pos指的就是
: ball_0~3.pos
: 再次感谢。