今週は、運動方程式に復元力$F=-k{x}$(${x}=0$に向けて戻そうとする力${x}>0$なら負の向きの力、${x}<0$なら正の向きの力が加わる。つまりどっちにしても、${x}=0$に向かうような力である。よってこれを「復元力」と呼ぶ。)を加えた、 \begin{equation} m\left({\mathrm d\over\mathrm dt}\right)^2 {x}= -K{\mathrm d\over\mathrm dt}{x}-k{x}\label{Fkvkx} \end{equation} を解いてみよう(重力は考えないことにする)。
この運動方程式があらわす現象は「空気抵抗を受けながら運動するバネ振り子」である。では微分方程式を解く前に、どんな運動ができるか想像してみよう。
というわけで、今日はまず下のプログラムをしばらく動かして運動の様子をみた。
初期位置:y0=0
初速度:vy0=4
バネ定数:k=4
空気抵抗:K=1
やってみて見つかったかな?
バネ定数$k=4$、空気抵抗の係数$K=6$で初速度0、初期位置2で運動させると、
のような運動が起こる。これは確かに「振動しながらどんどん遅くなる」(ちなみに「減衰振動」と呼ぶ)。
バネ定数$k=4$、空気抵抗の係数$K=2$で初速度0、初期位置2で運動させると、
のように、全く「振動」ではない、ただ「原点に止まる」という運動になる(これは「過減衰」と呼ぶ)。
「どのあたりに二つの境目があるか?」とやってみるのも面白い。
この違いがどのように出るのか?
もちろん、空気抵抗が大きくなれば振動できなくなるというのは気持ちとしてはわかるだろうが、そこを微分方程式から納得したい。また、実はこの方程式の解は2種類ではなく3種類にわかることが、実際に方程式を解いてみるとわかる。
では、方程式を解きながら考えてみよう。
運動方程式 \begin{equation} m\left({\mathrm d\over\mathrm dt}\right)^2 {x}= -K{\mathrm d\over\mathrm dt}{x}-k{x} \end{equation} から、例によって特性方程式を作ると、$m\lambda^2 +K\lambda +k=0$となる。これの解は(二次方程式の解の公式を用いて) \begin{equation} \lambda_\pm={-K\pm\sqrt{K^2-4mk}\over 2m}=-{K\over 2m}\mp{{\sqrt{K^2-4mk}}\over 2m} \end{equation} となる。ここから、$K^2-4mk$が負の場合、0の場合、正の場合の三つに分けて考える。
このように「場合分け」をしなければいけないというのが、運動に種類がある(予想では2種類だったが実は3種類だった)理由である。
$K^2-4mk<0$の場合この場合は$\lambda$は複素数になる。$\omega={\sqrt{4mk-K^2}\over 2m}$という定数($\omega$は実数である)を定義し、 $\lambda_\pm=-{K\over 2m}\mp\mathrm i\omega$と書くことにする。
こうして解を \begin{equation} {x}= C_+ \mathrm e^{-{K\over 2m}{t}+\mathrm i\omega{t}} +C_- \mathrm e^{-{K\over 2m}{t}-\mathrm i\omega{t}}\label{gensuisindou} \end{equation} と表すことができる。一見複素数であるが例によって係数を操作して、$C_+=C,C_-=C^*$とすることで \begin{equation} {x}= \mathrm e^{-{K\over 2m}{t}}\left( C\mathrm e^{\mathrm i\omega{t}} +C^* \mathrm e^{-\mathrm i\omega{t}}\right) \end{equation} が実数解となる。
つまり、この式をちゃんと計算すると、実は虚数$\mathrm i$はどこにもいないのだ!
そのことを確認するために、$C=A+\mathrm iB,C^*=A-\mathrm iB$などと置いて計算をしてみる。
オイラーの公式$\mathrm e^{\mathrm i\theta}=\cos\theta+\mathrm i\sin\theta$などを使って、結果を三角関数で表現すると以下のようになる($A,B,D,\alpha$は実数の定数)。 \begin{equation} {x}=2 \mathrm e^{-{K\over 2m}{t}}\left( A\cos \omega{t} +B\sin\omega{t}\right) = D\mathrm e^{-{K\over 2m}{t}} \cos \left( \omega{t}+\alpha\right)\label{gensui} \end{equation}
下の動くグラフに$K^2-4mk<0$になる値(ちなみに$m=1$で固定であるから、$K^2<4k$ということ)を代入してみよう。二つの解(青と赤)がプロットされる。線形結合を取ることでいろんな解が出ることもわかる。
$K^2-4mk=0$の場合この場合、$\lambda_+=\lambda_-$となる。よって添字は取って、$\lambda=-{K\over 2m}$と書こう。重解が出る場合であるから、解は以下の通り。 \begin{equation} {x}=\left(C_1{t}+C_0\right)\mathrm e^{\lambda{t}}\label{rinkai} \end{equation} この解「臨界振動」は、$K^2-4mk<0$と$K^2-4mk>0$のちょうど境目にあたる。
下の動くグラフに$K^2-4mk=0$になる値($K^2=4k$ということ)を代入してみよう。
$K^2-4mk>0$の場合この場合は単純に、 \begin{equation} {x}= C_+ \mathrm e^{\lambda_+{t}} +C_- \mathrm e^{\lambda_-{t}}\label{kagensui} \end{equation} が解である。$\lambda_\pm$はどちらも負の数になるから、この二つの解のどちらも「指数関数的に減衰する解」である。
これも下のグラフにいろんな数値を入れてみよう。
A=1
B=1
K=1
k=1
次に「定数係数」という条件を外して考えることにする。
一般的な一階線形非斉次微分方程式は、$p({x})$と$q({x})$を既知の${x}$の関数として、 \begin{equation} {{\mathrm d\over\mathrm dx}}f({x})+p({x})f({x})=q({x})\label{DEpq} \end{equation} と書くことができる。$f({x})$が今から求めようとしている「未知の関数」である。より一般的には \begin{equation} r({x}) {{\mathrm d\over\mathrm dx}}f({x})+p({x})f({x})=q({x}) \end{equation} という形も考えられるが、この式の両辺を$r({x})$で割って整理したものがこの式だと思えばよい(もちろんこの計算は$r({x})\neq0$の領域でのみ可)。この方程式は \begin{equation} \left( {\mathrm d\over\mathrm dx} +p({x})\right)f({x})=q({x}) \end{equation} とも書ける。ここで、 $$ \left({\mathrm d\over\mathrm dx}-A\right)\left(\mathrm e^{A{x}}F({x})\right)= \mathrm e^{A{x}}{\mathrm d\over\mathrm dx} F({x}) $$ を思い出す。ここでは${\mathrm d\over\mathrm dx}$の後には数$-A$ではなく関数$p({x})$がついているわけだが、その場合でも上の真似をして、$f({x})=\mathrm e^{-P({x})}F({x})$と置き直すことで \begin{equation} \left( {\mathrm d\over\mathrm dx} +p({x})\right)\left(\mathrm e^{-P({x})}F({x})\right) =\mathrm e^{-P({x})}{\mathrm d\over\mathrm dx} F({x}) \end{equation} とできないか(微分演算子と$\mathrm e^{-P({x})}$を交換することで$p({x})$を「消去」できないか)と考える。
ここで$\int \mathrm dx p({x})=P({x})+C$すなわち$P({x})$が$p({x})$の原始関数の一つであるとすれば、微分${\mathrm d\over\mathrm dx}$の結果が${\mathrm d\over\mathrm dx}\mathrm e^{-P({x})}=-p({x})\mathrm e^{-P({x})}$となって、ちょうど$p({x})$の項を打ち消す項が出てきて、 \begin{equation} \begin{array}{rl} \overbrace{-p({x})\mathrm e^{-P({x})}F({x})}^{\tiny\left({\mathrm d\over\mathrm dx}\mathrm e^{-P({x})}\right)F({x})} + \mathrm e^{-P({x})} {\mathrm d\over\mathrm dx} F({x}) +p({x})\mathrm e^{-P({x})}F({x}) =&q({x}) \\[3mm] \mathrm e^{-P({x})}{\mathrm d\over\mathrm dx} F({x})=&q({x}) \end{array}\label{ddxkoukan} \end{equation} が解くべき方程式となる。
こうして、$p({x})$の原始関数$P({x})$を使うことで \begin{equation}\left( {\mathrm d\over\mathrm dx} +p({x})\right)f({x})=q({x})~~~\to~~~ {\mathrm d\over\mathrm dx} F({x})=q({x})\mathrm e^{P({x})} \end{equation} と式を書き直せたので、後はこれを解く。右辺を${x}$で積分することができれば \begin{equation} F({x})= \int \mathrm dx \left(q({x})\mathrm e^{P({x})}\right) \end{equation} となり、 \begin{equation} f({x})=\mathrm e^{-P({x})}\int \mathrm dx \left(q({x})\mathrm e^{P({x})}\right) \end{equation} が一般解である。
と、ここまでで時間切れなので、具体的計算と、このやり方の別の方法は来週に。
青字は受講者からの声、赤字は前野よりの返答です。