背景
傳統(tǒng)的雙擺系統(tǒng)通常由兩個長度分別為l1 ，l2 的細桿和兩個固定在細桿末端，質(zhì)量為m1 ，m2 的小球組成。內(nèi)擺和垂直線之間的夾角為θ1 。外擺與垂直線之間的夾角為θ2 。細桿的質(zhì)量和球的形狀對于不同的研究人員有所不同。本文研究理想狀態(tài)，即系桿的質(zhì)量、系統(tǒng)阻尼和阻力忽略不計，小球看做質(zhì)點。雙擺系統(tǒng)簡圖如圖1所示。
圖1 雙擺結(jié)構(gòu)簡圖
模型建立
此處模型建立利用拉格朗日力學(xué)原理，拉格朗日力學(xué)（Lagrangian mechanics）是分析力學(xué)中的一種，于1788年由約瑟夫·拉格朗日所創(chuàng)立。拉格朗日力學(xué)是對經(jīng)典力學(xué)的一種的新的理論表述，著重于數(shù)學(xué)解析的方法，并運用最小作用量原理，是分析力學(xué)的重要組成部分。經(jīng)典力學(xué)最初的表述形式由牛頓建立，它著重于分析位移，速度，加速度，力等矢量間的關(guān)系，又稱為矢量力學(xué)。拉格朗日引入了廣義坐標的概念，又運用達朗貝爾原理，求得與牛頓第二定律等價的拉格朗日方程。不僅如此，拉格朗日方程具有更普遍的意義，適用范圍更廣泛。還有，選取恰當(dāng)?shù)膹V義坐標，可以大大地簡化拉格朗日方程的求解過程。
如圖2所示，建立直角坐標系：
圖2 建立直角坐標系
求解方法
最后得到的表達式很長，就不展示了。
上面的計算均在Matlab上完成，因為其計算很強大，下面就開始用JavaScript來實現(xiàn)計算。
采用四階龍格-庫塔算法進行迭代計算：
實現(xiàn)的關(guān)鍵代碼如下：
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; //這里最終得到原始數(shù)據(jù)
this.omega1 = y2;
this.theta2 = y3;
this.omega2 = y4;
// 由于原始數(shù)據(jù)過度密集，導(dǎo)致繪制曲線困難，下面將繪制曲線的數(shù)據(jù)進行稀釋，加速繪圖速度。
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])];
};
可能排版很亂，很多符號實在不好打，就截圖，請諒解。
在線體驗連接：點擊在線體驗
