今回は、強制振動の微分方程式の解法について解説します。また、共振（共鳴）現象についても解説します。
まず、減衰振動では、失った運動エネルギーが外界から供給されないために振幅は徐々に小さくなり、最終的に振幅はゼロになります。典型的な減衰振動の運動は以下の図のようになります。
さて、減衰振動の微分方程式は、次のように表せました。
\begin{eqnarray}
m\ff{\diff^2 x}{\diff t^2} + c \ff{\diff x}{\diff t} + kx &=& 0 \\
\,
\end{eqnarray}
減衰振動の運動の様子については、実際に計算しなくとも感覚的に理解できるでしょう。では、外界からエネルギーが供給され続ける場合、どんな現象が起きるでしょうか？
今回はこの疑問に対して、強制振動の微分方程式を解くことで、考察していきます。強制振動は、機械や建物の設計、通信にも深い関係を持ちます。
強制振動の微分方程式
時間的に変動する外力の影響を受け、強制的に引き起こされる振動を強制振動と呼びます。
この外力のことを強制外力と呼びます。減衰振動では外力を加えませんでしたが、強制振動では外部から力を加える点が異なります。今回は、強制外力が三角関数で表せる場合を考えます。
この強制振動を微分方程式でモデル化すると、次のようになります。
\begin{eqnarray}
m\ff{\diff^2 x}{\diff t^2} + c \ff{\diff x}{\diff t} + kx &=& f\sin \omega t \\
\,
\end{eqnarray}
物体の質量を$m$、ばね定数を$k$、外力の角振動数を$\omega$とします。質点とばねの系に、追加でダンパーを取り付けます。（ダンパーとは、車などに付けられる衝撃を吸収する部品です）
ダンパーは、速度に比例した減衰力を発揮する特徴があります。比例定数を$c$とすると、減衰力は$\DL{c\ff{\diff x}{\diff t}}$と表せます。
外力は$f \sin\omega t$のため、この系の運動方程式を立てると、
\begin{split}
&m\ff{\diff^2 x}{\diff t^2} = -\, c \ff{\diff x}{\diff t} -\, kx + f\sin \omega t \EE
\therefore\,\,&m\ff{\diff^2 x}{\diff t^2} + c \ff{\diff x}{\diff t} + kx = f\sin \omega t
\end{split}
となります。前述の強制振動の微分方程式が導出できました。今回は、減衰のある（ダンパー付きの）強制振動について考えます。
$c=0$のとき、何が起きるのかについては後ほど考察します。今後の展開の都合上、この微分方程式を少し変形します。すなわち、両辺を$m$で割って、
\begin{eqnarray}
\ff{\diff^2 x}{\diff t^2} + \ff{c}{m} \ff{\diff x}{\diff t} + \ff{k}{m}x &=& \ff{f}{m}\sin \omega t \\
\,
\end{eqnarray}
とし、$\omega_n = \DL{\sqrt{\ff{k}{m}}}$, $2\gamma=\DL{\ff{c}{m}}$, $f_0 = \DL{\ff{f}{m}}$とおくと、
\begin{eqnarray}
\ff{\diff^2 x}{\diff t^2} + 2\gamma \ff{\diff x}{\diff t} + \omega_n^2x &=& f_0\sin \omega t \\
\,
\end{eqnarray}
と簡単にできます。※ $\omega_n$を固有角振動数と呼びます。今後の議論では、以下の微分方程式について考えていきます。
\begin{eqnarray}
\ff{\diff^2 x}{\diff t^2} + 2\gamma \ff{\diff x}{\diff t} + \omega_n^2x &=& f_0\sin \omega t \\
\,
\end{eqnarray}
今回、$c^2 \,- 4mk < 0$を仮定します。
強制振動の微分方程式の性質
強制振動の微分方程式を解いていきましょう。まずは特性方程式を使って、この微分方程式が解けるか考えましょう。特性方程式は以下のようになります。
\begin{eqnarray}
\lambda^2 + 2\gamma \lambda + \omega_n^2 &=& f_0\sin \omega t \\
\end{eqnarray}
この式を$\lambda$について解いたとして、その結果から解を計算できるでしょうか？
ちょっと手を動かせば分かりますが、この方法では微分方程式の解を求めることはできません。ではどうするのか？
少し考えてみましょう。
斉次方程式と非斉次方程式
$r(t)$を与えれた関数として、次の二階線形微分方程式を考えます。（$p, q$は定数です）
\begin{eqnarray}
\ff{\diff x^2}{\diff t^2} + p\ff{\diff x}{\diff t} + qx = r(t)
\end{eqnarray}
この微分方程式を定数係数の二階線形微分方程式と呼びます。ここで、$r(t)=0$とした、以下のような微分方程式を斉次方程式と呼びます。
\begin{eqnarray}
\ff{\diff x^2}{\diff t^2} + p\ff{\diff x}{\diff t} + qx = 0
\end{eqnarray}
一方で、最初の微分方程式を非斉次方程式と呼びます。
$r(t)=0$：斉次方程式
$r(t)\neq 0$：非斉次方程式
ということで、減衰振動の微分方程式は斉次方程式、強制振動の微分方程式は非斉次方程式に分類されます。
ここで、重要なことは、非斉次方程式では特性方程式を利用した解法が使えないということです。
ということで、強制振動方程式の解を求めるためには工夫が必要になります。
非斉次方程式を解く方法として、定数変化法やラプラス変換を利用する方法もありますが、今回は最も標準的かつ簡単な代入法を利用します。
基本解と特解
さて、非斉次方程式の一般解は、基本解と特解の線形結合として表せるという性質があります。すなわち、基本解を$x_h(t)$、特解を$x_f(t)$とすると、一般解$X(t)$は、
\begin{eqnarray}
X(t) = x_h(t) + x_f(t)
\end{eqnarray}
と表されるのです。
さて、基本解とは$r(t)=0$とした斉次方程式の解のことです。つまり、今回の強制振動の微分方程式では、減衰振動の解が基本解になります。
問題になるのは特解の$x_f(t)$です。$x_f(t)$については、地道に解く方法もありますが、偉大な先人達の計算により三角関数となることが知られています。
ショートカットのため、この結果をフルに活用していきます。具体的な計算に取り掛かる前に、一般解$X(t)$を強制振動の微分方程式に代入するとどうなるかを観察します。$X(t)$は、強制振動の微分方程式の解のため、以下の式が成立します。
\begin{eqnarray}
\ff{\diff^2}{\diff t^2}X(t) + 2\gamma\ff{\diff}{\diff t}X(t) + \omega_n^2 X(t) = f_0\sin \omega t
\end{eqnarray}
$X(t) = x_h(t)+x_f(t)$を代入すると、次のように変形でき、
\begin{split}
\ff{\diff^2}{\diff t^2}\Big( x_h(t)+x_f(t) \Big) + 2\gamma\ff{\diff}{\diff t}\Big( x_h(t)+x_f(t) \Big) \quad &\qquad \EE
+ \omega_0^2\Big( x_h(t)+x_f(t) \Big) &= f_0\sin \omega t \EE
\left( \ff{\diff^2}{\diff t^2}x_h(t) + 2\gamma\ff{\diff}{\diff t}x_h(t) + \omega_n^2x_h(t) \right) \quad &\qquad \EE
+\left( \ff{\diff^2}{\diff t^2}x_f(t) + 2\gamma\ff{\diff}{\diff t}x_f(t) + \omega_n^2x_f(t) \right) &= \EE
0+\left( m\ff{\diff^2}{\diff t^2}x_f(t) + 2\gamma\ff{\diff}{\diff t}x_f(t) + \omega_n^2x_f(t) \right) &= \EE
\therefore\, m\ff{\diff^2}{\diff t^2}x_f(t) + 2\gamma\ff{\diff}{\diff t}x_f(t) + \omega_n^2x_f(t) &= f_0\sin \omega t \EE
\end{split}
となります。$x_h(t)$は、斉次方程式の解であるため、
\begin{eqnarray}
\ff{\diff^2}{\diff t^2}x_h(t) + 2\gamma\ff{\diff}{\diff t}x_h(t) + \omega_n^2x_h(t) = 0
\end{eqnarray}
となります。さて、減衰振動の一般解は、
\begin{eqnarray}
x_h(t) = Ae^{-\gamma t}\sin\left( \sqrt{\omega_n^2\,-\gamma^2}\, t + \alpha \right)
\end{eqnarray}
でした。このことから、$x_f(t)$を単独で求め、その後で$x_h(t)$とドッキングしてやればよいことが分かります。
強制振動の一般解
$x_f(t)$が三角関数で表せることと、$x_f(t)$は単独で計算できることが分かりました。このことから、$x_f(t) = Ce^{i\omega t}$として微分方程式に代入すると、
\begin{eqnarray}
-C\omega^2e^{i\omega t} + 2i\gamma C\omega e^{i\omega t} + \omega_n^2Ce^{i\omega t} &=& f_0\sin \omega t \EE
\left( -\omega^2+2i\gamma\omega + \omega_n^2 \right)Ce^{i\omega t} &=& f_0 \sin \omega t \EE
\end{eqnarray}
オイラーの公式を利用して、左辺を三角関数に変換すると以下のようになり、
\begin{eqnarray}
C\left( \omega_n^2\,-\omega^2+2i\gamma\omega \right)(\cos \omega t + i\sin \omega t) &=& f_0 \sin \omega t \EE
\end{eqnarray}
両辺の三角関数の係数を比較すると、次のような連立方程式が得られます。
$$
\left\{
\begin{split}
&C\left( \omega_n^2\,-\omega^2+2i\gamma\omega \right)\cos \omega t = 0 \EE \EE
&iC\left( \omega_n^2\,-\omega^2+2i\gamma\omega \right)\sin \omega t = f_0 \sin \omega t
\end{split}
\right.
$$
正弦関数の係数を比較すると、次のように計算でき、係数$C$を求められます。
\begin{eqnarray}
C &=& \ff{-if_0}{\omega_n^2\,-\omega^2+2i\gamma\omega} \EE
&=& \ff{2\gamma\omega \,- (\omega_n^2\,-\omega^2)i}{(\omega_n^2\,- \omega^2)^2+(2\gamma\omega)^2}f_0
\end{eqnarray}
したがって、$x_f(t)$は、
\begin{split}
x_f(t) &= \ff{2\gamma\omega \,- (\omega_n^2\,-\omega^2)i}{(\omega_n^2\,- \omega^2)^2+(2\gamma\omega)^2}f_0\cdot e^{i\omega t} \\[10pt]
\therefore\, x_f(t) &= \ff{2\gamma\omega f_0}{(\omega_n^2\,- \omega^2)^2+(2\gamma\omega)^2}\cos\omega t \, \EE
&\quad\qquad +\, \ff{(\omega_n^2\,-\omega^2)f_0}{(\omega_n^2\,- \omega^2)^2+(2\gamma\omega)^2}\sin\omega t
\end{split}
となります。以上より、強制振動の微分方程式の一般解は、次のように表せます。
\begin{split}
X(t) &= Ae^{-\gamma t}\sin\left( \sqrt{\omega_n^2 \,- \gamma^2}\,t+ \alpha \right) \EE
&\qquad +\,\ff{2\gamma\omega}{(\omega_n^2\,- \omega^2)^2+
(2\gamma\omega)^2}f_0\cos\omega t \EE
&\qquad +\,\ff{\omega_n^2\,-\omega^2}{(\omega_n^2\,- \omega^2)^2+(2\gamma\omega)^2}f_0\sin\omega t \\
\,
\end{split}
共振の物理学
強制振動の一般解を求めることができたので、一般解の振る舞いをグラフ化して調べてみましょう。いきなりグラフ化しても意味を掴みにくいため、変数の数を減らしましょう。
非減衰強制振動
具体的には$\gamma = 0$すなわち、$c=0$の場合を考えます。（ダンパー無し）このような振動を非減衰強制振動と呼びます。非減衰強制振動では、一般解は次のようになります。
\begin{split}
X(t) &= A\sin\left(\omega_n t+ \alpha \right)+\ff{f_0}{\omega_n^2\,- \omega^2}\sin\omega t \EE
\end{split}
ここで、$\eta = \DL{\ff{\omega}{\omega_n}}$とおくと、上式は次のように整理できます。
\begin{split}
X(t)&= A\sin\left(\omega_n t+ \alpha \right)+\ff{f_0}{ \omega_n^2(1 \,-\eta^2) }\sin\omega t \EE
\end{split}
また、$\eta$（イータ）を振動数比と呼びます。以下に、$\eta = 0.7, 0.8, 0.9$それぞれの振動の様子を示します。
振動数比が$1$に近づくにつれ、整った波形になることが分かります。さらに、$\eta=0.9, 0.99, 0.999$としたときの波形を以下に示します。
グラフから、$\eta$が$1$に近づくほど振幅が急激に大きくなることが分かります。
また、基本解については係数に指数関数が含まれるため、時間経過により振幅は急激に減少していきます。
つまり、長時間のスケールで強制振動の振幅を考えると基本解の振幅は無視できるのです。よって、強制振動の振幅で興味の対象になるのは、特解$x_f(t)$の部分になります。
そこで、より詳しく強制外力により生じる振動（$x_f(t)$）の振幅の大きさと、振動数比の関係について考えてみましょう。まず、強制外力による振動の振幅は次のように変形できます。
\begin{eqnarray}
x_f &=& \ff{f_0}{ \omega_n^2(1 \,-\eta^2) }\sin\omega t\EE
&=& \ff{f}{k}\cdot\ff{1}{1 \,-\eta^2 }\sin\omega t \EE
&=& x_{st}\cdot\ff{1}{1 \,-\eta^2 }\sin\omega t \EE
\therefore \, \ff{x_f}{x_{st}} &=& \ff{1}{1 \,-\eta^2 }\sin\omega t
\end{eqnarray}
$\DL{\ff{f}{k}}$については、$f$の力を静かに加えたときのばねの変位に等しく、これを静変位と呼び、$x_{st}$とします。
このとき、$\DL{\ff{x_f}{x_{st}}}$の振幅を振幅倍率と呼びます。振幅倍率の絶対値をグラフ化すると次のようになります。
グラフより、周波数比が$1$に近づくほど振幅倍率は増加し、外力の角周波数$\omega$が固有角振動数$\omega_n$と一致したとき、強制振動の振幅が無限大になることが見て取れます。
この現象を共振（共鳴）と呼びます。そして、共振を発生させる外力の周波数を共振周波数と呼びます。
減衰強制振動
次に、$c\neq 0$の場合での振幅の様子を見てみましょう。このような振動を減衰強制振動と呼びます。（今回考えている微分方程式です）
先ほどと同様に振幅倍率を計算します。強制変位の振幅$x_f$は周波数比$\eta$と静変位$x_{st}$を使い、三角関数の合成を利用すると、次のようになります。
\begin{eqnarray}
x_f &=& \ff{f_0}{(\omega_n^2\,- \omega^2)^2+
(2\gamma\omega)^2}\Big\{ (\omega_n^2\,-\omega^2)\cos\omega t + 2\gamma\omega \sin\omega t \Big\} \EE
&=& \ff{f_0}{\omega_n^4(1\,- \eta^2)^2+\omega_n^4
(2\zeta\eta)^2}\Big\{ \omega_n^2(1\,-\eta^2)\cos\omega t + \omega_n^2\cdot 2\zeta\eta \sin\omega t \Big\} \EE
&=& \ff{f_0}{\omega_n^2}\cdot\ff{1}{(1\,- \eta^2)^2+
(2\zeta\eta)^2}\Big\{ (1\,-\eta^2)\cos\omega t + 2\zeta\eta \sin\omega t \Big\} \EE
&=& \ff{f}{k}\cdot\ff{1}{\sqrt{(1\,- \eta^2)^2+
(2\zeta\eta)^2}}\,\cos(\omega t \,- \phi) \EE
&=& x_{st}\cdot\ff{1}{\sqrt{(1\,- \eta^2)^2+
(2\zeta\eta)^2}}\,\cos(\omega t \,- \phi) \EE
\end{eqnarray}
とできて、振幅倍率の絶対値は、
\begin{eqnarray}
\left| \ff{x_f}{x_{st}} \right|
&=& \ff{1}{\sqrt{(1\,- \eta^2)^2+
(2\zeta\eta)^2}}
\end{eqnarray}
とできます。
※計算の過程で、減衰比$\zeta$を活用して式を整理しています。$\zeta = \DL{\ff{c}{2\sqrt{mk}}}$です。
振幅倍率のグラフは次のようになります。
グラフから、減衰比$\zeta$が小さくなるほど、より鋭い山になることが見て取れます。また、このグラフを共振曲線と呼びます。分母に関して平方完成を行うと、
\begin{eqnarray}
(1- \eta^2)^2+(2\zeta\eta)^2 &=& \big\{ \eta^2 -(1-2\zeta^2) \big\}^2 +4\zeta^2 (1-\zeta^2)
\end{eqnarray}
とできて、
\begin{eqnarray}
\eta^2 -(1-2\zeta^2) &=& 0
\end{eqnarray}
すなわち、
\begin{eqnarray}
\eta &=& \sqrt{1-2\zeta^2}
\end{eqnarray}
のとき、拡大振幅倍率は$\DL{\ff{1}{2\zeta\sqrt{1-\zeta^2}}}$の最大値を取ることが分かります。
※ 振幅倍率が最大になる周波数比が、$\eta = 1$では無いことに注意してください。
したがって、共振周波数$\omega_r$は、次のようになります。
\begin{eqnarray}
\omega_r &=& \omega_n \sqrt{1-2\zeta^2} \\
\,
\end{eqnarray}
$Q$値
共振現象は特定の周波数を強く増幅する現象とも言えます。共振が起きた時、どの程度の範囲の周波数を強く増幅するかの指標として、Q値と呼ばれる値があります
Q値が小さいほど、より狭い範囲の周波数がより大きく増幅されることを意味します。つまり、Q値が小さいほどより鋭い共振曲線になるということです。
まず、振幅倍率が$1$より大きくなるためには、計算から
\begin{split}
&1 < \ff{1}{2\zeta\sqrt{1-\zeta^2}} \EE
&\therefore \zeta < \ff{1}{\sqrt{2}}
\end{split}
$\zeta$を$\DL{\ff{1}{\sqrt{2}}}$より小さくする必要があることが分かります。
逆に、共振を起こしたくない場合、減衰比$\zeta$が$\DL{\ff{1}{\sqrt{2}}}$より大きくなるように設計すれば良いことが分かります。
共振の工学との関係について
理論的な話は一旦脇に置いて、工学的な話に移りましょう。共振現象は特定の周波数を強く増幅する現象とも言えます。
共振現象には、厄介な側面と有用な側面があります。通信の分野では特定の周波数を取り出すため、共振現象が利用されています。
現代のIT社会の発展に共振現象は大きく貢献しているのです。一方で、共振現象が厄介な事態を引き起こすのは、建造物や回転機械のような建築や機械に関わる分野です。
これまで見てきたように、外力の振動数が建物や機械の固有振動数に一致すると共振現象が発生します。運悪く共振現象が発生すると、建物や機械が非常に大きな振幅で揺さぶられることになります。
最悪の場合、建物の倒壊や機械の破壊につながります。そのため、設計者は共振が起きないよう事前に対策を施す必要があります。
東日本大震災時に共振により、特定の高層ビルが周りの建物より大きく揺さぶられた事例や、イギリスのミレニアムブリッジの固有振動数が通行者の歩く周期と一致し共振が発生した事例などがあります。
一方、機械で共振現象が問題になるのは、回転機械です。機械に回転機構が良く使われるため、共振が起きないように注意する必要があります。回転機械が共振により破壊された例として、発電所のタービンが共振により破壊に至った事例があります。
回転機械にて共振が発生する回転速度（角速度）を危険速度と呼び、設計にあたっては危険速度を避けるように設計するか、迅速に危険速度を通過するように心がける必要があります。
危険速度については、こちらの記事で解説しています。
