振幅の大きな振り子

運動方程式

上図のようにパラメターを設定する．運動方程式は、

m L (((d^2)θ)/(d(t^2)))==-m g sin(θ)

と書ける．（教科書を参考）

問１　上の方程式をニュートンの第二法則から導き出せ．

θ が十分小さい

まず、θ が十分小さいとして、この方程式を解いてみよう．

Series[Sin[θ], {θ, 0, 10}]

であるから、第一項だけをとり微分方程式を解く．ただし、初期条件としてt=0のときθ[0]==Θ, θ'[0]==0であるとする．

DSolve[{ D[θ[t], {t, 2}]== - g/L θ[t],θ[0]==Θ,θ'[0]==0}, θ[t], t]

問２
（１）上の解より、振りの角度θが小さい時、振り子はどのような運動をするといえるか．
（２）また、その周期は何か．

問３　Taylor 展開は上の例で分かるようにSeriesを用いればよい．例えば、

Series[f[t], {t, 0,10}]

上の例では0の周りで展開したが、x=1の周りで展開したいときは、

Series[Log[x], {x,1,5}]

とすればよい．

次の関数を0の周りでx^10のオーダーまで求めよ・
1.  log(1-x)
2.  cos(x)
3.  tan(x)
4 . (e^x)
5. tanh(x)
6. (((sin^2)(x))/(x^2))

θ が一般の場合

運動の周期に関する議論は教科書をみてください．ここでは微分方程式を数値的に解いて、相空間図を描くことにする．まず

g=9.8;
L=1.0;

と仮定し（もちろん、数値計算では全てのパラメターに数字を代入しなければならない）、

Θ=Pi/3;

まずΘ=π/3として計算してみよう．

Θ=Pi/3

solution=NDSolve[{ D[θ[t], {t, 2}]== - g/L Sin[θ[t]],θ[0]==Θ,θ'[0]==0}, θ[t],{t, 0, 10}]

注意：

?NDSolve

角度と振幅の図を描く．

Plot[Evaluate[θ[t]/.solution[[1]]], {t, 0, 10} , AxesLabel->{angle, amplitude}]

角速度と振幅の図を描く．

Plot[Evaluate[D[θ[t]/.solution[[1]], t]], {t, 0, 10}, AxesLabel->{"angular velocity", amplitude}]

角度と角速度の相空間図を描きg1と名付ける．

g1=ParametricPlot[Evaluate[{θ[t]/.solution[[1]],D[θ[t]/.solution[[1]],t]}], {t, 0, 10}, PlotStyle->Hue[1]]

この結果と上の単振動の結果を相空間でどのように違うか、見てみよう．同じ初期条件でプロットすると、

g2=ParametricPlot[Evaluate[{Θ Cos[(((Sqrt[g]) t)/(Sqrt[L]))], D[Θ Cos[(((Sqrt[g]) t)/(Sqrt[L]))], t]}], {t, 0, 10}, PlotStyle->Hue[.6]]

この二つの結果を一つのグラフにする．

Show[g1, g2,AxesLabel->{"θ","((\[DifferentialD]θ)/(\[DifferentialD]t))" }]

問４上のグラフからどのようなことがいえるか．(なぜ上の二つの軌道は一致しないのか？）

問５　　NDSolve を使って以下の微分方程式を解け．
例

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

を

y(0)==1

y'(0)==1

として、0<x<10 の間で解いてみよう．

solution=NDSolve[{y''[x]==-Sin[y[x]], y[0]==1, y'[0]==1},y[x], {x, 0, 10}]

この結果を使うには、例えば

Plot[y[x]/.solution, {x, 0, 10}]

とすれば、グラフが得られるし、x=0.5の時のyの値がほしければ

y[x]/.solution/.x->0.5

というようにできる．

(1)

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

を

y(0)==1

y'(0)==1

として、-3Pi<x<3Pi の間で数値的に解いてグラフを描け．（解答用紙にプリントアウトを切り取ってはっても良いし、スケッチ的に書いても良い．）

(2)

(((d^2)y(x))/(d(x^2)))==-cos(10 y(x))(e^x)

を

y(0)==1

y'(0)==1

として、0<x<2Pi の間で数値的に解いてグラフを描け．

Θ=Pi-.01;

Θ=πの位置は特異点である．その近くでの振り子の振る舞いを数値計算して見てみることにする．

Θ=Pi-.001;

solution=NDSolve[{ D[θ[t], {t, 2}]== - g/L Sin[θ[t]],θ[0]==Θ,θ'[0]==0}, θ[t],{t, 0, 40}]

Plot[Evaluate[θ[t]/.solution[[1]]], {t, 0, 40}]

Θ=n πのところ（nは整数）ではあまり激しく動いていないことが読みとれる．相空間図は

g3=ParametricPlot[Evaluate[{θ[t]/.solution[[1]],D[θ[t]/.solution[[1]],t]}], {t, 0, 40}, PlotStyle->Hue[.4], PlotRange->All]

となる．上の図と一緒に描く．

Show[g1, g2,g3,AxesLabel->{"θ","((\[DifferentialD]θ)/(\[DifferentialD]t))" }, PlotRange->All]

問６ Θ=Pi のときは不安定な平衡点と呼ばれる．何が起こるか？
問７　 Θ=Pi 、  Θ=Pi+ 0.00000001 や　 Θ=Pi- 0.00000001 のときなどの数値計算をして運動の様子を調べ考察し、自分の意見（不安定な平衡点に対する認識とコンピュータによる数値計算法の良さや限界についてを述べよ．

周期

振り子の運動方程式は、

m L (((d^2)θ)/(d(t^2)))==-m g sin(θ)

と書ける．両辺に((dθ)/(dt)) をかけt について積分すると
(1/2) (((((dθ)/(dt))))^2)==-((g ∫sin(θ)dθ)/L)
t=0のとき，θ==(Subscript[θ, 0])，((dθ)/(dt))=0 とすると
((dθ)/(dt))==(Sqrt[((2 g(cos(θ)-cos((Subscript[θ, 0]))))/L)])
したがって周期τ は
τ==4 (∫_0^(Subscript[θ, 0]))(1/(Sqrt[((2 g(cos(θ)-cos((Subscript[θ, 0]))))/L)]))dθ
を計算すればよい。

g=9.8;
L=1.0;

として(Subscript[θ, 0])の値を0からπまで変化させて周期τを計算する。時間がかかるので数度に分けて計算してみよう．

list=Table[{Θ,4*NIntegrate[1/Sqrt[(2*g(Cos[θ] - Cos[Θ]))/L], {θ, 0, Θ}]}, {Θ, 0, Pi, Pi/10000}];

g4=ListPlot[list, PlotRange->All]

これではπに近づいたときの様子が詳しくわからないので，

listextra=Table[{Θ,4*NIntegrate[1/Sqrt[(2*g(Cos[θ] - Cos[Θ]))/L], {θ, 0, Θ}]}, {Θ, Pi-Pi/1000000, Pi, Pi/50000000}];

listextra1=Table[{Θ,4*NIntegrate[1/Sqrt[(2*g(Cos[θ] - Cos[Θ]))/L], {θ, 0, Θ}]}, {Θ, Pi-Pi/100000, Pi, Pi/5000000}];

g5=ListPlot[listextra, PlotRange->All]

g6=ListPlot[listextra1, PlotRange->All]

Show[g4,g5, g6]

問８　初期角度(Subscript[θ, 0])がπに近づいたとき周期τはどのように振る舞うか？