強制振動

DSolveを用いてとにかく力任せに解いてみる(brute-force method)

まず外力が正弦関数的に変化する場合の強制振動の例から論じることにしよう．減衰
振動をする系にこのような外力 F =(Subscript[F , 0])cos(Ω t)がかかったとすると、運動方程式は次のように書ける．

x(t) k+2b x'(t)+m (x')'(t)==(Subscript[F , 0])cos(Ω t)

または，

x(t) (Subscript[ω, 0]^2)+2 β x'(t)+(x')'(t)==A cos(Ω t) 

ここで， A = (Subscript[F, 0])/m で，その他の記号の定義はこれまでどうりである．

ここでは力任せに解を求めてみよう．

次の方程式を解く．初期条件はx[0]==(Subscript[X, 0]),x'[0]==(Subscript[V, 0])とする．

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

解を定義しておく．その前に，解を少し見やすくしておく．

solutionNew=FullSimplify[solution]

x[t_]=x[t]/.solutionNew[[1]];

相当に複雑な解であるが（教科書は振幅、位相等を導入して見やすく書いてある．）、とにかく、変位、速度を適当な初期条件でプロットしてみよう．

変位は

Plot[Evaluate[x[t]/.{(Subscript[ω, 0])->1,  β->.2, Ω->2.5, (Subscript[X, 0])->0, (Subscript[V, 0])->1, A->1}], {t, 0, 50}]

となり、速度は、

Plot[Evaluate[D[x[t],t ]/.{(Subscript[ω, 0])->1,  β->.2, Ω->2.5, (Subscript[X, 0])->0, (Subscript[V, 0])->1, A->1}], {t, 0, 50}]

となる．時間tが小さいところは複雑に振動するが、時間がたつと余関数が指数関数的に小さくなり（教科書参照）、特解のみが生き残る．

問１　上の図で時間が十分たったときの角振動数はいくらか．

問２　減衰パラメターbetaをもっと大きくすると系はどのように振る舞うか．

hint:以下のようにベターの値を変えて変位を時間の関数としてプロットしてみなさい．アウトプットはアニメーションで見て下さい．

y[t_,β_]=x[t]/.{(Subscript[ω, 0])->1, Ω->2.5, (Subscript[X, 0])->0, (Subscript[V, 0])->1, A->1}

Manipulate[Plot[y[t,β],{t,0,50},PlotRange->{-1,1.2},PlotLabel->{"β=",β},PlotStyle->Hue[β]],{β,0.1001,3,0.1}]

位相空間図

ここでParametricPlotを使って位相空間図を描いてみよう．

ParametricPlot[Evaluate[Chop[{ x[t], D[x[t],t ]}/.{(Subscript[ω, 0])->1,  β->.2, Ω->2.5, (Subscript[X, 0])->0, (Subscript[V, 0])->1, A->1}]], {t, 0, 50}, AxesLabel->{変位,速度}]

問３　上の位相空間図を解釈せよ．（時間の経過に運動がどのように変わるか．）また、ベターの値を変えると位相空間図はどのように変わるか．

Hint: ベータを変えながら、速度と位置の相空間図を描いてみる．

z[t_,β_]={x[t],(∂_t)x[t]}/.{(Subscript[ω, 0])->1,Ω->2.5`,(Subscript[X, 0])->0, (Subscript[V, 0])->1,A->1};

Manipulate[ParametricPlot[z[t,β],{t,0,50},PlotLabel->{"(β)" β},PlotRange->{{-1,1.2`},{-1,1.2`}},AspectRatio->1,PlotStyle->Hue[β]],{β,0.05`,3,0.1`}]

複素数を用い特解を求める

x(t) (Subscript[ω, 0]^2)+2 β x'(t)+(x')'(t)==A Re[Exp[Ω t)] 

であることに注意して，

X[t_]=B Exp[I Ω t]

と仮定して，与えられた方程式に代入し，B を求める．

Solve[X''[t] +2β X'[t]+(Subscript[ω, 0])^2 X[t]==A Exp[I Ω t], B]

B を結果のように定義する．

B=B/.%[[1]]

B の絶対値 |B| を求めるには、

ComplexExpand[Abs[B]]

Aが正であるとして，

Simplify[%, A>0]

また，偏角をαとすると，αは

α = ComplexExpand[Arg[B],TargetFunctions->{Re,Im}]

または，

FullSimplify[Tan[α]]

以上より，X''[t] +2β X'[t]+(Subscript[ω, 0])^2 X[t]==A Exp[I ω t]　の特解は(Subscript[X, p])は

(Subscript[X, p])[t_]:= (A/(Sqrt[4 (β^2) (Ω^2)+(Ω^4)-2 (Ω^2) (Subscript[ω, 0]^2)+(Subscript[ω, 0]^4)])) Cos[Ω t +α]

で，

Tan[α]=((2 β Ω)/((Ω^2)-(Subscript[ω, 0]^2)))
である．

共鳴

とにかく力任せに解いてみる(brute-force method)の続き

上の結果を使って議論しよう．共鳴とは自然振動数(Subscript[ω , 0])に強制振動の振動数Ωが近づいたとき振幅が大きくなり、位相が大きく変わる現象を言う．（教科書を参照）．

まず、Ωを変えながら位置を時間の関数として３次元的に描いてみよう．

Plot3D[Evaluate[x[t]/.{(Subscript[ω, 0])->1,  β->.1, (Subscript[X, 0])->0, (Subscript[V, 0])->1, A->1}], {t, 0, 50}, {Ω, 0.1, 2},  PlotRange->All]

上の図は少しわかりにくのでrotate して真正面から見てみよう．

Plot3D[Evaluate[x[t]/.{(Subscript[ω, 0])->1,  β->.1,  (Subscript[X, 0])->0, (Subscript[V, 0])->1, A->1}], {t, 0, 50}, {Ω, 0.1, 2},  PlotRange->All,ViewPoint->{3.382, 0.019, 0.109}]

Ω の値が１に近づくと振幅が大きくなる様子が読みとれる．βの値をもう少し小さくすると，

Plot3D[Evaluate[x[t]/.{(Subscript[ω, 0])->1,  β->.005,  (Subscript[X, 0])->0, (Subscript[V, 0])->1, A->1}], {t, 0, 50}, {Ω, 0.1, 2},  PlotRange->All,ViewPoint->{3.382, 0.019, 0.109}]

さて、振幅の大きさを見つけよう．時間がたつと系は外力と同じ周期で変化する．時間が十分たったとしてFindMaximumを使って振幅の最大値を探そう．

FindMaximum[Chop[Evaluate[x[t]/.{(Subscript[ω, 0])->1,  β->.2, Ω->2.5,(Subscript[X, 0])->0, (Subscript[V, 0])->1, A->1}]], {t, 100}]

この結果はt =100の近くで最大値を探すと0.187112であるとアウトプットは示している．

OMEGA の値を変えつつ最大値の数列を作ってみよう．（この計算は少し時間がかかる．また，計算が収束しない場合もあるが無視することにする．）

p[t_, Ω_]=x[t]/.{(Subscript[ω, 0])->1,  β->.2,(Subscript[X, 0])->0, (Subscript[V, 0])->1, A->1};

listup= Table[{Ω, First[FindMaximum[Abs[p[t,Ω]], {t,100.}]]}, {Ω, .1, 2, .06}]//Timing

ListPlot[listup[[2]],PlotRange->{{0,2},{0,3}},Joined->True,PlotStyle->Hue[1], PlotRange->All]

Ωが(Subscript[ω, 0])の値（＝１）に近づいたとき振幅が最大値をとると推測できる．

問４　Ωが(Subscript[ω, 0])の値（＝１）に近づいたとき振幅を大きくするためにはどのような減衰パラメターを選べばよいか．

つぎにΩの変化に従いどのように変化するかみよう．

listupPhase= Table[{Ω,t/.FindMinimum[ Chop[x[t]/.{(Subscript[ω, 0])->1,  β->.2,(Subscript[X, 0])->0, (Subscript[V, 0])->1, A->1}], {t, 101}][[2]]}, {Ω, .1, 2, .06}]

これをプロットすると

ListPlot[listupPhase]

上の図は極小値を与える時間t をOMEGAの関数としてプロットしたものである．従って位相はOMEGAの関数として次のように描ける．

ListPlot[Table[{listupPhase[[n,1]], Mod[listupPhase[[n,1]] listupPhase[[n,2]], Pi]}, {n, Length[listupPhase]}]]

この図は共鳴の起こるところで(OMEGA ~1)で位相が大きく変わる様子が読みとれる．

注：ここで位相が2Piでもとに戻るのでMod を用いた．例としては

Mod[3.6, 2]

Mod[9.2, 2]

のようにMod[m,n]はmをnで割った残りを示す．

”複素数を用い特解を求める”の続き

興味のあるのは特解のみであるので

(Subscript[X, p])[t_]:= (A/(Sqrt[4 (β^2) (Ω^2)+(Ω^4)-2 (Ω^2) (Subscript[ω, 0]^2)+(Subscript[ω, 0]^4)])) Cos[Ω t +α]

と定義する．

(Subscript[X, p])[t]/.{(Subscript[ω, 0])->1,  β->.1, (Subscript[X, 0])->0, (Subscript[V, 0])->1, A->1}

Plot3D[Evaluate[(Subscript[X, p])[t]/.{(Subscript[ω, 0])->1,  β->.1, (Subscript[X, 0])->0, (Subscript[V, 0])->1, A->1}], {t, 0, 50}, {Ω, 0.1, 2},  PlotRange->All]

振幅Ｂのみをプロットすると，

Plot3D[Evaluate[(Subscript[X, p])[t]/Cos[Ω t +α]/.{(Subscript[ω, 0])->1,   (Subscript[X, 0])->0, (Subscript[V, 0])->1, A->1}], {Ω, 0, 2}, {β, 0.01, .3},  PlotRange->All]

振幅Abs[Ｂ]を最大とする振動数 (Subscript[Ω, R]) を振幅共鳴振動数 (amplitude resonancefrequency) という．即ち、

Solve[ D[(Subscript[X, p])[t]/Cos[Ω t +α],Ω]==0, Ω]

より,正の解のみひろうと

(Subscript[Ω, R]) = Ω/.%[[3]]

となる．．これからわかるように β が大きくなれば共鳴振動数 (Subscript[Ω, R]) は小さくなり，もっと大きくなると(Subscript[Ω, R]) は虚数となり 振幅 は Ω が大きくなるに従い減少する．
β が (Subscript[ω, 0])に比べて小さいとき、共鳴振動数 (Subscript[Ω, R]) はほとんど減少しないが、Ω ∼(Subscript[ω, 0])のとき当然振幅 は大きな変化をみせる．

位相のみをプロットすると

Plot3D[Evaluate[α/.{(Subscript[ω, 0])->1,   (Subscript[X, 0])->0, (Subscript[V, 0])->1, A->1}], {Ω, 0, 2}, {β, 0.05, 1},  PlotRange->All]

β が小さくなるに従って位相αは自然振動数 (Subscript[ω, 0]) の近傍にあるときより急激に変化する．
β = 0 ではその変化が π になることに注意する．

Q値

強制振動では、よく Q 値というものが引用される．これは共鳴の鋭さをあらわすパラメターで
Q≡((Subscript[Ω, R])/(2 β))
で定義される． Q 値が大きいほど共鳴は鋭い．ここでは、緩やかな減衰振動の場合、すなわち減衰パラメター β が自然振動数 (Subscript[ω, 0]) 比べて小さいとき、
Q≈((Subscript[ω, 0])/(Δω))
であることを示す．ここでΔω は半値幅 (FWHM: Full Width at Half Maximum) と呼ばれ、振幅がその最大値の丁度(1/(Sqrt[2]))になった２点をむすんだ振動数の幅に対応する．すなわち，

Solve[(Subscript[X, p])[t]/Cos[Ω t +α]==Evaluate[1/Sqrt[2] (Subscript[X, p])[t]/Cos[Ω t +α]/.Ω->(Sqrt[-2 (β^2)+(Subscript[ω, 0]^2)])],Ω]

を得る．負の解は物理的に意味をなさないので第２番目と第４番目の解のみひろい、 x1と x2 とする．即ち、

x1=Ω/.%[[2]]

x2=Ω/.%%[[4]]

と定義し、半値幅を Δω とよべば、それは、（ b について１次のオーダーのみ計算す
る）

X2-X1

Δω=Normal[ Series[x2-x1,{β,0,1}]]

となる．また，(Subscript[Ω, R])は

Simplify[Normal[Series[(Subscript[Ω, R]),{β,0,1}]],(Subscript[ω, 0])>0]

であり，従って
Q≡((Subscript[Ω, R])/(2 β))≈　((Subscript[ω, 0])/(Δω))
であることが証明された．

問５　以上の結果を踏まえ，強制振動とは何か、分かり易く説明せよ．

問６　共鳴とは何か、分かり易く説明せよ．

問７　共強制振動する振動子の力学的エネルギーはどのような時間変化をするか．計算してプロットしてみよ．