双摆是物理中的一个概念，先给下单摆与双摆的定义:
单摆：由一根不可伸长、质量不计的绳子，上端固定，下端系一个质点的装置。
双摆：是一个摆的支点装在另一摆的下部所形成的组合物体。双摆有两个摆角，所以有两个自由度。双摆是多自由度振动系统的最简单的力学模型之一，它也是一种混沌实例。
这里对双摆的实现程序与上一篇文章中三体的实现很相似，要实现它的程序需要一定的数学基础。
下面两个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