減衰振動

教科書に説明してあるように、抗力が速度に比例するとし、その比例定数をk とすると、運動方程式は

m (((d^2)x)/(d(t^2)))==-k x-b ((dx)/(dt))

と書ける．

β==(b/m);
(Subscript[ω, 0])==(Sqrt[(k/m)]);

とし、初期条件

x[0]==(Subscript[X, 0]);
x'[0]==(Subscript[V, 0]);

として、運動方程式を解いてみよう．

DSolve[{((x'))'[t]+2 β x'[t]+(Subscript[ω, 0]^2) x[t]==0,x[0]==(Subscript[X, 0]),x'[0]==(Subscript[V, 0])},x[t],t]

もう少し見やすくしておこう．

FullSimplify[%]//Timing

これが上の初期条件を満たす解である．　x[t]をこの解であるとすると定義する．

x[t_] =x[t]/.%[[2,1]]

問１　ここで双曲線関数について述べておく．

TrigToExp[Sinh[x]]

TrigToExp[Cosh[x]]

TrigToExp[Tanh[x]]

である．まとめると，、

Sinh[x]==-((e^(-x))/2)+((e^x)/2);
Cos[x]==((e^(-x))/2)+((e^x)/2);
Tanh[x]==((-(e^(-x))+(e^x))/((e^(-x))+(e^x)))

さて、Mathematicaを使わず 以下を証明せよ．
(a) Sinh[i x]==i Sin[x]
(b) ((Cosh[x])^2)-((Sinh[x])^2)==1

減衰パラメターがゼロの場合

β =0であるから、

β=0;

として、

x[t]

若干見にくいので、

FullSimplify[%]

としてみると、分かり易い．初期位置を変化させながら、位置を時間の関数としてプロットしてみよう．

Plot[Evaluate[Table[{x[t]}/.{(Subscript[ω, 0])->1,(Subscript[V, 0])->1},{(Subscript[X, 0]),0,10,.5}]],{t,0,20 π},PlotStyle->Table[Hue[(Subscript[X, 0])/10],{(Subscript[X, 0]),0,10,.5}],AxesLabel->{時間,位置}]

または，

ParametricPlot[Evaluate[Table[{t,x[t]}/.{(Subscript[ω, 0])->1,(Subscript[V, 0])->1},{(Subscript[X, 0]),0,10,.5}]],{t,0,20 π},PlotStyle->Table[Hue[(Subscript[X, 0])/10],{(Subscript[X, 0]),0,10,.5}],AxesLabel->{時間,位置}]

相空間図を描いてみよう．

ParametricPlot[Evaluate[Table[{x[t],(∂_t)x[t]}/.{(Subscript[ω, 0])->1,(Subscript[V, 0])->1},{(Subscript[X, 0]),0,10,.5}]],{t,0,2 π},PlotStyle->Table[Hue[(Subscript[X, 0])/10],{(Subscript[X, 0]),0,10,.5}],AxesLabel->{位置,速度}]

問２　この場合、系の力学的エネルギーはどのような時間変化をするか．実際に運動エネルギーと位置エネルギーの和を計算してみよ．

減衰パラメターが小さい場合(減衰振動)　　　β<(Subscript[ω, 0])

以下、減衰パラメターは、β　は

β<(Subscript[ω, 0])

と仮定する．

β =0と上で定義していたので、

β=.

として．その定義を消しておく．位置は

x[t]/.{(Sqrt[(β^2)-(Subscript[ω, 0]^2)])->I(Sqrt[(Subscript[ω, 0]^2)-(β^2)]),1/(Sqrt[(β^2)-(Subscript[ω, 0]^2)])->1/I(Sqrt[(Subscript[ω, 0]^2)-(β^2)]) }

Simplify[%]

である．

減衰パラメターβ　を変化させながら、位置を時間の関数としてプロットしてみよう．

Plot[Evaluate[Table[{x[t]}/.{ (Subscript[ω, 0])->1,(Subscript[V, 0])->1,(Subscript[X, 0])->0}, {β,0, 1-0.01, .1}]],{t, 0, 10Pi}, PlotStyle->Table[Hue[α],{α, 0, 1, .1}], AxesLabel->{ 時間,位置}, PlotRange->All]

減衰パラメターを変化させながら、相空間図を描いてみよう．

graph1=ParametricPlot[Evaluate[Table[{x[t],D[x[t],t]}/.{ (Subscript[ω, 0])->1,(Subscript[V, 0])->1,(Subscript[X, 0])->0}, {β,0, 1-0.01, .1}]],{t, 0, 10Pi}, PlotStyle->Table[Hue[α],{α, 0, 1, .1}], AxesLabel->{ 位置,速度}, PlotRange->All]

問３　上の図は何を表しているのか．物理的な解釈を行え．

アニメーション

速度v[t]を以下のように定義しておく．

v[t_]=(∂_t)x[t]

または，

FullSimplify[%]

さて、上でつかったgraph1を用いて、減衰パラメターを変化させながら、物体の位置の動きがわかるアニメーションをつくる．

Manipulate[Show[graph1,Chop[Table[Graphics[{Hue[.7],Disk[Evaluate[{x[t],v[t]}/.{(Subscript[ω, 0])->1,(Subscript[V, 0])->1,(Subscript[X, 0])->0}],.03]}],{β,.05,1-.1,.1}]]],{t, 0, 60,.2}]

問４　上のアニメーションから、何がわかるか、物理的に解釈して分かりやすく説明せよ．減衰パラメターβ　がゼロに近いとき，また，(Subscript[ω, 0])に近いときどのように軌跡は振る舞うか．

減衰パラメターが大きい場合(過減衰)　　　β>(Subscript[ω, 0])

以下、減衰パラメターは、β　は

β>(Subscript[ω, 0])

位置は

x[t]

である．

減衰パラメターβ　を変化させながら、位置を時間の関数としてプロットしてみよう．

Plot[Evaluate[Table[{x[t]}/.{ (Subscript[ω, 0])->1,(Subscript[V, 0])->1,(Subscript[X, 0])->0}, {β, 1+0.1,3, .2}]],{t, 0, 10Pi}, PlotStyle->Table[Hue[α],{α, 0, 1, .1}], AxesLabel->{ 時間,位置}, PlotRange->All]

減衰パラメターを変化させながら、相空間図を描いてみよう．

graph2=ParametricPlot[Evaluate[Table[{x[t],D[x[t],t]}/.{ (Subscript[ω, 0])->1,(Subscript[V, 0])->1,(Subscript[X, 0])->0}, {β, 1+0.1,3, .2}]],{t, 0, 10Pi}, PlotStyle->Table[Hue[α],{α, 0, 1, .1}], AxesLabel->{ 位置,速度}, PlotRange->All]

アニメーション

Manipulate[Show[graph2,Chop[Table[Graphics[{Hue[.7],Disk[Evaluate[{x[t],v[t]}/.{(Subscript[ω, 0])->1,(Subscript[V, 0])->1,(Subscript[X, 0])->0}],.01]}],{β,1.01,3,.1}]]],{t,0,10,.5}]

β < (Subscript[ω, 0]) の時は軌跡は螺旋を描き原点に収束する． β ≥ (Subscript[ω, 0]) 0 の場合は螺旋は描かず直接原点に収束する．この収束した点を 固定点 (fixed point) またはポイントアトラクター (point
attractor) と呼ぶ．このケースでは ポイントアトラクターは一つであるが、カオスを起こしうるの場合では ポイントアトラクター が有数個であったり無い事がある．

問５上のアニメーションを解釈して分かりやすく説明せよ．βの値が何の時一番早く固定点に収束するか．

graph2a=ParametricPlot[Evaluate[Table[{x[t],D[x[t],t]}/.{ (Subscript[ω, 0])->1,(Subscript[V, 0])->Cos[θ],β->1.5,(Subscript[X, 0])->Sin[θ]}, {θ,0,2Pi, Pi/30}]],{t, 0, 10Pi}, PlotStyle->Table[Hue[α],{α, 0, 1, .1}], AxesLabel->{ 位置,速度}, PlotRange->All]

Manipulate[Show[graph2a,Chop[Table[Graphics[{Hue[.7],Disk[Evaluate[{x[t],v[t]}/.{ (Subscript[ω, 0])->1,(Subscript[V, 0])->Cos[θ],β->1.5,(Subscript[X, 0])->Sin[θ]}],.01]}],{θ,0,2Pi, Pi/30}]]],{t,0,10,.1}]

減衰パラメターが自然角振動数に等しい場合(臨界減衰)　　　β=(Subscript[ω, 0])

以下、減衰パラメターは、β　は

β=(Subscript[ω, 0])

位置は

x[t_]=Limit[x[t],β->(Subscript[ω, 0])]

である．

初期位置を変化させながら、位置を時間の関数としてプロットしてみよう．

Plot[Evaluate[Table[{x[t]}/.{(Subscript[ω, 0])->1,(Subscript[V, 0])->1},{(Subscript[X, 0]),0,10,.5}]],{t,0,3 π},PlotRange->All,PlotStyle->Table[Hue[(Subscript[X, 0])/10],{(Subscript[X, 0]),0,10,.5}],AxesLabel->{時間,位置}]

初期位置を変化させながら、相空間図を描いてみよう．

graph3=ParametricPlot[Evaluate[Table[{x[t],D[x[t],t]}/.{ (Subscript[ω, 0])->1,(Subscript[V, 0])->1}, {(Subscript[X, 0]),0,3, .2}]],{t, 0, 10Pi}, PlotStyle->Table[Hue[α],{α, 0, 1, .1}], AxesLabel->{ 位置,速度}, PlotRange->All]

アニメーション

v[t_]=(∂_t)x[t]

Simplify[%]

x[t]

Manipulate[Show[graph3,Chop[Table[Graphics[{Hue[.7],Disk[Evaluate[{x[t],v[t]}/.{(Subscript[ω, 0])->1,(Subscript[V, 0])->1}],.04]}],{(Subscript[X, 0]),0,3, .2}]]],{t, 0,7,.2}]

graph3a=ParametricPlot[Evaluate[Table[{x[t],D[x[t],t]}/.{ (Subscript[ω, 0])->1,(Subscript[V, 0])->Cos[θ],β->1,(Subscript[X, 0])->Sin[θ]}, {θ,0,2Pi, Pi/30}]],{t, 0, 10Pi}, PlotStyle->Table[Hue[α],{α, 0, 1, .1}], AxesLabel->{ 位置,速度}, PlotRange->All]

Manipulate[Show[graph3a,Chop[Table[Graphics[{Hue[.7],Disk[Evaluate[{x[t],v[t]}/.{ (Subscript[ω, 0])->1,(Subscript[V, 0])->Cos[θ],β->1,(Subscript[X, 0])->Sin[θ]}],.01]}],{θ,0,2Pi, Pi/30}]]],{t,0,10,.1}]

問６上のアニメーションを解釈して臨界減衰について分かりやすく説明せよ．

問７　過減衰においては、臨界減衰の時に比べて減衰パラメターの値が大きいのに固定点（原点）にたどり着くのに時間が多くかかるのはなぜか．また、どうして時間がたつと、直線に漸近して原点に収束するが、なぜか（教科書を参照）．

問８上のアニメ-ションはわかりやすくするためにβ,(Subscript[ω, 0])の二つのパラメターの大きさの違いによって場合分けをしたが、適当にパラメターを選び、全ての場合を一つの図にして、上と同じようにアニメーションを作るにはどうインプットをすれば良いか．