双摆的程序实现

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

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

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

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

下面两个GIF动画图像为双摆的录屏:

      

      

先写了一个双摆的抽象基类：

// --------------------------------------------------------------------------------------

inline float     rand2(float a, float b)
{
    const float r = (float)(::rand() / ((float)RAND_MAX + 1));
    return r*(b-a) + a;
}

// --------------------------------------------------------------------------------------

class IDoublePendulum
{
public:

    IDoublePendulum()
    {
        m_fGravity = 9.8f;

        m_m1 = 1.0f;
        m_m2 = 2.0f;
    
        m_l1 = 1.0f;
        m_l2 = 2.0f;
    }

    virtual ~IDoublePendulum() {}

    // 重置
    virtual void        Reset()
    {
        m_m1 = rand2(1.0f, 5.0f);
        m_m2 = rand2(1.0f, 5.0f);

        m_l1 = rand2(1.0f, 5.0f);
        m_l2 = rand2(1.0f, 5.0f);
    }

    // 按时更新
    virtual void        Update(float deltaTime) = NULL;

    // 重力系数
    virtual void        SetGravity(float g)
    {
        m_fGravity = g;
    }

    // 质量
    virtual void        SetMass1(float m)
    {
        m_m1 = m;
    }

    virtual void        SetMass2(float m)
    {
        m_m2 = m;
    }

    // 长度
    virtual void        SetLength1(float l)
    {
        m_l1 = l;
        UpdatePosition();
    }

    virtual void        SetLength2(float l)
    {
        m_l2 = l;
        UpdatePosition();
    }

    float               GetGravity() const
    {
        return m_fGravity;
    }

    float               GetMass1() const
    {
        return m_m1;
    }

    float               GetMass2() const
    {
        return m_m2;
    }

    float               GetLength1() const
    {
        return m_l1;
    }

    float               GetLength2() const
    {
        return m_l2;
    }

    const Vec3&         GetPosition1() const
    {
        return m_pos1;
    }

    const Vec3&         GetPosition2() const
    {
        return m_pos2;
    }

protected:
    // 更新两球位置
    virtual void        UpdatePosition() = NULL;

protected:
    float m_fGravity;   // 重力系数

    float m_m1, m_m2;   // 两球质量
    float m_l1, m_l2;   // 两球距离

    Vec3 m_pos1, m_pos2;// 两球位置
};

从IDoublePendulum中可以看到双摆的几个输入数据是：重力系数，两摆的距离和两球质量。而计算的数据是：每一个时刻两个球的位置。

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

.h

// --------------------------------------------------------------------------------------

#ifndef _DoublePendulum01_H_
#define _DoublePendulum01_H_

// --------------------------------------------------------------------------------------

#include "IDoublePendulum.h"

// --------------------------------------------------------------------------------------

class DoublePendulum01 : public IDoublePendulum
{
public:
    DoublePendulum01()
    {
        m_a1 = 1.0f;
        m_a2 = 2.0f;    

        m_v1 = 0.0f;
        m_v2 = 0.0f;

        m_w1 = m_a1;
        m_w2 = m_a2;
    }

    // 重置
    void        Reset();

    // 按时更新
    void        Update(float deltaTime);

protected:
    void        UpdatePosition()
    {
        m_pos1.x = m_l1*sinf(m_w1);
        m_pos1.y = -m_l1*cosf(m_w1);
        m_pos1.z = 0.0f;

        m_pos2.x = m_pos1.x + m_l2*sinf(m_w2);
        m_pos2.y = m_pos1.y - m_l2*cosf(m_w2);
        m_pos2.z = 0.0f;
    }

private:
    float m_a1, m_a2;   // 两球初始角度
    float m_w1, m_w2;   // 两球当前角度
    float m_v1, m_v2;   // 两球的角加速度
};

// --------------------------------------------------------------------------------------

#endif

CPP

// --------------------------------------------------------------------------------------

#include "DoublePendulum01.h"

// --------------------------------------------------------------------------------------

/*
求解线型方程组
a*x + b*y + c = 0
d*x + e*y + f = 0
*/
inline bool SolveLinalg(float a, float b, float c, float d, float e, float f, float& x, float& y)
{
    float t = b*d - a*e;
    if (fabs(t) < FLT_EPSILON)
    {
        x = 0.0f;
        y = 0.0f;

        return false;
    }

    x = (c*e - b*f)/t;
    y = (a*f - c*d)/t;

    return true;
}

// --------------------------------------------------------------------------------------

void        DoublePendulum01::Reset()
{
    IDoublePendulum::Reset();

    m_a1 = rand2(-2.0f, 2.0f);
    m_a2 = rand2(-2.0f, 2.0f); 

    m_v1 = 0.0f;
    m_v2 = 0.0f;

    m_w1 = m_a1;
    m_w2 = m_a2;

    UpdatePosition();
}

void        DoublePendulum01::Update(float deltaTime)
{
    float a = m_l1*m_l1*(m_m1+m_m2);
    float b = m_l1*m_m2*m_l2*cos(m_w1-m_w2);
    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));

    float d = m_m2*m_l2*m_l1*cos(m_w1-m_w2);
    float e = m_m2*m_l2*m_l2;
    float f = m_m2*m_l2*(-m_l1*sin(m_w1-m_w2)*m_v1*m_v1 + m_fGravity*sin(m_w2));

    // 角加速度
    float dv1;
    float dv2;
    SolveLinalg(a, b, c, d, e, f, dv1, dv2);

    // 角速度
    m_v1 += dv1*deltaTime;
    m_v2 += dv2*deltaTime;

    // 角度
    m_w1 += m_v1*deltaTime;
    m_w2 += m_v2*deltaTime;

    UpdatePosition();
}

这里使用的是角度变化实现双摆，参考的资料是：

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

它的理论有些难，我看得也似懂非懂的。好在公式就在那里，只要会用，不求会懂，程序就能实现。

这个程序是在2D的,双摆可以出现在三维空间中.我也试着写了下,使用力学来进行模拟,效果不太精确:

.h

// --------------------------------------------------------------------------------------

#ifndef _DoublePendulum02_H_
#define _DoublePendulum02_H_

// --------------------------------------------------------------------------------------

#include "IDoublePendulum.h"

// --------------------------------------------------------------------------------------

class DoublePendulum02 : public IDoublePendulum
{
public:
    DoublePendulum02()
    {
        m_velocity1 = Vec3(0.0f, 0.0f, 0.0f);
        m_velocity2 = Vec3(0.0f, 0.0f, 0.0f);
    }

    // 重置
    void        Reset();

    // 按时更新
    void        Update(float deltaTime);

protected:
    void        UpdatePosition()
    {
        Vec3 v = m_pos2 - m_pos1;
        v.Normalize();

        m_pos1.Normalize();
        m_pos1 *= m_l1;

        m_pos2 = m_pos1 + v*m_l2;
    }

private:
    Vec3 m_velocity1, m_velocity2;  // 两球当前速度
};

// --------------------------------------------------------------------------------------

#endif

.cpp

// --------------------------------------------------------------------------------------

#include "DoublePendulum02.h"

// --------------------------------------------------------------------------------------

void        DoublePendulum02::Reset()
{
    IDoublePendulum::Reset();

    m_pos1.x = rand2(-1.0f, 1.0f);
    m_pos1.y = rand2(-1.0f, 0.0f);
    m_pos1.z = rand2(-1.0f, 1.0f);
    m_pos2.x = rand2(-1.0f, 1.0f);
    m_pos2.y = rand2(-0.5f, 0.5f);
    m_pos2.z = rand2(-1.0f, 1.0f);
    m_pos2 += m_pos1;

    m_velocity1 = Vec3(0.0f, 0.0f, 0.0f);
    m_velocity2 = Vec3(0.0f, 0.0f, 0.0f);

    UpdatePosition();
}

void        DoublePendulum02::Update(float t)
{
    float dot;

    Vec3 line2 = m_pos2 - m_pos1;
    line2.Normalize();

    float xzsq = line2.x*line2.x + line2.z*line2.z;
    float xz = sqrtf(xzsq);

    // 下面的物体当前的加速度
    Vec3 dir2(-line2.x, xzsq/line2.y,-line2.z);
    if (line2.y > 0.0f)
    {
        dir2 = -dir2;
    }
    float d = dir2.Magnitude();
    if (d > 0.00001f)
    {
        dir2 /= d;
    }
    else
    {
        dir2 = m_velocity2;
        dir2.Normalize();
    }
    Vec3 acc2 = dir2*(m_fGravity*xz/m_l2);

    dot = dir2|m_velocity2;
    m_velocity2 = dir2*m_velocity2.Magnitude();
    if (dot < 0.0f)
    {
        m_velocity2 = -m_velocity2;
    }
    m_pos2 += m_velocity2*t + acc2*(0.5f*t*t);
    m_velocity2 += acc2*t;

    // 上面的物体受到下面绳子的拉力
    Vec3 force = line2*(-m_m2*m_fGravity*line2.y/m_l2);
    // 加上重力
    force.y -= m_m1*m_fGravity;

    // 上面的物体
    Vec3 line1 = m_pos1;
    line1.Normalize();
    xzsq = line1.x*line1.x + line1.z*line1.z;
    xz = sqrtf(xzsq);
    Vec3 dir1(-line1.x, xzsq/line1.y,-line1.z);
    if (line1.y > 0.0f)
    {
        dir1 = -dir1;
    }
    d = dir1.Magnitude();
    if (d > 0.00001f)
    {
        dir1 /= d;
    }
    else
    {
        dir1 = m_velocity1;
        dir1.Normalize();
    }

    dot = dir1|force;
    Vec3 acc1 = dir1*(dot/m_m1);

    dot = dir1|m_velocity1;
    m_velocity1 = dir1*m_velocity1.Magnitude();
    if (dot < 0.0f)
    {
        m_velocity1 = -m_velocity1;
    }
    m_pos1 += m_velocity1*t + acc1*(0.5f*t*t);
    m_velocity1 += acc1*t;

    UpdatePosition();

}

软件截图:

软件下载地址：

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

使用说明:

程序启动后,会出现三个随机大小的球体在运动.

鼠标右键用于控制视角
键盘U用于开关UI用户界面.
通过UI用户界面可以设置两个球体的质量,连线长度,设置引力系数,设置时间缩放速度,设置球体的显示大小.

键盘1,2用于开关两个球体运动轨迹的显示

键盘G,用于开关三维网格的显示
键盘C,用于开关坐标轴的显示
键盘P,用于暂停
键盘R,用于重置,这时会随机为两个球体设置质量与初速度.

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

相关文章:三体三体