二次元の調和振動子

（教科書:線形・非線形力学とカオスへの入門(７ページから)を参照）
以下、２次元の調和振動子の運動をMathematicaを使って調べてみよう．微分方程式の解き方、初期条件の入力の仕方、リサージュ図形の描き方などを学ぶ．

等方的な場合

運動方程式がF=- k rとかける場合：ここでk はバネ定数でrは距離を表す．ω=(Sqrt[(k/m)])とすれば、

solution=DSolve[{x''[t]==-(ω^2) x[t],y''[t]==-(ω^2) y[t]},{x[t],y[t]},t]

となり、解が求まる．二階の微分方程式で２次元であるから四つの未知数が表れる．初期条件を時間t=0 の時、
x[0]==A Cos[α],  y[0]== BCos[β], 
x'[0]==-A ω Sin[α],
y'[0]==-B ω Sin[β]
であるとすると、（この初期条件は以下の計算において式が簡単になるように選んだ．このように選んでも一般性が失われないことに注意）

x[t_]=x[t]/.solution[[1]]

y[t_]=y[t]/.solution[[1]]

Solve[{x[0]==A Cos[α],y[0]==B Cos[β],x'[0]==-A ω Sin[α],y'[0]==-B ω Sin[β]},{C[1], C[2], C[3], C[4]}]

x[t], y[t]を以下のように定義してリサージュ図形を描いてみよう．

{x[t_], y[t_]}={x[t], y[t]}/.First[%]

または，

TrigFactor[%]

補足）もちろん上のような面倒なことをしなくても直接DSolveにおいて初期条件を入れてもよい．
DSolve[{x''[t]==-(ω^2) x[t],y''[t]==-(ω^2) y[t],x[0]==A Cos[α],y[0]==B Cos[β], x'[0]==-A ω Sin[α], y'[0]==-B ω Sin[β] },{x[t],y[t]},t]

まず、A=B=1, α=β=0, ω=1の場合

graph =ParametricPlot[{x[t], y[t]}/.{A->1,B->1,α->0,β->Pi/2,ω->1}, {t, 0, 10}]

となる．（/.は置き換えを意味する．）

実際、振動子がどのような運動をするのか、アニメーションを作ってみよう．

Manipulate[Show[graph,Graphics[{{Hue[1],PointSize[.1],Point[Evaluate[{x[t], y[t]}/.{A->1,B->1,α->0,β->Pi/2,ω->1}]]}}],PlotRange->{{-1.5,1.5},{-1.5,1.5}},AspectRatio->1],{t,0,2 Pi,.1}]

アニメーションのスピードやスッテップ間隔を変えてみよう．

問１以下のMathematica を使って微分方程式を解け．

ClearAll[x,y]

(例) ((dy(x))/(dx))==x (e^(-x))

DSolve[D[y[x],x]==x Exp[-x],y[x],x]

(1)((dy(x))/(dx))==sin(x) (e^(-x))

(2)((dy(x))/(dx))==(sin^2)(x) (e^(-x))

(3)(((d^2)y(x))/(d(x^2)))==sin(x) (e^(-x))

(4)(((d^2)y(x))/(d(x^2)))==(sin^2)(x) (e^(-x))

(5)((dy(x))/(dx))==(x^2) cos(x)

(6)(((d^2)y(x))/(d(x^2)))==(x^2) cos(x)

問2 以下の質問に実際にMathematica を使ってグラフを描き答えよ．

{x[t_], y[t_]}={A Cos[α+t ω], B Cos[β+t ω]}

とする．

(1) {A->1,B->1,α->0,β->π/2,ω->1} のときどうのようなリサージュ図形がえられるか．

(2) {A->1,B->1,α->0,β->π,ω->1} のときどうのようなリサージュ図形がえられるか．

(3) {A->1,B->1,α->0,β->π/4,ω->1} のときどうのようなリサージュ図形がえられるか．

(4) {A->1,B->1,α->0,β->π/4,ω->1} のときの調和振動子のアニメーションをつくってみよ．（解答用紙にはできたかどうかを書けばよい．）

(5) {A->1,B->1,α->0,β->π/2,ω->1} と {A->1,B->1,α->0,β->3π/2,ω->1} の のときの調和振動子の運動では何が違うか．

非等方的な場合

運動方程式がF=- (({{(Subscript[k, 1]), 0}, {0, (Subscript[k, 2])}})) rとかける場合：
(Subscript[ω, 1])=(Sqrt[((Subscript[k, 1])/m)]),(Subscript[ω, 2])=(Sqrt[((Subscript[k, 2])/m)])とすれば解は

{x[t_], y[t_]}={A Cos[α+t (Subscript[ω, 1])], B Cos[β+t (Subscript[ω, 2])]}

となる．いろいろなパラメター値にたいするリサージュ図を描いてみよう．

ParametricPlot[Evaluate[{x[t], y[t]}/.{A->1,B->1,α->0,β->0,(Subscript[ω, 1])->1,(Subscript[ω, 2])->2}], {t, 0, 10}]

問３　上の図でβ->(π/2)とするとどのような図形がえられるか．

角振動数の比が無理数の場合は軌跡は閉じない．例えば

ParametricPlot[Evaluate[{x[t], y[t]}/.{A->1,B->1,α->0,β->0,(Subscript[ω, 1])->1,(Subscript[ω, 2])->(Sqrt[2])}], {t, 0, 50}]

となるが,

問４　上の図で時間tをもっと大きくとるとどのような図が得られるか、結果を解釈せよ．

となるが,

問５　（optional) 教科書、１５ページにある図３において、初期角θ=60°で静止していた振り子が描く角度と各速度の相空間図を描け（２３ページの図を参照）．