ゼロ重力2次元空間における粒子の重力計算の最適化

Question

私は、粒子の小さな視覚化をパイソンで作成しました。 私はゼロ重力の2D空間で粒子の動きを解明しています。 各粒子は、粒子の質量および距離に基づいて、他のすべての粒子を引き付ける。

私はpygameでビジュアライゼーションを行いましたが、すべてがプランとして機能しますが（caluclationあり）、extreamly計算を最適化する必要があります。今日、システムは、横方向のフレームレートにおいて約100〜150個の粒子を計算することができる。すべての計算を別のスレッドに入れてくれました。それは私に何かを与えてくれましたが、私が望むものはほとんどありませんでした。

私はscipyとnumpyを見ましたが、私は科学者でもmathguruでもないので、ちょっと混乱します。私は正しい軌道にいるように見えますが、私は手がかりを持っていません。

私はループ内のすべてのパーティクル上のすべてのアトラクションを計算する必要があります。 そして、私が衝突したかどうかを知る必要があるので、私は同じことをやり直す必要があります。

それはコードのようなものを書くために私の心を壊し....

numpyのは、配列の配列を計算する能力を持っている、しかし、私は見つかっていない全てのアイテムを配列内のすべての項目を計算するためにどのような任意の同じ/別の配列から。 1つはありますか？ ので、私は作成して、配列のカップルとはるかに高速に計算し、その値が一致（Collitiondetect IOW）ここで

は、今日の魅力/ collsion計算で2つのアレイへのインデックスを取得する機能が存在しなければならない可能性がある場合：

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 

を持っていますが、

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

を行い、その後、そのように位置を更新してしまうdelta_Tの時刻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)に劣化します。

もう1つの方法は、粒子をOctreeの内側に置くことです。 1つを構築することはO(n)のようなものです。それからそれを照会して、お互いが近くにあるかどうかを調べることができます。その時点でペア上で衝突検出を行うだけです。これは私が信じているO(n log n)です。

それだけでなく、Octreeを使って重力計算を高速化することもできます。 O(n^2)の動作の代わりに、O(n log n)にも低下します。ほとんどのOctreeの実装には、速度と正確なトレードオフを制御する「開始パラメータ」が含まれています。だからOctreeは直接ペアワイズ計算よりも精度が低く、コード化が複雑ですが、大規模なシミュレーションも可能です。

Octreeをこのように使用すると、Barnes-Hut Simulationと呼ばれることが起こります。

注：2Dで作業しているため、Octreeの2Dアナログは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は、粒子のx、y位置を記憶する5x2アレイである。各ペア間の距離を計算するには、使用することができます：

[[ 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.  ]] 

か、scipyのダウンロードからcdistを使用することができます：

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

を出力します

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 '値の行列は、1/abs(z-z')**3を計算するときには、Z-Z[:, numpy.newaxis]の座標の1D配列Zから得られます（対角項[z = z']は特別な注意が必要です。

時間発展に関しては、確かにSciPyの高速微分方程式ルーチンを使用することができます。これは、ステップオイラー積分よりもはるかに正確です。

いずれにしても、特にNumPyが非常に高速であるため、科学計算を計画している場合は、NumPyを掘り下げることは非常に便利です。

Answer 5

これがあなたを助けてくれるのかどうかはわかりませんが、同じ問題で私が取り組んでいる解決策の一部です。このようにしてパフォーマンスを大幅に向上させることはできませんでしたが、まだ200個ほどの粒子が停止するようになり始めましたが、多分それはあなたにいくつかのアイデアを与えるでしょう。2D平面上に引力のxとy成分を計算するための

C++モジュール：あなたはPYDファイルとしてこれをコンパイルする場合

#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); 
} 

はあなたがインポートすることができPythonオブジェクトを取得する必要があります。あなたはすべてのあなたのコンパイラとリンカーのオプションを正しく取得する必要があります。私はdev-C++を使用しています。コンパイラオプションを-shared -o gforces.pydに設定し、リンカを-lpython27に設定します（インストールした同じバージョンを使用してください）。そして、pythonディレクトリパスをインクルードとライブラリに追加しますタブ。

オブジェクトは引数（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 

・ホープこのことができます少し。それ以上の助言や私の部分の全体的なエラーを指摘することは歓迎します、私はこれにも新しいです。