二重振り子

はじめに

二重振り子の運動は振幅が十分小さいときには線形な微分方程式で記述できるが、一般には複雑な非線形運動を行う．この系は散逸系ではないカオスを出現できる系として知られている．ここではその運動を記述する微分方程式を導きだし、運動のアニメーションをつくることにする．

ClearAll["Global`*"]

次の様な座標とパラメターを下図のように設定する．

x1=0.426024;y1=-0.904712;x2=1.40871;y2=-0.719444;
Show[Graphics[{Hue[.8],Line[{{x1,y1}, {x2,y2}}]}],Graphics[{Hue[.6],PointSize[.06],Point[{x2,y2}]}],Graphics[{Hue[.8],Line[{{0,0}, {x1,y1}}]}],Graphics[{Hue[1],PointSize[.06],Point[{x1,y1}]}],Graphics[{Text["(Subscript[m, 1])",{0.417078,-0.881541}],Text["(Subscript[m, 2])",{1.39758,-0.697567}],{Hue[.4],Line[{{x1,-2},{x1,2}}]},Text[(Subscript[θ, 1]),{ .1,-0.5}],Text[(Subscript[θ, 2]),{0.59269,-1.11151}],Text["(Subscript[l, 1])",{0.329272,-0.344951}],Text["(Subscript[l, 2])",{0.841473,-0.666905}]}],Axes->True,Ticks->None,AspectRatio->1,PlotRange->{{-2.1,2.1},{-2.2,2.2}}]

それぞれの質量の位置を表す座標を次のようにする．

ClearAll[θ1,θ2]

x1=l1 Sin[θ1[t]];y1=-(l1 Cos[θ1[t]]);
x2=l1 Sin[θ1[t]]+l2 Sin[θ2[t]];
y2=-(l1 Cos[θ1[t]])-l2 Cos[θ2[t]];

次にラグランジェアン(運動エネルギーと位置エネルギーの差)を求める．

Lagrange=(1/2) m1 (((((∂_t)x1))^2)+((((∂_t)y1))^2))+(1/2) m2 (((((∂_t)x2))^2)+((((∂_t)y2))^2))-(m1 g y1+m2 g y2)

ここで、次のパラメターがコンスタントであることを宣言しておく．

SetAttributes[{l1, l2, m1, m2, g}, Constant]

次の計算でそれぞれの結果(eq1とeq2 )をゼロとおけば、系の運動方程式が得られる．

eq1=Simplify[Dt[(∂_((θ1)'[t]))Lagrange,t]-(∂_(θ1[t]))Lagrange]

eq2=Simplify[Dt[(∂_((θ2)'[t]))Lagrange,t]-(∂_(θ2[t]))Lagrange]

運動は次の微分方程式に従う．

l2 m2 sin(θ1(t)-θ2(t)) (((θ2)'(t))^2)+g m1 sin(θ1(t))+g m2 sin(θ1(t))+l1 (m1+m2) (θ1)''(t)+l2 m2 cos(θ1(t)-θ2(t)) (θ2)''(t)==0

-l1 sin(θ1(t)-θ2(t)) (((θ1)'(t))^2)+g sin(θ2(t))+l1 cos(θ1(t)-θ2(t)) (θ1)''(t)+l2 (θ2)''(t)==0

実際の計算

ここで数値計算をすることにする．次のようにNDSolveを用いればいろいろなパラメターや初期値に対応する解が得られる．

{l1,l2,m1,m2,g}={1,1,1,1.0,9.8};
Timing[result=NDSolve[{g m1 Sin[θ1[t]]+g m2 Sin[θ1[t]]+l2 m2 Sin[θ1[t]-θ2[t]] (((θ2)'[t])^2)+l1 (m1+m2) (θ1)''[t]+l2 m2 Cos[θ1[t]-θ2[t]] (θ2)''[t]==0,g Sin[θ2[t]]-l1 Sin[θ1[t]-θ2[t]] (((θ1)'[t])^2)+l1 Cos[θ1[t]-θ2[t]] (θ1)''[t]+l2 (θ2)''[t]==0,θ1[0]==π,(θ1)'[0]==-0.1,θ2[0]==π,(θ2)'[0]==0.11},{θ1[t],θ2[t]},{t,0,40},MaxSteps->9000]]

θ1[t_]=θ1[t]/.First[result];
θ2[t_]=θ2[t]/.First[result];

Manipulate[Evaluate[Show[Graphics[{Hue[.8],Thick,Line[{{x1,y1}, {x2,y2}}]}],Graphics[{Hue[.6],PointSize[.06],Point[{x2,y2}]}],Graphics[{Hue[.8],Thick,Line[{{0,0}, {x1,y1}}]}],Graphics[{Hue[1],PointSize[.06],Point[{x1,y1}]}],Axes->True,Ticks->None,AspectRatio->1,PlotRange->{{-2.1,2.1},{-2.2,2.2}}]],{t,0,40,0.005}]

list=Table[Evaluate[Show[Graphics[{Hue[.8],Thick,Line[{{x1,y1}, {x2,y2}}]}],Graphics[{Hue[.6],PointSize[.06],Point[{x2,y2}]}],Graphics[{Hue[.8],Thick,Line[{{0,0}, {x1,y1}}]}],Graphics[{Hue[1],PointSize[.06],Point[{x1,y1}]}],Axes->True,Ticks->None,AspectRatio->1,PlotRange->{{-2.1,2.1},{-2.2,2.2}}]],{t,0,20,0.1}];

Export["doublependulum.gif",list]