N体运动的程序模拟

时间：2014-09-24 14:13:47 阅读：413 评论：0 收藏：0 [点我收藏+]

标签：des style blog http color io os 使用 ar

这依然是与《三体》有关的一篇文章。空间中三个星体在万有引力作用下的运动被称之为三体问题，参见我的上一篇文章：三体运动的程序模拟。而这一节，对三体问题进行了扩展，实现了空间中N个星体在万有引力下的运动模拟。

程序中使用了两个物理定律：

（1）万有引力定律

这是牛顿发现的:任意两个质点有通过连心线方向上的力相互吸引。该引力的大小与它们的质量乘积成正比，与它们距离的平方成反比，与两物体的化学本质或物理状态以及中介物质无关。

以数学表示为： $bubuko.com,布布扣$

其中：

$bubuko.com,布布扣$ : 两个物体之间的引力
$bubuko.com,布布扣$ : 万有引力常数
$bubuko.com,布布扣$ : 物体1的质量
$bubuko.com,布布扣$ : 物体2的质量
$bubuko.com,布布扣$ : 两个物体之间的距离

依照国际单位制，F的单位为牛顿(N)，m₁和m₂的单位为千克(kg)，r 的单位为米(m)，常数G近似地等于6.67 × 10⁻¹¹ N m² kg⁻²（牛顿米的平方每千克的平方）。当然程序中G不可能设置为这么小的一个数,用户可以自己设置万有引力常数,以实现合理的运动.

(2)势能与动能守恒定律

引力位能或引力势能是指物体因为大质量物体的万有引力而具有的位能，其大小与其到大质量的距离有关。

$bubuko.com,布布扣$

其中G为万有引力常数，M、m分别为两物体质量，r为两者距离。

动能，简单的说就是指物体因运动而具有的能量。数值上等于（1/2）Mv^2. 动能是能量的一种，它的国际单位制下单位是焦耳(J)，简称焦。

势能与动能守恒即在整个运动过程中,其总能量是不变的,即N体的动能与势能之和为固定值.

本来我只打算使用万有引力定律,没有考虑势能动能守恒定律.但程序运行后发现,由于没有使用微积分,即使采样间隔时间再小,误差也很大.这让我想起当年学大学物理的时候,大一上学期开高数,然后下学期开了物理.现在对大学物理讲的东西,只记得是用高数中的微积分解物理问题.可那时,我高数学得就是个渣,所以物理也没办法学好.当然就算物理学好了,现在也早忘光了.没有用微积分,这个程序中算法误差不小,模拟一下吧,不要太计较.

核心代码也发一下,希望有兴趣的人可以对其做出优化:

  1 /****************************************************************
  2 
  3   File name   :  NBodyKernel.h
  4   Author      :  叶峰
  5   Version     :  1.0a
  6   Create Date :  2014/09/19   
  7   Description :  
  8 
  9 *****************************************************************/
 10 
 11 // --------------------------------------------------------------------------------------
 12 
 13 #ifndef _NBodyKernel_H_
 14 #define _NBodyKernel_H_
 15 
 16 // INCLUDES -----------------------------------------------------------------------------
 17 
 18 #include "..\..\..\Why\YCommon_h\WhyMath.h"
 19 #include "..\..\..\Why\YInterface\Others\YdD3D9Math.h"
 20 
 21 // --------------------------------------------------------------------------------------
 22 
 23 #define YD_ORBIT_VERTICES_COUNT         1024
 24 
 25 // --------------------------------------------------------------------------------------
 26 
 27 class Star
 28 {
 29 public:
 30     Vector3 vPosition;          // 位置
 31     Vector3 vVelocity;          // 速度
 32     Vector3 vForce;             // 受力
 33     Yreal   fWeight;            // 质量
 34     Yreal   fRadiusScale;       // 缩放
 35     Yreal   fKineticEnergy;     // 动能
 36     Ycolor  color;              // 颜色
 37     Ybool   bOrbitVisible;      // 轨迹可见
 38     Yuint   id;
 39     Vector3 verticesOrbit[YD_ORBIT_VERTICES_COUNT];
 40     Star*   pNext;
 41 
 42     Star()
 43     {
 44         fWeight = 0.0f;
 45         fRadiusScale = 1.0f;  
 46         fKineticEnergy = 0.0f;
 47         color = 0xffffffff;
 48         bOrbitVisible = YD_FALSE;
 49         id = -1;
 50         memset(verticesOrbit, 0, sizeof(verticesOrbit));
 51         pNext = NULL;
 52     }
 53 
 54     // 根据速度计算动能
 55     void CalculateKineticEnergyByVelocity()
 56     {
 57         Yreal sqv = D3DXVec3LengthSq(&vVelocity);
 58         fKineticEnergy = 0.5f*fWeight*sqv;
 59     }
 60 
 61     // 根据动能计算速度
 62     void CalculateVelocityByKineticEnergy()
 63     {
 64         float v = sqrt(2*fKineticEnergy/fWeight);
 65         float w = D3DXVec3Length(&vVelocity);
 66         vVelocity *= v/w;
 67     };
 68 };
 69 
 70 // --------------------------------------------------------------------------------------
 71 
 72 class CNBodyKernel
 73 {
 74 public:
 75     CNBodyKernel();
 76 
 77     ~CNBodyKernel();
 78 
 79     // 更新N体
 80     void        UpdateNBody(Yreal deltaTime);
 81 
 82     // 清空N体
 83     void        Clear();
 84 
 85     // 设置引力系数
 86     void        SetGravitationCoefficient(Yreal c);
 87 
 88     // 获得引力系数
 89     Yreal       GetGravitationCoefficient() const
 90     {
 91         return m_fGravitationCoefficient;
 92     }
 93 
 94     // 返回星体数目
 95     Yuint       GetStarsCount() const
 96     {
 97         return m_starsCount;
 98     }
 99 
100     // 返回星体链表
101     const Star* GetStarsList() const
102     {
103         return m_starsList;
104     }
105 
106     // 返回星体对象
107     const Star* GetStar(Yuint index) const;
108 
109     // 移除星体
110     bool        RemoveStar(Yuint index);
111 
112     // 添加星体
113     void        AddStar(
114         const Vector3& vPosition,
115         const Vector3& vVelocity,
116         const Yreal fWeight,
117         Ycolor  color);
118 
119 
120     // 设置星体的重量
121     void        SetStarWeight(Yuint index, Yreal weight);
122 
123     // 设置星体的轨迹是否可见
124     void        SetStarOrbitVisible(Yuint index, bool visible);
125 
126     // 设置星体的颜色
127     void        SetStarColor(Yuint index, Ycolor color);
128 
129     Yuint       GetCurrentOrbit() const
130     {
131         return m_currentOrbit;
132     }
133 
134     const Vector3* GetStarPosition(Yuint index) const;
135 
136     Yreal       GetStarWeight(Yuint index) const;
137 
138     bool        GetStarOrbitVisible(Yuint index) const;
139 
140     Ycolor      GetStarColor(Yuint index) const;
141 
142     Yuint       GetStarID(Yuint index) const;
143 
144     void        GetSpaceCenter(Vector3& vCenter) const;
145 
146     void        GetWeightCenter(Vector3& vCenter) const; 
147 
148 private:
149     // 更新星体的位置,速度,动能
150     void        UpdateStar(Star& star, Yreal t);
151 
152     // 根据动能重新计算速度
153     void        RecalculateStarVelocity(Star& star);
154 
155     // 计算每个星体的当前受力
156     void        CalculateStarsForce();
157 
158     // 计算总势能
159     void        CalculateAmountOfPotentialEnergy();
160 
161     // 计算总动能
162     void        CalculateAmountOfKineticEnergy();
163 
164     // 计算总能量
165     void        CalculateAmountEnergy()
166     {
167         CalculateAmountOfPotentialEnergy();
168         CalculateAmountOfKineticEnergy();
169         m_fAmountEnergy = m_fAmountOfPotentialEnergy + m_fAmountOfKineticEnergy;
170     }
171 
172 private:
173     Yreal m_fGravitationCoefficient;    // 引力系数
174     Yreal m_fAmountOfPotentialEnergy;   // 当前总势能
175     Yreal m_fAmountOfKineticEnergy;     // 当前总动能
176     Yreal m_fAmountEnergy;              // 当前能量m_fAmountOfPotentialEnergy + m_fAmountOfKineticEnergy;
177 
178     Star* m_starsList;                  // N体链表
179     Yuint m_starsCount;                 // N体数目
180 
181     Yuint m_currentOrbit;
182 };
183 
184 // --------------------------------------------------------------------------------------
185 
186 #endif

  1 /****************************************************************
  2 
  3   File name   :  NBodyKernel.cpp
  4   Author      :  叶峰
  5   Version     :  1.0a
  6   Create Date :  2014/09/19
  7   Description :  
  8 
  9 *****************************************************************/
 10 
 11 // --------------------------------------------------------------------------------------
 12 
 13 #include "..\..\..\Why\YCommon_h\WhyMemoryDefines.h"
 14 #include "NBodyKernel.h"
 15 #include <assert.h>
 16 
 17 // --------------------------------------------------------------------------------------
 18 
 19 #define YD_ESP                          0.00001f
 20 
 21 // --------------------------------------------------------------------------------------
 22 
 23 CNBodyKernel::CNBodyKernel()
 24 {
 25     m_fGravitationCoefficient = 100.0f;
 26     m_fAmountOfPotentialEnergy = 0.0f;
 27     m_fAmountOfKineticEnergy = 0.0f;
 28     m_fAmountEnergy = 0.0f;
 29 
 30     m_starsList = NULL;
 31     m_starsCount = 0;
 32 
 33     m_currentOrbit = 0;
 34 } 
 35 
 36 CNBodyKernel::~CNBodyKernel()
 37 {
 38     Clear();
 39 }
 40 
 41 // 清空N体
 42 void        CNBodyKernel::Clear()
 43 {
 44     Star* pDel = m_starsList;
 45     Star* pNext;
 46     while (pDel)
 47     {
 48         pNext = pDel->pNext;
 49         YD_DELETE(pDel);
 50         pDel = pNext;
 51     }
 52 
 53     m_starsList = NULL;
 54     m_starsCount = 0;
 55 }
 56 
 57 // 获得引力系数
 58 void        CNBodyKernel::SetGravitationCoefficient(Yreal c)
 59 {
 60     m_fGravitationCoefficient = c;
 61 
 62     // 计算总能量
 63     CalculateAmountEnergy();
 64 }
 65 
 66 // 计算每个星体的当前受力
 67 void        CNBodyKernel::CalculateStarsForce()
 68 {
 69     // 清空所有引力
 70     Star* pStar = m_starsList;
 71     while (pStar)
 72     {
 73         pStar->vForce = Vector3(0.0f, 0.0f, 0.0f);
 74         pStar = pStar->pNext;
 75     }
 76 
 77     // 计算引力
 78     Star* pCurrent = m_starsList;
 79     Star* pNext;
 80     while (pCurrent)
 81     {
 82         pNext = pCurrent->pNext;
 83         while (pNext)
 84         {
 85             Vector3 vAB = pNext->vPosition - pCurrent->vPosition;
 86             float disSqAB = D3DXVec3LengthSq(&vAB);
 87             if (disSqAB < YD_ESP)
 88             {
 89                 disSqAB = YD_ESP;
 90             }
 91             float disAB = sqrt(disSqAB);
 92             Vector3 vForceDir = vAB/disAB;
 93             Yreal fForce = m_fGravitationCoefficient*pCurrent->fWeight*pNext->fWeight/disSqAB;
 94             Vector3 vForce = vForceDir*fForce;
 95 
 96             pCurrent->vForce += vForce;
 97             pNext->vForce -= vForce;
 98 
 99             pNext = pNext->pNext;
100         }
101         pCurrent = pCurrent->pNext;
102     }
103 }
104 
105 void        CNBodyKernel::UpdateStar(Star& star, Yreal t)
106 {
107     // 加速度
108     const Vector3 acc = star.vForce/star.fWeight;
109 
110     star.vPosition.x = star.vPosition.x + star.vVelocity.x*t + 0.5f*acc.x*t*t;
111     star.vPosition.y = star.vPosition.y + star.vVelocity.y*t + 0.5f*acc.y*t*t;
112     star.vPosition.z = star.vPosition.z + star.vVelocity.z*t + 0.5f*acc.z*t*t;
113 
114     star.vVelocity.x += acc.x*t;
115     star.vVelocity.y += acc.y*t;
116     star.vVelocity.z += acc.z*t;
117 
118     // 重新计算动能
119     star.CalculateKineticEnergyByVelocity();
120 }
121 
122 void        CNBodyKernel::UpdateNBody(Yreal deltaTime)
123 {
124     // 计算每个星体的当前受力
125     CalculateStarsForce();
126 
127     // 更新星体的位置与速度
128     Star* pStar = m_starsList;
129     while (pStar)
130     {
131         UpdateStar(*pStar, deltaTime);
132         pStar = pStar->pNext;
133     }
134 
135     // 能量守恒
136     CalculateAmountOfKineticEnergy();
137     CalculateAmountOfPotentialEnergy();
138 
139     // 理论上的动能
140     Yreal fKineticEnergy = m_fAmountEnergy - m_fAmountOfPotentialEnergy;
141     if (fKineticEnergy < 0.0f)
142     {
143         fKineticEnergy = 0.0f;
144     }
145 
146     Yreal r = fKineticEnergy/m_fAmountOfKineticEnergy;
147     pStar = m_starsList;
148     while (pStar)
149     {
150         if (r > YD_ESP)
151         {
152             pStar->fKineticEnergy *= r;
153         }
154         else
155         {
156             pStar->fKineticEnergy = 1.0f;
157         }
158         // 根据动能计算速度
159         pStar->CalculateVelocityByKineticEnergy();
160         pStar = pStar->pNext;
161     }
162 
163     // 更新轨迹点
164     pStar = m_starsList;
165     while (pStar)
166     {
167         pStar->verticesOrbit[m_currentOrbit] = pStar->vPosition;
168         pStar = pStar->pNext;
169     }
170     m_currentOrbit++;
171     if (m_currentOrbit >= YD_ORBIT_VERTICES_COUNT)
172     {
173         m_currentOrbit = 0;
174     }
175 }
176 
177 // 返回星体对象
178 const Star* CNBodyKernel::GetStar(Yuint index) const
179 {
180     if (index >= m_starsCount)
181     {
182         return NULL;
183     }
184 
185     const Star* pStar = m_starsList;
186     while (index && pStar)
187     {
188         index--;
189         pStar = pStar->pNext;
190     }
191 
192     return pStar;
193 }
194 
195 // 移除星体
196 bool        CNBodyKernel::RemoveStar(Yuint index)
197 {
198     if (index >= m_starsCount)
199     {
200         return false;
201     }
202 
203     if (index == 0)
204     {
205         Star* pStar = m_starsList->pNext;
206         YD_DELETE(m_starsList);
207         m_starsList = pStar;
208     }
209     else
210     {
211         Star* pStar = m_starsList;
212         Yuint t = index - 1;
213         while (t && pStar)
214         {
215             t--;
216             pStar = pStar->pNext;
217         }
218         assert(pStar);
219 
220         Star* pDel = pStar->pNext;
221         pStar->pNext = pDel->pNext;
222         YD_DELETE(pDel);
223     }
224 
225     m_starsCount--;
226 
227     Yuint id = 0;
228     Star* pStar = m_starsList;
229     while (pStar)
230     {
231         pStar->id = id;
232         id++;
233         pStar = pStar->pNext;
234     }
235 
236     // 重算能量
237     CalculateAmountEnergy();
238 
239     return true;
240 }
241 
242 // 添加星体
243 void        CNBodyKernel::AddStar(
244                     const Vector3& vPosition,
245                     const Vector3& vVelocity,
246                     const Yreal fWeight,
247                     Ycolor  color)
248 {
249     Star* pNewStar = YD_NEW(Star);
250     pNewStar->vPosition = vPosition;
251     pNewStar->vVelocity = vVelocity;
252     pNewStar->fWeight = fWeight;
253     pNewStar->fRadiusScale = yf_pow(fWeight, 1.0f/3.0f);
254     pNewStar->CalculateKineticEnergyByVelocity();
255     pNewStar->color = color;
256     pNewStar->id = m_starsCount;
257     for (Yuint i = 0; i < YD_ORBIT_VERTICES_COUNT; i++)
258     {
259         pNewStar->verticesOrbit[i] = vPosition;
260     }
261 
262     if (!m_starsList)
263     {
264         m_starsList = pNewStar;
265     }
266     else
267     {
268         Star* pStar = m_starsList;
269         while (pStar->pNext)
270         {
271             pStar = pStar->pNext;
272         }
273         pStar->pNext = pNewStar;
274     }
275 
276     m_starsCount++;
277 
278     // 重算能量
279     CalculateAmountEnergy();
280 }
281 
282 // 设置星体的重量
283 void        CNBodyKernel::SetStarWeight(Yuint index, Yreal weight)
284 {
285     if (index >= m_starsCount)
286     {
287         return;
288     }
289 
290     Star* pStar = m_starsList;
291     while (index && pStar)
292     {
293         index--;
294         pStar = pStar->pNext;
295     }
296 
297     pStar->fWeight = weight;
298     pStar->fRadiusScale = yf_pow(weight, 1.0f/3.0f);
299     pStar->CalculateKineticEnergyByVelocity();
300 
301     // 重算能量
302     CalculateAmountEnergy();
303 }
304 
305 // 设置星体的轨迹是否可见
306 void        CNBodyKernel::SetStarOrbitVisible(Yuint index, bool visible)
307 {
308     if (index >= m_starsCount)
309     {
310         return;
311     }
312 
313     Star* pStar = m_starsList;
314     while (index && pStar)
315     {
316         index--;
317         pStar = pStar->pNext;
318     }
319 
320     pStar->bOrbitVisible = visible;
321 }
322 
323 // 设置星体的颜色
324 void        CNBodyKernel::SetStarColor(Yuint index, Ycolor color)
325 {
326     if (index >= m_starsCount)
327     {
328         return;
329     }
330 
331     Star* pStar = m_starsList;
332     while (index && pStar)
333     {
334         index--;
335         pStar = pStar->pNext;
336     }
337 
338     pStar->color = color;
339 }
340 
341 // 计算总动能
342 void        CNBodyKernel::CalculateAmountOfKineticEnergy()
343 {
344     // 动能
345     m_fAmountOfKineticEnergy = 0.0f;
346     Star* pStar = m_starsList;
347     while (pStar)
348     {
349         //pStar->CalculateKineticEnergyByVelocity();
350         m_fAmountOfKineticEnergy += pStar->fKineticEnergy;
351         pStar = pStar->pNext;
352     }
353 }
354 
355 // 计算总势能
356 void       CNBodyKernel::CalculateAmountOfPotentialEnergy()
357 {
358     m_fAmountOfPotentialEnergy = 0.0f;
359 
360     Star* pCurrent = m_starsList;
361     Star* pNext;
362     while (pCurrent)
363     {
364         pNext = pCurrent->pNext;
365         while (pNext)
366         {
367             Vector3 vAB = pCurrent->vPosition - pNext->vPosition;
368             float disAB = D3DXVec3Length(&vAB);
369             float fPotentialEnergyAB = -m_fGravitationCoefficient*pCurrent->fWeight*pNext->fWeight/disAB;
370             m_fAmountOfPotentialEnergy += fPotentialEnergyAB;
371             pNext = pNext->pNext;
372         }
373         pCurrent = pCurrent->pNext;
374     }
375 }
376 
377 const Vector3* CNBodyKernel::GetStarPosition(Yuint index) const
378 {
379     if (index >= m_starsCount)
380     {
381         return 0;
382     }
383 
384     const Star* pStar = m_starsList;
385     while (index && pStar)
386     {
387         index--;
388         pStar = pStar->pNext;
389     }
390 
391     assert(pStar);
392     return &pStar->vPosition;
393 }
394 
395 Yreal       CNBodyKernel::GetStarWeight(Yuint index) const
396 {
397     if (index >= m_starsCount)
398     {
399         return 0.0f;
400     }
401 
402     const Star* pStar = m_starsList;
403     while (index && pStar)
404     {
405         index--;
406         pStar = pStar->pNext;
407     }
408 
409     assert(pStar);
410     return pStar->fWeight;
411 }
412 
413 bool        CNBodyKernel::GetStarOrbitVisible(Yuint index) const
414 {
415     if (index >= m_starsCount)
416     {
417         return false;
418     }
419 
420     const Star* pStar = m_starsList;
421     while (index && pStar)
422     {
423         index--;
424         pStar = pStar->pNext;
425     }
426 
427     assert(pStar);
428     return pStar->bOrbitVisible != YD_FALSE;
429 }
430 
431 Ycolor      CNBodyKernel::GetStarColor(Yuint index) const
432 {
433     if (index >= m_starsCount)
434     {
435         return 0;
436     }
437 
438     const Star* pStar = m_starsList;
439     while (index && pStar)
440     {
441         index--;
442         pStar = pStar->pNext;
443     }
444 
445     assert(pStar);
446     return pStar->color;
447 }
448 
449 Yuint      CNBodyKernel::GetStarID(Yuint index) const
450 {
451     if (index >= m_starsCount)
452     {
453         return 0;
454     }
455 
456     const Star* pStar = m_starsList;
457     while (index && pStar)
458     {
459         index--;
460         pStar = pStar->pNext;
461     }
462 
463     assert(pStar);
464     return pStar->id;
465 }
466 
467 void        CNBodyKernel::GetSpaceCenter(Vector3& vCenter) const
468 {
469     vCenter = Vector3(0.0f, 0.0f, 0.0f);
470     if (!m_starsCount)
471     {
472         return;
473     }
474 
475     Star* pStar = m_starsList;
476     while (pStar)
477     {
478         vCenter += pStar->vPosition;
479         pStar = pStar->pNext;
480     }
481 
482     vCenter /= (Yreal)m_starsCount;
483 }
484 
485 void        CNBodyKernel::GetWeightCenter(Vector3& vCenter) const
486 {
487     vCenter = Vector3(0.0f, 0.0f, 0.0f);
488     if (!m_starsCount)
489     {
490         return;
491     }
492 
493     Yreal weights = 0.0f;
494     Star* pStar = m_starsList;
495     while (pStar)
496     {
497         vCenter += pStar->vPosition*pStar->fWeight;
498         weights += pStar->fWeight;
499         pStar = pStar->pNext;
500     }
501 
502     vCenter /= weights;
503 }