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://webptt.com/m.aspx?n=bbs/Python/M.1453395445.A.6EB.html