雙擺的程序實現

本文轉載自查看原文 2015-02-27 12:52 7037 物理

雙擺是物理中的一個概念，先給下單擺與雙擺的定義:

單擺：由一根不可伸長、質量不計的繩子，上端固定，下端系一個質點的裝置。

雙擺：是一個擺的支點裝在另一擺的下部所形成的組合物體。雙擺有兩個擺角，所以有兩個自由度。雙擺是多自由度振動系統的最簡單的力學模型之一，它也是一種混沌實例。

這里對雙擺的實現程序與上一篇文章中三體的實現很相似，要實現它的程序需要一定的數學基礎。

下面兩個GIF動畫圖像為雙擺的錄屏:

先寫了一個雙擺的抽象基類：

  1 // --------------------------------------------------------------------------------------
  2 
  3 inline float     rand2(float a, float b)
  4 {
  5     const float r = (float)(::rand() / ((float)RAND_MAX + 1));
  6     return r*(b-a) + a;
  7 }
  8 
  9 // --------------------------------------------------------------------------------------
 10 
 11 class IDoublePendulum
 12 {
 13 public:
 14 
 15     IDoublePendulum()
 16     {
 17         m_fGravity = 9.8f;
 18 
 19         m_m1 = 1.0f;
 20         m_m2 = 2.0f;
 21     
 22         m_l1 = 1.0f;
 23         m_l2 = 2.0f;
 24     }
 25 
 26     virtual ~IDoublePendulum() {}
 27 
 28     // 重置
 29     virtual void        Reset()
 30     {
 31         m_m1 = rand2(1.0f, 5.0f);
 32         m_m2 = rand2(1.0f, 5.0f);
 33 
 34         m_l1 = rand2(1.0f, 5.0f);
 35         m_l2 = rand2(1.0f, 5.0f);
 36     }
 37 
 38     // 按時更新
 39     virtual void        Update(float deltaTime) = NULL;
 40 
 41     // 重力系數
 42     virtual void        SetGravity(float g)
 43     {
 44         m_fGravity = g;
 45     }
 46 
 47     // 質量
 48     virtual void        SetMass1(float m)
 49     {
 50         m_m1 = m;
 51     }
 52 
 53     virtual void        SetMass2(float m)
 54     {
 55         m_m2 = m;
 56     }
 57 
 58     // 長度
 59     virtual void        SetLength1(float l)
 60     {
 61         m_l1 = l;
 62         UpdatePosition();
 63     }
 64 
 65     virtual void        SetLength2(float l)
 66     {
 67         m_l2 = l;
 68         UpdatePosition();
 69     }
 70 
 71     float               GetGravity() const
 72     {
 73         return m_fGravity;
 74     }
 75 
 76     float               GetMass1() const
 77     {
 78         return m_m1;
 79     }
 80 
 81     float               GetMass2() const
 82     {
 83         return m_m2;
 84     }
 85 
 86     float               GetLength1() const
 87     {
 88         return m_l1;
 89     }
 90 
 91     float               GetLength2() const
 92     {
 93         return m_l2;
 94     }
 95 
 96     const Vec3&         GetPosition1() const
 97     {
 98         return m_pos1;
 99     }
100 
101     const Vec3&         GetPosition2() const
102     {
103         return m_pos2;
104     }
105 
106 protected:
107     // 更新兩球位置
108     virtual void        UpdatePosition() = NULL;
109 
110 protected:
111     float m_fGravity;   // 重力系數
112 
113     float m_m1, m_m2;   // 兩球質量
114     float m_l1, m_l2;   // 兩球距離
115 
116     Vec3 m_pos1, m_pos2;// 兩球位置
117 };

從IDoublePendulum中可以看到雙擺的幾個輸入數據是：重力系數，兩擺的距離和兩球質量。而計算的數據是：每一個時刻兩個球的位置。

下一步是對該基類進行繼承實現。

 1 // --------------------------------------------------------------------------------------
 2 
 3 #ifndef _DoublePendulum01_H_
 4 #define _DoublePendulum01_H_
 5 
 6 // --------------------------------------------------------------------------------------
 7 
 8 #include "IDoublePendulum.h"
 9 
10 // --------------------------------------------------------------------------------------
11 
12 class DoublePendulum01 : public IDoublePendulum
13 {
14 public:
15     DoublePendulum01()
16     {
17         m_a1 = 1.0f;
18         m_a2 = 2.0f;    
19 
20         m_v1 = 0.0f;
21         m_v2 = 0.0f;
22 
23         m_w1 = m_a1;
24         m_w2 = m_a2;
25     }
26 
27     // 重置
28     void        Reset();
29 
30     // 按時更新
31     void        Update(float deltaTime);
32 
33 protected:
34     void        UpdatePosition()
35     {
36         m_pos1.x = m_l1*sinf(m_w1);
37         m_pos1.y = -m_l1*cosf(m_w1);
38         m_pos1.z = 0.0f;
39 
40         m_pos2.x = m_pos1.x + m_l2*sinf(m_w2);
41         m_pos2.y = m_pos1.y - m_l2*cosf(m_w2);
42         m_pos2.z = 0.0f;
43     }
44 
45 private:
46     float m_a1, m_a2;   // 兩球初始角度
47     float m_w1, m_w2;   // 兩球當前角度
48     float m_v1, m_v2;   // 兩球的角加速度
49 };
50 
51 // --------------------------------------------------------------------------------------
52 
53 #endif

CPP

 1 // --------------------------------------------------------------------------------------
 2 
 3 #include "DoublePendulum01.h"
 4 
 5 // --------------------------------------------------------------------------------------
 6 
 7 /*
 8 求解線型方程組
 9 a*x + b*y + c = 0
10 d*x + e*y + f = 0
11 */
12 inline bool SolveLinalg(float a, float b, float c, float d, float e, float f, float& x, float& y)
13 {
14     float t = b*d - a*e;
15     if (fabs(t) < FLT_EPSILON)
16     {
17         x = 0.0f;
18         y = 0.0f;
19 
20         return false;
21     }
22 
23     x = (c*e - b*f)/t;
24     y = (a*f - c*d)/t;
25 
26     return true;
27 }
28 
29 // --------------------------------------------------------------------------------------
30 
31 void        DoublePendulum01::Reset()
32 {
33     IDoublePendulum::Reset();
34 
35     m_a1 = rand2(-2.0f, 2.0f);
36     m_a2 = rand2(-2.0f, 2.0f); 
37 
38     m_v1 = 0.0f;
39     m_v2 = 0.0f;
40 
41     m_w1 = m_a1;
42     m_w2 = m_a2;
43 
44     UpdatePosition();
45 }
46 
47 void        DoublePendulum01::Update(float deltaTime)
48 {
49     float a = m_l1*m_l1*(m_m1+m_m2);
50     float b = m_l1*m_m2*m_l2*cos(m_w1-m_w2);
51     float c = m_l1*(m_m2*m_l2*sin(m_w1-m_w2)*m_v2*m_v2 + (m_m1+m_m2)*m_fGravity*sin(m_w1));
52 
53     float d = m_m2*m_l2*m_l1*cos(m_w1-m_w2);
54     float e = m_m2*m_l2*m_l2;
55     float f = m_m2*m_l2*(-m_l1*sin(m_w1-m_w2)*m_v1*m_v1 + m_fGravity*sin(m_w2));
56 
57     // 角加速度
58     float dv1;
59     float dv2;
60     SolveLinalg(a, b, c, d, e, f, dv1, dv2);
61 
62     // 角速度
63     m_v1 += dv1*deltaTime;
64     m_v2 += dv2*deltaTime;
65 
66     // 角度
67     m_w1 += m_v1*deltaTime;
68     m_w2 += m_v2*deltaTime;
69 
70     UpdatePosition();
71 }

這里使用的是角度變化實現雙擺，參考的資料是：

http://sebug.net/paper/books/scipydoc/double_pendulum.html

它的理論有些難，我看得也似懂非懂的。好在公式就在那里，只要會用，不求會懂，程序就能實現。

這個程序是在2D的,雙擺可以出現在三維空間中.我也試着寫了下,使用力學來進行模擬,效果不太精確:

 1 // --------------------------------------------------------------------------------------
 2 
 3 #ifndef _DoublePendulum02_H_
 4 #define _DoublePendulum02_H_
 5 
 6 // --------------------------------------------------------------------------------------
 7 
 8 #include "IDoublePendulum.h"
 9 
10 // --------------------------------------------------------------------------------------
11 
12 class DoublePendulum02 : public IDoublePendulum
13 {
14 public:
15     DoublePendulum02()
16     {
17         m_velocity1 = Vec3(0.0f, 0.0f, 0.0f);
18         m_velocity2 = Vec3(0.0f, 0.0f, 0.0f);
19     }
20 
21     // 重置
22     void        Reset();
23 
24     // 按時更新
25     void        Update(float deltaTime);
26 
27 protected:
28     void        UpdatePosition()
29     {
30         Vec3 v = m_pos2 - m_pos1;
31         v.Normalize();
32 
33         m_pos1.Normalize();
34         m_pos1 *= m_l1;
35 
36         m_pos2 = m_pos1 + v*m_l2;
37     }
38 
39 private:
40     Vec3 m_velocity1, m_velocity2;  // 兩球當前速度
41 };
42 
43 // --------------------------------------------------------------------------------------
44 
45 #endif

View Code

.cpp

  1 // --------------------------------------------------------------------------------------
  2 
  3 #include "DoublePendulum02.h"
  4 
  5 // --------------------------------------------------------------------------------------
  6 
  7 void        DoublePendulum02::Reset()
  8 {
  9     IDoublePendulum::Reset();
 10 
 11     m_pos1.x = rand2(-1.0f, 1.0f);
 12     m_pos1.y = rand2(-1.0f, 0.0f);
 13     m_pos1.z = rand2(-1.0f, 1.0f);
 14     m_pos2.x = rand2(-1.0f, 1.0f);
 15     m_pos2.y = rand2(-0.5f, 0.5f);
 16     m_pos2.z = rand2(-1.0f, 1.0f);
 17     m_pos2 += m_pos1;
 18 
 19     m_velocity1 = Vec3(0.0f, 0.0f, 0.0f);
 20     m_velocity2 = Vec3(0.0f, 0.0f, 0.0f);
 21 
 22     UpdatePosition();
 23 }
 24 
 25 void        DoublePendulum02::Update(float t)
 26 {
 27     float dot;
 28 
 29     Vec3 line2 = m_pos2 - m_pos1;
 30     line2.Normalize();
 31 
 32     float xzsq = line2.x*line2.x + line2.z*line2.z;
 33     float xz = sqrtf(xzsq);
 34 
 35     // 下面的物體當前的加速度
 36     Vec3 dir2(-line2.x, xzsq/line2.y,-line2.z);
 37     if (line2.y > 0.0f)
 38     {
 39         dir2 = -dir2;
 40     }
 41     float d = dir2.Magnitude();
 42     if (d > 0.00001f)
 43     {
 44         dir2 /= d;
 45     }
 46     else
 47     {
 48         dir2 = m_velocity2;
 49         dir2.Normalize();
 50     }
 51     Vec3 acc2 = dir2*(m_fGravity*xz/m_l2);
 52 
 53     dot = dir2|m_velocity2;
 54     m_velocity2 = dir2*m_velocity2.Magnitude();
 55     if (dot < 0.0f)
 56     {
 57         m_velocity2 = -m_velocity2;
 58     }
 59     m_pos2 += m_velocity2*t + acc2*(0.5f*t*t);
 60     m_velocity2 += acc2*t;
 61 
 62     // 上面的物體受到下面繩子的拉力
 63     Vec3 force = line2*(-m_m2*m_fGravity*line2.y/m_l2);
 64     // 加上重力
 65     force.y -= m_m1*m_fGravity;
 66 
 67     // 上面的物體
 68     Vec3 line1 = m_pos1;
 69     line1.Normalize();
 70     xzsq = line1.x*line1.x + line1.z*line1.z;
 71     xz = sqrtf(xzsq);
 72     Vec3 dir1(-line1.x, xzsq/line1.y,-line1.z);
 73     if (line1.y > 0.0f)
 74     {
 75         dir1 = -dir1;
 76     }
 77     d = dir1.Magnitude();
 78     if (d > 0.00001f)
 79     {
 80         dir1 /= d;
 81     }
 82     else
 83     {
 84         dir1 = m_velocity1;
 85         dir1.Normalize();
 86     }
 87 
 88     dot = dir1|force;
 89     Vec3 acc1 = dir1*(dot/m_m1);
 90 
 91     dot = dir1|m_velocity1;
 92     m_velocity1 = dir1*m_velocity1.Magnitude();
 93     if (dot < 0.0f)
 94     {
 95         m_velocity1 = -m_velocity1;
 96     }
 97     m_pos1 += m_velocity1*t + acc1*(0.5f*t*t);
 98     m_velocity1 += acc1*t;
 99 
100     UpdatePosition();
101 
102 }

View Code

軟件截圖:

軟件下載地址：

http://files.cnblogs.com/files/WhyEngine/DoublePendulum.7z

使用說明:

程序啟動后,會出現三個隨機大小的球體在運動.

鼠標右鍵用於控制視角
鍵盤U用於開關UI用戶界面.
通過UI用戶界面可以設置兩個球體的質量,連線長度,設置引力系數,設置時間縮放速度,設置球體的顯示大小.

鍵盤1,2用於開關兩個球體運動軌跡的顯示

鍵盤G,用於開關三維網格的顯示
鍵盤C,用於開關坐標軸的顯示
鍵盤P,用於暫停
鍵盤R,用於重置,這時會隨機為兩個球體設置質量與初速度.

最后發一幅通過雙擺繪制的圖畫:

免責聲明！

本站轉載的文章為個人學習借鑒使用，本站對版權不負任何法律責任。如果侵犯了您的隱私權益，請聯系本站郵箱yoyou2525@163.com刪除。

猜您在找 uniapp -- 實現雙聯動微信小程序實現圖片雙滑縮放大小 Duilib的雙緩沖實現，附帶GDI、WTL的雙緩沖實現壓擺率 echarts實現雙y軸簡單雙序列比對實現 echarts 雙y軸實現 (原創)Notepad++怎么實現雙視圖/雙窗口？ centos7添加新網卡實現雙IP雙網關北上或者南下之被擺了一道