振り子の運動を題材に、微分方程式の解法を考えていきます。
今回は、四通りの方法で振り子の微分方程式の解法を解説します。早速、次のような振り子を考えましょう。
糸の長さを$l$、小球の質量を$m$、糸が鉛直方向と成す角を$\theta$とします。また、重力加速度を$g$とし、糸の重さは無視できるものとします。
このとき、振り子の周期$T$は、
\begin{eqnarray}
T &=& 2\pi \sqrt{\ff{l}{g}} \\
\,
\end{eqnarray}
と表せると、高校時代に天下り的に教えられましたが、その導出過程については追及しませんでした。今回は、振り子の微分方程式を解くことで、この謎について考えていきます。
※微分の表記法は以下の記事で解説しているので、参考にしてください。
振り子の運動方程式の導出
振り子の運動に関する問題なので、運動方程式を考えなければなりません。まずは座標軸を定めましょう。
今回は、周方向に沿った軸を設定し、これを $\theta$（シータ）方向とします。
次に、加速度を$a$として、$\theta$方向での運動方程式を考えます。
\begin{eqnarray}
ma &=& -mg\sin\theta(t) \\
\end{eqnarray}
最初に加速度について考えます。
今、$\diff t$の間に$\diff \theta(t)$だけ角度が変化したとします。すると、この間に小球は、$l\diff \theta(t)$だけ移動します。
したがって、小球の加速度$a$は、$\displaystyle l\ff{\diff^2 \theta(t)}{\diff t^2}$となります。
以上より振り子の運動方程式は、
\begin{eqnarray}
m l\ff{\diff^2 \theta(t)}{\diff t^2} &=& -mg\sin\theta(t) \EE
\ff{\diff^2 \theta(t)}{\diff t^2} &=& -\ff{g}{l}\sin\theta(t) \tag{1}
\end{eqnarray}
と求められます。
振り子の運動方程式の標準的な解法
振り子の運動方程式を解くことを考えましょう。具体的には、$\theta$の時間変化を表す関数、$\theta(t)$を求めることがゴールになります。
さて、式(1)には三角関数が含まれていますが、このような微分方程式は非線形微分方程式に分類されます。
一般に、非線形微分方程式を解くことはかなり難しいため、より簡単な微分方程式にします。
微分方程式を簡単にするため、角度$\theta$は微小であるという条件を課すことにします。
すると、$\sin\theta \NEQ \theta$と近似できるので、式(1)を次のように簡単にできます。（簡略化のため、$\theta(t)$を$\theta$と表しています。）
\begin{eqnarray}
\ff{\diff^2 \theta}{\diff t^2} &=& -\ff{g}{l}\theta \tag{2}
\end{eqnarray}
式(2)のような微分方程式は線形微分方程式に分類されます。特に、最高階数が二階微分なので、二階線形微分方程式と呼ばれます。
線形微分方程式は比較的簡単に解けるため、式(2)について解いて行きます。
指数関数の利用
さて、線形微分方程式の場合、解が指数関数として表せるという性質があります。
線形微分方程式の解は、指数関数により表せる。
この性質を利用しましょう。
$\theta (t) = Ce^{\lambda t}$と仮定し、式(2)に代入します。（$C, \lambda$は複素数の範囲で動きます。）
すると、$\lambda$は、
\begin{eqnarray}
\lambda^2 Ce^{\lambda t} &=& -\ff{g}{l}Ce^{\lambda t} \EE
\lambda^2 &=& -\ff{g}{l} \tag{3} \EE
\lambda &=& \pm i\sqrt{\ff{g}{l}}
\end{eqnarray}
と求められます。
$\lambda$の答えとして、虚数解が二つ出てきました。
線形結合の利用
線形微分方程式のもう一つの性質を使います。
その性質とは、微分方程式の解は線形結合により表せるというものです。
線形微分方程式の解は、線形結合により表せる。
線形結合をざっくり説明すると、足し算で表せるということです。
今回の場合、$\theta (t)$はこのように表せます。（$C, C’$は積分定数です。）
\begin{eqnarray}
\theta (t) &=& Ce^{ i\sqrt{\ff{g}{l}}t } + C’e^{ -i\sqrt{\ff{g}{l}}t } \tag{4} \EE
\end{eqnarray}
これを式(2)に代入すると、両辺が等しくなることが確かめられます。
オイラーの公式の利用
指数関数の肩に虚数が入っていてますが、このままでは理解が難しいので、オイラーの公式というものを利用します。
オイラーの公式は虚数を含む指数関数と三角関数を繋げる公式で、線形微分方程式を解く上で最重要の公式の一つです。オイラーの公式を利用して、式(4)を変形すると、
\begin{split}
\theta (t) &= Ce^{ i\sqrt{\ff{g}{l}} } + C’e^{ -i\sqrt{\ff{g}{l}} } \EE
&= C_1 \cos\sqrt{\ff{g}{l}}t + iC_2\sin\sqrt{\ff{g}{l}}t \\
&\qquad+C_3 \cos\sqrt{\ff{g}{l}}t \,-\, iC_4\sin\sqrt{\ff{g}{l}}t \\
&= (C_1+C_3)\cos\sqrt{\ff{g}{l}}t + i(C_2-C_4)\sin\sqrt{\ff{g}{l}}t
\end{split}
となります。
さて、物理的に考えて、$\theta$が虚数になることはあり得ないので、$C_2 \,-\, C_4 = 0$とし、虚部を消去します。
また、$C=C_1+C_3$とすると、
\begin{eqnarray}
\theta (t) = C\cos\sqrt{\ff{g}{l}}t
\end{eqnarray}
となり、$\theta$の時間変化が求められました。
境界条件の検討
積分定数を決定します。
振り子の初期の角度を$\theta_0$とします。
つまり、$t=0$のとき、$\theta(0) = \theta_0$なので、$C=\theta_0$となります。
従って$\theta(t)$は、以下のよう表せます。
\begin{eqnarray}
\theta (t) = \theta_0\cos\sqrt{\ff{g}{l}}t \tag{5}
\end{eqnarray}
周期$T$の計算
振り子の周期$T$を求めましょう。
振り子が反対側の折り返し地点に来た時、角度は$-\theta_0$となり、ここまで$\displaystyle \ff{T}{2}$で到達します。
このことを利用し、周期を求めます。
式(5)に適用すると、
\begin{eqnarray}
\theta \left( \ff{T}{2} \right) &=& \theta_0\cos \sqrt{\ff{g}{l}}\ff{T}{2} \EE
-\theta_0 &=& \theta_0\cos \sqrt{\ff{g}{l}}\ff{T}{2} \EE
-1 &=& \cos \sqrt{\ff{g}{l}}\ff{T}{2}
\end{eqnarray}
となります。
さて、左辺が$-1$であるため、右辺は$\cos \pi$となります。
したがって、
\begin{eqnarray}
\pi &=& \sqrt{\ff{g}{l}}\ff{T}{2} \EE
T &=& 2\pi \sqrt{\ff{l}{g}}
\end{eqnarray}
となります。
周期$T$が求められました。先ほどの振り子の周期の公式と一致することが分かります。
特性方程式を用いた解法
微分方程式の標準的な解法を説明しましたが、なかなか面倒でした。
式(3)に注目してください。
微分方程式を解く過程で、二次方程式がでてきました。これは偶然でしょうか？
今回の場合に限った話かどうかを確かめるため、次のような一般的な二階微分方程式を考えましょう。
\begin{eqnarray}
\ff{\diff^2 x}{\diff t^2} + \alpha\ff{\diff x}{\diff t} + \beta x = 0 \tag{6}
\end{eqnarray}
ただし、$\alpha, \beta$は実数定数とします。
線形微分方程式の性質①を使い、この微分方程式の解を$x(t) = Ce^{\lambda t}$と置きます。
式(6)に代入すると、
\begin{eqnarray}
\lambda^2 Ce^{\lambda t} + \alpha\lambda Ce^{\lambda t} + \beta Ce^{\lambda t} &=& 0 \EE
(\lambda^2 + \alpha\lambda + \beta)Ce^{\lambda t} &=& 0
\end{eqnarray}
となります。
$Ce^{\lambda t} \neq 0$なので、$\lambda^2 + \alpha\lambda + \beta = 0$となるはずです。
一般の二階線形微分方程式の場合でも、二次方程式が現れました。
この二次方程式を$\lambda$について解けば、微分方程式が解けそうです。
そこで、この二次方程式を特性方程式と呼ぶことにしましょう。
特性方程式は、指数関数の指数に注目した方程式です。
微分方程式に指数関数を代入して、式を整理せずとも、特性方程式を立式して計算してやれば指数を求めることができます。
その後は、先ほどと同様に微分方程式を解いて行きます。
特性方程式による解法は式(6)の右辺が$0$となる場合のみ使えます。右辺が$0$とならない場合には違うアプローチが必要になります。
演算子法による解法
特性方程式による解法に似た解法として、演算子法を紹介します。
演算子$D$というものを考えます。
ここで、$D$を次のように定義します。
\begin{eqnarray}
D \equiv \ff{\diff }{\diff t}\, , \qquad D^2 \equiv \ff{\diff^2 }{\diff t^2}
\end{eqnarray}
演算子法を式(6)に適用すると、
\begin{eqnarray}
D^2 x + \alpha Dx + \beta x = 0
\end{eqnarray}
となり、$x$についてまとめると、
\begin{eqnarray}
(D^2 + \alpha D + \beta) x = 0
\end{eqnarray}
となります。
さて、振り子の微分方程式に演算子法を適用すると、
\begin{eqnarray}
\left(D^2 + \ff{g}{l}\right) \theta = 0
\end{eqnarray}
となります。
この式は恒等式となるので、$\displaystyle \left(D^2 + \ff{g}{l}\right) = 0 $となる必要があります。
$D$について解くと、
\begin{eqnarray}
D = \pm i \sqrt{\ff{g}{l}}
\end{eqnarray}
となります。
両辺に$\theta$を掛けて、演算子を微分の形に直すと、
\begin{split}
&D\theta = \pm i \sqrt{\ff{g}{l}} \theta \EE
&\ff{\diff \theta}{\diff t} = \pm i \sqrt{\ff{g}{l}} \theta \EE
\therefore\, &\ff{\diff \theta}{\pm i \sqrt{\ff{g}{l}} \theta } = \diff t
\end{split}
となります。
両辺を積分して$\theta$について計算すると、
\begin{split}
&\int \ff{\diff \theta}{\pm i \sqrt{\ff{g}{l}} \theta } = \int \diff t \EE
&\pm i \sqrt{\ff{l}{g}}\ln \theta + C = t \EE
\therefore\, &\ln \theta = \pm i\sqrt{\ff{g}{l}}(t+C)
\end{split}
とできます。
$\theta$を指数関数として表すと、
\begin{split}
&\theta = e^{\pm i\sqrt{\ff{g}{l}}(t+C) } \EE
\therefore\, &\theta(t) = Ce^{\pm i\sqrt{\ff{g}{l}}}
\end{split}
となります。（$C=e^{i\sqrt{\ff{g}{l}}C}$としています。）
後は、オイラーの公式を利用して式変形を進めていきます。
このように、演算子法は、二階微分を一階微分に変形し、そこから微分方程式を解く手法です。
ラプラス変換による解法
特性方程式による解法も、演算子法による解法でも計算が面倒なため、もっと簡単に解きたいと思うのが、人情というものです。
その際に利用されるのが、ラプラス変換を使った解法です。
発展的な内容になりますが、ラプラス変換による微分方程式の解法を紹介します。
ラプラス変換を簡単に言うと、線形微分方程式の解が指数関数によって表せるという性質を利用した解法です。
微分方程式について学習が進んだ後にこの節を復習すれば、より理解が深まると思います。
前置きが長くなりましたが、ラプラス変換とは次のような無限積分を用いた変換のことです。
$f(t)$を$0 \leq t$で定義された関数とするとき、$s$の関数
\begin{eqnarray}
\mathscr{L}[f(t)] = F(s) = \int_0^{\infty} e^{-st} f(t) \diff t
\end{eqnarray}
を$f(t)$のラプラス変換と呼び、$\mathscr{L}[f(t)]$と書く。
早速、式(2)をラプラス変換してみましょう。（初めての方は、眺めるだけで大丈夫です。）
\begin{split}
&\mathscr{L}\left[ \ff{\diff^2 \theta}{\diff t^2} +\ff{g}{l}\theta \right] = \mathscr{L}[ 0 ] \EE
&\mathscr{L}\left[ \ff{\diff^2 \theta}{\diff t^2}\right] + \ff{g}{l}\mathscr{L}\left[ \theta \right] = 0 \EE
&\left( s^2\mathscr{L}\left[ \theta \right] \,-\, s\theta(0) \,-\, \dot{\theta}(0) \right) + \ff{g}{l}\mathscr{L}\left[ \theta \right] = 0
\end{split}
ここで、$\mathscr{L}\left[ \theta \right] \equiv \Theta(s) ,\, \theta(0)=\theta_0 $とし、$\dot{\theta}(0)=0$を代入します。
$\Theta (s)$について先ほどの式を整理すると、以下のようになります。
\begin{split}
s^2&\Theta(s) \,-\, s\theta_0+ \ff{g}{l}\Theta(s) = 0 \EE
\therefore \, & \Theta(s) = \ff{s}{s^2 + \left(\sqrt{\ff{g}{l}}\right)^2 }\cdot \theta_0
\end{split}
ところで、ラプラス変換の結果がいくつかのパターンに分類できることを利用し、あらかじめ計算結果をまとめた表が存在します。
この表をラプラス変換表と呼びます。
ラプラス変換表を見ると、下のような対応関係を見つけることができます。
\begin{array}{|c|c|}
\mathscr{L}[f(t)] & \mathscr{L}^{-1}[f(t)] \EE
\hline
\displaystyle\ff{s}{s^2+a^2} & \cos at \EE
\end{array}
ラプラス変換表の係数を比較することで、微分方程式の解を得ることができます。
この操作をラプラス逆変換と呼びます。
$\Theta(s)$をラプラス逆変換してやると、
\begin{split}
& \mathscr{L}^{-1}[\Theta(s)] = \theta(t) = \theta_0\cos\sqrt{\ff{g}{l}}t
\end{split}
このようになり、見事に$\theta(t)$を求めることに成功しました。
境界条件とラプラス変換表が手元にあれば、線形微分方程式を式変形のみで解くことができます。
特性方程式や演算子を使用する手法よりも簡単に微分方程式を解くことができます。
ところで、非線形微分方程式である式(1)の解はどのように表せるのでしょうか？
近似を使わない、振り子の周期の一般解についてはこちらの記事にて解説しています。
関連記事・参考記事
