您的位置：首页 > 其它

双摆的程序实现

2016-07-29 00:00 387 查看

双摆是物理中的一个概念，先给下单摆与双摆的定义:

单摆：由一根不可伸长、质量不计的绳子，上端固定，下端系一个质点的装置。

双摆：是一个摆的支点装在另一摆的下部所形成的组合物体。双摆有两个摆角，所以有两个自由度。双摆是多自由度振动系统的最简单的力学模型之一，它也是一种混沌实例。

这里对双摆的实现程序与上一篇文章中三体的实现很相似，要实现它的程序需要一定的数学基础。

下面两个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中可以看到双摆的几个输入数据是：重力系数，两摆的距离和两球质量。而计算的数据是：每一个时刻两个球的位置。

下一步是对该基类进行继承实现。

.h

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的,双摆可以出现在三维空间中.我也试着写了下,使用力学来进行模拟,效果不太精确:

.h

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,用于重置,这时会随机为两个球体设置质量与初速度.

最后发一幅通过双摆绘制的图画: