背景
傳統(tǒng)的雙擺系統(tǒng)通常由兩個(gè)長(zhǎng)度分別為l1 ，l2 的細(xì)桿和兩個(gè)固定在細(xì)桿末端怀跛，質(zhì)量為m1 拣宰，m2 的小球組成轩端。內(nèi)擺和垂直線之間的夾角為θ1 放典。外擺與垂直線之間的夾角為θ2 。細(xì)桿的質(zhì)量和球的形狀對(duì)于不同的研究人員有所不同基茵。本文研究理想狀態(tài)奋构，即系桿的質(zhì)量、系統(tǒng)阻尼和阻力忽略不計(jì)拱层，小球看做質(zhì)點(diǎn)弥臼。雙擺系統(tǒng)簡(jiǎn)圖如圖1所示。
圖1 雙擺結(jié)構(gòu)簡(jiǎn)圖
模型建立
此處模型建立利用拉格朗日力學(xué)原理根灯，拉格朗日力學(xué)（Lagrangian mechanics）是分析力學(xué)中的一種径缅，于1788年由約瑟夫·拉格朗日所創(chuàng)立。拉格朗日力學(xué)是對(duì)經(jīng)典力學(xué)的一種的新的理論表述烙肺，著重于數(shù)學(xué)解析的方法纳猪，并運(yùn)用最小作用量原理，是分析力學(xué)的重要組成部分桃笙。經(jīng)典力學(xué)最初的表述形式由牛頓建立氏堤，它著重于分析位移，速度搏明，加速度鼠锈，力等矢量間的關(guān)系，又稱(chēng)為矢量力學(xué)星著。拉格朗日引入了廣義坐標(biāo)的概念购笆，又運(yùn)用達(dá)朗貝爾原理，求得與牛頓第二定律等價(jià)的拉格朗日方程虚循。不僅如此同欠，拉格朗日方程具有更普遍的意義，適用范圍更廣泛邮丰。還有行您，選取恰當(dāng)?shù)膹V義坐標(biāo)，可以大大地簡(jiǎn)化拉格朗日方程的求解過(guò)程剪廉。
如圖2所示娃循，建立直角坐標(biāo)系：
圖2 建立直角坐標(biāo)系
求解方法
最后得到的表達(dá)式很長(zhǎng)，就不展示了斗蒋。
上面的計(jì)算均在Matlab上完成捌斧，因?yàn)槠溆?jì)算很強(qiáng)大笛质，下面就開(kāi)始用JavaScript來(lái)實(shí)現(xiàn)計(jì)算。
采用四階龍格-庫(kù)塔算法進(jìn)行迭代計(jì)算：
實(shí)現(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ù)過(guò)度密集捞蚂，導(dǎo)致繪制曲線困難妇押，下面將繪制曲線的數(shù)據(jù)進(jìn)行稀釋?zhuān)铀倮L圖速度。
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])];
};
可能排版很亂姓迅，很多符號(hào)實(shí)在不好打敲霍，就截圖，請(qǐng)諒解丁存。
在線體驗(yàn)連接：點(diǎn)擊在線體驗(yàn)
