為什麼這篇numpy append速度鄉民發文收入到精華區：因為在numpy append速度這個討論話題中，有許多相關的文章在討論，這篇最有參考價值！作者painkiller (肚子餓~)看板Python標題Re: [問題] 請問如何產生N個名稱時間...
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
: 再次感謝。
--
※ 發信站: 批踢踢實業坊(ptt.cc), 來自: 130.203.95.180
※ 文章網址: https://www.ptt.cc/bbs/Python/M.1453395445.A.6EB.html
推 alen332l: 推 01/22 09:25
推 hsnusonic: 推 01/22 09:54
推 lefan: 感謝。 01/22 10:42