背景
传统的双摆系统通常由两个长度分别为l1 ，l2 的细杆和两个固定在细杆末端，质量为m1 ，m2 的小球组成。内摆和垂直线之间的夹角为θ1 。外摆与垂直线之间的夹角为θ2 。细杆的质量和球的形状对于不同的研究人员有所不同。本文研究理想状态，即系杆的质量、系统阻尼和阻力忽略不计，小球看做质点。双摆系统简图如图1所示。
图1 双摆结构简图
模型建立
此处模型建立利用拉格朗日力学原理，拉格朗日力学（Lagrangian mechanics）是分析力学中的一种，于1788年由约瑟夫·拉格朗日所创立。拉格朗日力学是对经典力学的一种的新的理论表述，着重于数学解析的方法，并运用最小作用量原理，是分析力学的重要组成部分。经典力学最初的表述形式由牛顿建立，它着重于分析位移，速度，加速度，力等矢量间的关系，又称为矢量力学。拉格朗日引入了广义坐标的概念，又运用达朗贝尔原理，求得与牛顿第二定律等价的拉格朗日方程。不仅如此，拉格朗日方程具有更普遍的意义，适用范围更广泛。还有，选取恰当的广义坐标，可以大大地简化拉格朗日方程的求解过程。
如图2所示，建立直角坐标系：
图2 建立直角坐标系
求解方法
最后得到的表达式很长，就不展示了。
上面的计算均在Matlab上完成，因为其计算很强大，下面就开始用JavaScript来实现计算。
采用四阶龙格-库塔算法进行迭代计算：
实现的关键代码如下：
var m1 = this.m1; //带this的是与输入绑定的变量
var m2 = this.m2;
var l1 = this.l1;
var l2 = this.l2;
var g = this.g;
var y1 = [];
var y2 = [];
var y3 = [];
var y4 = [];
y1[0] = [0, this.initialTheta1 / 180 * Math.PI]; //迭代初始条件
y2[0] = [0, this.initialOmega1 / 180 * Math.PI];
y3[0] = [0, this.initialTheta2 / 180 * Math.PI];
y4[0] = [0, this.initialOmega2 / 180 * Math.PI];
var k1_1 = 0;
var k1_2 = 0;
var k1_3 = 0;
var k1_4 = 0;
var k2_1 = 0;
var k2_2 = 0;
var k2_3 = 0;
var k2_4 = 0;
var k3_1 = 0;
var k3_2 = 0;
var k3_3 = 0;
var k3_4 = 0;
var k4_1 = 0;
var k4_2 = 0;
var k4_3 = 0;
var k4_4 = 0;
var h = this.stepLength;
var timeRange = this.timeRange;
for (var i = 0; i < timeRange / h; i++) {
k1_1 = y2[i][1];
k1_2 = y2[i][1] + h / 2 * k1_1;
k1_3 = y2[i][1] + h / 2 * k1_2;
k1_4 = y2[i][1] + h * k1_3;
y1[i + 1] = [(i + 1) * h, y1[i][1] + h / 6 * (k1_1 + 2 * k1_2 + 2 * k1_3 + k1_4)];
k2_1 = (m2 * Math.cos(y1[i][1] - y3[i][1]) * (g * Math.sin(y3[i][1]) - l1 *
Math.sin(y1[i][1] - y3[i][1]) * Math.pow(y2[i][1], 2)) / (l1 * m1 + l1 *
m2 - l1 * m2 * Math.pow(Math.cos(y1[i][1] - y3[i][1]), 2)) - (g * Math.sin(
y1[i][1]) * (m1 + m2) + l2 * m2 * Math.sin(y1[i][1] - y3[i][1]) * Math
.pow(y4[i][1], 2)) / (l1 * m1 + l1 * m2 - l1 * m2 * Math.pow(Math.cos(
y1[i][1] - y3[i][1]), 2)));
k2_2 = (m2 * Math.cos((y1[i][1] + h / 2 * k2_1) - (y3[i][1] + h / 2 * k2_1)) * (g *
Math.sin(y3[i][1] + h / 2 * k2_1) - l1 * Math.sin((y1[i][1] + h / 2 * k2_1) - (y3[i][1] +
h / 2 * k2_1)) * Math.pow((y2[i][1] + h / 2 * k2_1), 2)) / (l1 * m1 + l1 * m2 -
l1 * m2 * Math.pow(Math.cos((y1[i][1] + h / 2 * k2_1) - (y3[i][1] + h / 2 * k2_1)), 2)
) - (g * Math.sin(y1[i][1] + h / 2 * k2_1) * (m1 + m2) + l2 * m2 * Math.sin((
y1[i][1] + h / 2 * k2_1) - (y3[i][1] + h / 2 * k2_1)) * Math.pow((y4[i][1] + h / 2 * k2_1),
2)) / (l1 * m1 + l1 * m2 - l1 * m2 * Math.pow(Math.cos((y1[i][1] + h /
2 * k2_1) - (y3[i][1] + h / 2 * k2_1)), 2)));
k2_3 = (m2 * Math.cos((y1[i][1] + h / 2 * k2_2) - (y3[i][1] + h / 2 * k2_2)) * (g *
Math.sin(y3[i][1] + h / 2 * k2_2) - l1 * Math.sin((y1[i][1] + h / 2 * k2_2) - (y3[i][1] +
h / 2 * k2_2)) * Math.pow((y2[i][1] + h / 2 * k2_2), 2)) / (l1 * m1 + l1 * m2 -
l1 * m2 * Math.pow(Math.cos((y1[i][1] + h / 2 * k2_2) - (y3[i][1] + h / 2 * k2_2)), 2)
) - (g * Math.sin(y1[i][1] + h / 2 * k2_2) * (m1 + m2) + l2 * m2 * Math.sin((
y1[i][1] + h / 2 * k2_2) - (y3[i][1] + h / 2 * k2_2)) * Math.pow((y4[i][1] + h / 2 * k2_2),
2)) / (l1 * m1 + l1 * m2 - l1 * m2 * Math.pow(Math.cos((y1[i][1] + h /
2 * k2_2) - (y3[i][1] + h / 2 * k2_2)), 2)));
k2_4 = (m2 * Math.cos((y1[i][1] + h * k2_3) - (y3[i][1] + h * k2_3)) * (g * Math.sin(
y3[i][1] + h * k2_3) - l1 * Math.sin((y1[i][1] + h * k2_3) - (y3[i][1] + h * k2_3)) *
Math.pow((y2[i][1] + h * k2_3), 2)) / (l1 * m1 + l1 * m2 - l1 * m2 *
Math.pow(Math.cos((y1[i][1] + h * k2_3) - (y3[i][1] + h * k2_3)), 2)) - (g * Math.sin(y1[
i][1] + h * k2_3) * (m1 + m2) + l2 * m2 * Math.sin((y1[i][1] + h * k2_3) - (
y3[i][1] + h * k2_3)) * Math.pow((y4[i][1] + h * k2_3), 2)) / (l1 * m1 + l1 * m2 -
l1 * m2 * Math.pow(Math.cos((y1[i][1] + h * k2_3) - (y3[i][1] + h * k2_3)), 2)));
y2[i + 1] = [(i + 1) * h, y2[i][1] + h / 6 * (k2_1 + 2 * k2_2 + 2 * k2_3 + k2_4)];
k3_1 = y4[i][1];
k3_2 = y4[i][1] + h / 2 * k3_1;
k3_3 = y4[i][1] + h / 2 * k3_2;
k3_4 = y4[i][1] + h * k3_3;
y3[i + 1] = [(i + 1) * h, y3[i][1] + h / 6 * (k3_1 + 2 * k3_2 + 2 * k3_3 + k3_4)];
k4_1 = (Math.cos(y1[i][1] - y3[i][1]) * (g * Math.sin(y1[i][1]) * (m1 + m2) +
l2 * m2 * Math.sin(y1[i][1] - y3[i][1]) * Math.pow(y4[i][1], 2))) / (l2 *
m1 + l2 * m2 - l2 * m2 * Math.pow(Math.cos(y1[i][1] - y3[i][1]), 2)) - (
(m1 + m2) * (g * Math.sin(y3[i][1]) - l1 * Math.sin(y1[i][1] - y3[i][1]) *
Math.pow(y2[i][1], 2))) / (l2 * m1 + l2 * m2 - l2 * m2 * Math.pow(Math.cos(
y1[i][1] - y3[i][1]), 2));
k4_2 = (Math.cos((y1[i][1] + h / 2 * k4_1) - (y3[i][1] + h / 2 * k4_1)) * (g * Math.sin((
y1[i][1] + h / 2 * k4_1)) * (m1 + m2) + l2 * m2 * Math.sin((y1[i][1] + h /
2 * k4_1) - (y3[i][1] + h / 2 * k4_1)) * Math.pow((y4[i][1] + h / 2 * k4_1), 2))) / (l2 *
m1 + l2 * m2 - l2 * m2 * Math.pow(Math.cos((y1[i][1] + h / 2 * k4_1) - (y3[
i][1] + h / 2 * k4_1)), 2)) - ((m1 + m2) * (g * Math.sin((y3[i][1] + h / 2 * k4_1)) -
l1 * Math.sin((y1[i][1] + h / 2 * k4_1) - (y3[i][1] + h / 2 * k4_1)) * Math.pow((y2[i]
[1] + h / 2 * k4_1), 2))) / (l2 * m1 + l2 * m2 - l2 * m2 * Math.pow(Math.cos(
(y1[i][1] + h / 2 * k4_1) - (y3[i][1] + h / 2 * k4_1)), 2));
k4_3 = (Math.cos((y1[i][1] + h / 2 * k4_2) - (y3[i][1] + h / 2 * k4_2)) * (g * Math.sin((
y1[i][1] + h / 2 * k4_2)) * (m1 + m2) + l2 * m2 * Math.sin((y1[i][1] + h /
2 * k4_2) - (y3[i][1] + h / 2 * k4_2)) * Math.pow((y4[i][1] + h / 2 * k4_2), 2))) / (l2 *
m1 + l2 * m2 - l2 * m2 * Math.pow(Math.cos((y1[i][1] + h / 2 * k4_2) - (y3[
i][1] + h / 2 * k4_2)), 2)) - ((m1 + m2) * (g * Math.sin((y3[i][1] + h / 2 * k4_2)) -
l1 * Math.sin((y1[i][1] + h / 2 * k4_2) - (y3[i][1] + h / 2 * k4_2)) * Math.pow((y2[i]
[1] + h / 2 * k4_2), 2))) / (l2 * m1 + l2 * m2 - l2 * m2 * Math.pow(Math.cos(
(y1[i][1] + h / 2 * k4_2) - (y3[i][1] + h / 2 * k4_2)), 2));
k4_4 = (Math.cos((y1[i][1] + h * k4_3) - (y3[i][1] + h * k4_3)) * (g * Math.sin((y1[i]
[1] + h * k4_3)) * (m1 + m2) + l2 * m2 * Math.sin((y1[i][1] + h * k4_3) - (y3[
i][1] + h * k4_3)) * Math.pow((y4[i][1] + h * k4_3), 2))) / (l2 * m1 + l2 * m2 -
l2 * m2 * Math.pow(Math.cos((y1[i][1] + h * k4_3) - (y3[i][1] + h * k4_3)), 2)) - ((
m1 + m2) * (g * Math.sin((y3[i][1] + h * k4_3)) - l1 * Math.sin((y1[i][1] +
h * k4_3) - (y3[i][1] + h * k4_3)) * Math.pow((y2[i][1] + h * k4_3), 2))) / (l2 * m1 +
l2 * m2 - l2 * m2 * Math.pow(Math.cos((y1[i][1] + h * k4_3) - (y3[i][1] + h *
k4_3)), 2));
y4[i + 1] = [(i + 1) * h, y4[i][1] + h / 6 * (k4_1 + 2 * k4_2 + 2 * k4_3 + k4_4)];
};
this.theta1 = y1; //这里最终得到原始数据
this.omega1 = y2;
this.theta2 = y3;
this.omega2 = y4;
// 由于原始数据过度密集，导致绘制曲线困难，下面将绘制曲线的数据进行稀释，加速绘图速度。
for(var j = 0; j < this.timeRange/this.stepLength; j += 100){
this.roughTheta1.push([j*this.stepLength, this.theta1[j][1]]);
this.roughOmega1.push([j*this.stepLength, this.omega1[j][1]]);
this.roughTheta2.push([j*this.stepLength, this.theta2[j][1]]);
this.roughOmega2.push([j*this.stepLength, this.omega2[j][1]]);
this.innerBallPosition[j/100] = [l1 * Math.sin(y1[j][1]), -l1 * Math.cos(y1[j][1])];
this.outerBallPosition[j/100] = [l1 * Math.sin(y1[j][1]) + l2 * Math.sin(y3[j][1]), -l1 * Math.cos(y1[j][1]) - l2 * Math.cos(y3[j][1])];
};
可能排版很乱，很多符号实在不好打，就截图，请谅解。
在线体验连接：点击在线体验
