matlab双摆程序,双摆运动仿真背后的原理

背景

传统的双摆系统通常由两个长度分别为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])];

};

可能排版很乱，很多符号实在不好打，就截图，请谅解。

在线体验连接：点击在线体验