在零重力二維空間中優化粒子的引力計算

Question

我已經在python中創建了一個粒子的小型可視化。 我是在零重力的二維空間中粒子的運動。由於每個粒子都基於粒子質量和距離吸引所有其他粒子。

我在pygame中做了一個可視化，並且所有工作都按照計劃進行（有了計算），但是我需要優化計算extreamly。今天，該系統可以計算大約100-150個顆粒的相鄰幀率。我把所有的計算都放在一個單獨的線程中，給了我更多但並不是我想要的東西。

我看着scipy和numpy，但因爲我不是科學家或mathguru我只是感到困惑。它看起來像我在正確的軌道上，但我不知道如何。

我需要計算所有粒子上的所有吸引力，我必須在一個循環中循環。 因爲我需要找到是否有任何碰撞，我必須再一次做同樣的事情。

它打破了我的心臟寫那種代碼....

numpy的具有計算與陣列陣列的能力，但是我還沒有發現任何東西來計算所有的項目在數組中的所有項目來自同一個/另一個陣列。有一個嗎？ 如果是這樣我可以創建和對數組和計算速度更快，而且必須是從2個陣列，其中它們的值匹配（Collitiondetect IOW）

下面是今天的吸引/ collsion計算得到指數函數：

class Particle: 
    def __init__(self): 
     self.x = random.randint(10,790) 
     self.y = random.randint(10,590) 
     self.speedx = 0.0 
     self.speedy = 0.0 
     self.mass = 4 

#Attraction  
for p in Particles: 
    for p2 in Particles: 
     if p != p2: 
      xdiff = P.x - P2.x 
      ydiff = P.y - P2.y 
      dist = math.sqrt((xdiff**2)+(ydiff**2)) 
      force = 0.125*(p.mass*p2.mass)/(dist**2) 
      acceleration = force/p.mass 
      xc = xdiff/dist 
      yc = ydiff/dist 
      P.speedx -= acceleration * xc 
      P.speedy -= acceleration * yc 
for p in Particles: 
    p.x += p.speedx 
    p.y += p.speedy 

#Collision 
for P in Particles: 
    for P2 in Particles: 
     if p != P2: 
      Distance = math.sqrt( ((p.x-P2.x)**2) + ((p.y-P2.y)**2) ) 
      if Distance < (p.radius+P2.radius): 
       p.speedx = ((p.mass*p.speedx)+(P2.mass*P2.speedx))/(p.mass+P2.mass) 
       p.speedy = ((p.mass*p.speedy)+(P2.mass*P2.speedy))/(p.mass+P2.mass) 
       p.x = ((p.mass*p.x)+(P2.mass*P2.x))/(p.mass+P2.mass) 
       p.y = ((p.mass*p.y)+(P2.mass*P2.y))/(p.mass+P2.mass) 
       p.mass += P2.mass 
       p.radius = math.sqrt(p.mass) 
       Particles.remove(P2) 

Answer 1

（這可能應該去評論，但我沒有必要的聲譽這樣做）

我不看你怎麼辦步進的時間。你有

P.speedx -= acceleration * xc 
P.speedy -= acceleration * yc 

但要獲得新的速度在時間t + delta_t你會做

P.speedx -= acceleration * xc * delta_t 
P.speedy -= acceleration * yc * delta_t 

，然後更新的位置，像這樣：

P.x = P.x + P.speedx * delta_t 
P.y = P.y + P.speedy * delta_t 

然後到你的速度關注。也許最好將粒子信息存儲在不是在類中，而是在numpy數組中存儲？但我不認爲你可以避免循環。

另外，你看看wikipedia，那裏描述了一些加快計算的方法。

（編輯由於麥克的評論）

Answer 2

我已經在這之前曾和我見過的過去，加速碰撞計算的事情之一是實際存儲附近粒子的列表。

基本上，這個想法是你的重心計算裏面你做這樣的事情：

for (int i = 0; i < n; i++) 
{ 
    for (int j = i + 1; j < n; j++) 
    { 
     DoGravity(Particle[i], Particle[j]); 
     if (IsClose(Particle[i], Particle[j])) 
     { 
      Particle[i].AddNeighbor(Particle[j]); 
      Particle[j].AddNeighbor(Particle[i]); 
     } 
    } 
} 

然後，你只需通過對所有粒子和你做的每個碰撞檢測上又將。在最好的情況下，這通常類似於O(n)，但在最壞的情況下，它很容易降級到O(n^2)。

另一種替代方法是嘗試將您的粒子放在Octree的內部。構建一個類似於O(n)的東西，然後您可以查詢它以查看是否有任何東西彼此靠近。在這一點上，你只需對對進行碰撞檢測。這樣做是O(n log n)我相信。

不僅如此，還可以使用八叉樹來加速重力計算。而不是O(n^2)行爲，它也下降到O(n log n)。大多數八叉樹實現包含一個「開放參數」，用於控制您將要進行的速度與準確性的平衡。所以Octrees往往不如直接成對計算準確，編碼複雜，但它們也可以進行大規模的模擬。

如果您以這種方式使用八叉樹，您將執行所謂的Barnes-Hut Simulation。

注意：由於您使用的是2D工作模式，所以2D模擬工具Octree被稱爲Quadtree。有關詳細信息，請參見以下Wikipedia文章：http://en.wikipedia.org/wiki/Quadtree

Answer 3

要做快速計算，您需要在numpy數組中存儲x，y，speedx，speedy，m。例如：

import numpy as np 

p = np.array([ 
    (0,0), 
    (1,0), 
    (0,1), 
    (1,1), 
    (2,2), 
], dtype = np.float) 

p是一個5x2陣列，它存儲粒子的x，y位置。要計算每對之間的距離，你可以使用：

print np.sqrt(np.sum((p[:, np.newaxis] - p[np.newaxis, :])**2, axis=-1)) 

輸出爲：

[[ 0.   1.   1.   1.41421356 2.82842712] 
[ 1.   0.   1.41421356 1.   2.23606798] 
[ 1.   1.41421356 0.   1.   2.23606798] 
[ 1.41421356 1.   1.   0.   1.41421356] 
[ 2.82842712 2.23606798 2.23606798 1.41421356 0.  ]] 

，或者您可以使用cdist從SciPy的：

from scipy.spatial.distance import cdist 
print cdist(p, p) 

Answer 4

你可以先嚐試與複數：相關的引力和動力學公式在這種形式中非常簡單，並且也可以很快（因爲NumPy可以進行計算在內部，而不是分別處理x和y座標）。例如，在Z和Z 2個particules之間的力「很簡單：

(z-z')/abs(z-z')**3 

NumPy的可以計算這樣一個數量非常迅速，對於所有Z/Z」對。例如，所有z-z'值的矩陣僅從座標的一維陣列Z獲得，如Z-Z[:, numpy.newaxis]（當計算1/abs(z-z')**3時，對角線項[z = z']確實需要特別小心：它們應被設置爲零）。

至於時間的演變，你當然可以使用SciPy的快速微分方程程序：它們比逐步歐拉積分更精確。

無論如何，鑽研NumPy將會非常有用，特別是如果您計劃進行科學計算，因爲NumPy速度非常快。

Answer 5

不知道這是否會對您有所幫助，但它是我一直致力於解決同一問題的解決方案的一部分。我並沒有注意到在這樣的表現中獲得了巨大的收益，仍然開始遏制200個粒子，但也許它會給你一些想法。

C++計算2D平面上的引力的x和y分量模塊：

#include <Python.h> 
#include <math.h> 

double _acceleration(double &Vxa, double &Vya, double &Vxb, double &Vyb, double xa, double ya, double xb, double yb, double massa, double massb) 
{ 
    double xdiff = xa - xb; 
    double ydiff = ya - yb; 
    double distance = sqrt(xdiff*xdiff + ydiff*ydiff) * pow(10, 5); 

    if (distance <= 0) 
     distance = 1; 

    double force = (6.674 * pow(10, -11))*(massa*massb)/(distance*distance); 

    double acca = force/massa; 
    double accb = force/massb; 
    double xcomponent = xdiff/distance; 
    double ycomponent = ydiff/distance; 

    Vxa -= acca * xcomponent; 
    Vya -= acca * ycomponent; 
    Vxb += accb * xcomponent; 
    Vyb += accb * ycomponent; 

    return distance; 
} 

static PyObject* gforces(PyObject* self, PyObject* args) 
{ 
    double Vxa, Vya, Vxb, Vyb, xa, ya, xb, yb, massa, massb, distance; 

    if (!PyArg_ParseTuple(args, "dddddddddd", &Vxa, &Vya, &Vxb, &Vyb, &xa, &ya, &xb, &yb, &massa, &massb)) 
     return NULL; 

    distance = _acceleration(Vxa, Vya, Vxb, Vyb, xa, ya, xb, yb, massa, massb); 

    return Py_BuildValue("ddddd", Vxa, Vya, Vxb, Vyb, distance); 
} 

static PyMethodDef GForcesMethods[] = { 
{"gforces", gforces, METH_VARARGS, "Calculate the x and y acceleration of two masses and the distance between the two."}, 
{NULL, NULL, 0, NULL} 
}; 

PyMODINIT_FUNC 
initgforces(void) 
{ 
(void) Py_InitModule("gforces", GForcesMethods); 
} 

如果編譯這是一個PYD文件，你應該得到一個Python對象，你可以導入。你必須讓所有的編譯器和鏈接器選項正確。我正在使用dev-C++並將我的編譯器選項設置爲-shared -o gforces.pyd並將鏈接器設置爲-lpython27（確保您使用的是您安裝的相同版本）並將Python目錄路徑添加到includes和庫標籤。

該對象需要參數（p1.speedx，p1.speedy，p2.speedx，p2.speedy，p1.x，p1.y，p2.x，p2.y，p1.mass，p2.mass）並返回新的p1.speedx，p1.speedy，p2.speedx，p2.speedy和p1 p2之間的距離。

使用上面的模塊，我也試着通過返回的距離，這樣比較到粒子的半徑之和切出的碰撞檢測的幾個步驟：

def updateForces(self):   #part of a handler class for the particles 
    prevDone = [] 
    for i in self.planetGroup: #planetGroup is a sprite group from pygame 
     prevDone.append(i.ID) 
     for j in self.planetGroup: 
      if not (j.ID in prevDone):    #my particles have unique IDs 
       distance = i.calcGForce(j)   #calcGForce just calls the above 
       if distance <= i.radius + j.radius: #object and assigns the returned 
        #collision handling goes here #values for the x and y speed 
                #components to the particles 

希望這有助於一點點。任何進一步的建議或指出我的重大錯誤是值得歡迎的，我也是新手。