微分方程式の数値解法(オイラー法)

ある種の微分方程式については， 解析解を得る具体的な方法が確立している。 しかし， そのような場合でも，手順が膨大で実質的には人が手(と紙と鉛筆)で 解くことが不可能であったり， 一般的には，解析解が得られないのが普通である。 最近ではコンピュータの能力が向上したために， 比較的簡単な微分方程式については 解法の手順を予めコンピュータに与えておき， 「数式処理システム」 ¹⁰を用いることでコンピュータによっても 解析解を得ることができるが， 一般的には，コンピュータを利用した数値解に頼らざるを得ない。

微分方程式を数値的に解くとはどのようなことだろうか。 また，どのようにして数値解を得れば良いのだろうか。

具体的な例として， 地表からの高さ $h$ の地点から，水平となす角 $\theta$ ， 速さ $v_0$ で投げられた質量 $m$ の物体の運動について 考えてみよう(図)。 ここで，物体には一様な重力と，空気の粘性による抵抗力がはたらくとしよう。 空気抵抗の力の大きさを見積もることは， 物体の速さや大きさなどに依存する複雑な要素が絡んで， 一般的には非常に難しい問題であるけれども， 小石程度の物体が空気中を自由落下する場合については， 抵抗力の大きさは概ね，物体の速さ $v$ の二乗に比例することが知られている。 抵抗力の向きは，常に速度 $\vec{v}$ と逆向きである。 そこで，空気中を運動する物体のモデルとして，物体にはたらく力 $\vec{f}$ が 次式で与えられるとしよう。

$$\vec{f} = m\vec{g} - \gamma v\vec{v}$$ (10)

ここで，右辺第一項は一様な重力， 第二項は大きさが速さの二乗に比例する空気抵抗の力で， $\vec{g}$ は重力加速度(ベクトル)，$\gamma$ は比例定数である。

図: 放物体の軌跡

図 のように $x$-座標と $y$-座標を設定すると ¹¹，()式の運動方程式は両辺を $m$ で割り， $x$-，$y$- の各成分に分けることで

$$\left. \begin{array}{rcl} a_x(t) &=& f_x(t)/m \\ a_y(t) &=& f_y(t)/m \end{array}\right\}$$ (11)

となる。$a_x,\;a_y$ は加速度ベクトル $\vec{a}$ の成分である。 また，()式の力 $\vec{f}$ の成分は，

$$\left. \begin{array}{rcl} f_x(t) &=& - \gamma v(t) v_x(t) \\ f_y(t) &=& - m g - \gamma v(t) v_y(t) \end{array}\right\}$$ (12)

となる。 ここに， $g \simeq 9.8 ({\rm m}/{\rm s}^2)$ は重力加速度の大きさ， $v_x,\;v_y$ は速度ベクトル $\vec{v}$ の成分，また，速さ $v$ は

$$v(t) = \sqrt{v_x(t)^2 + v_y(t)^2}$$ (13)

である。 $f_x$ と $f_y$ のそれぞれは， $v$ を通して $v_x$ と $v_y$ の両方に依存することに注意しよう。

題意より，「初期条件」は，時刻 $t = 0$ のとき，

$$\left. \begin{array}{rclrcl} x(0) &=& 0,& v_x(0) &=& v_0 \co... ...\ y(0) &=& h,& v_y(0) &=& v_0 \sin \theta \end{array}\right\}$$ (14)

である。 問題は，以上の条件より， 任意の時刻 $t$ における $x(t)$ と $y(t)$ を得ることである。 空気抵抗が無い場合，すなわち，$\gamma = 0$ のときには， この問題はコンピュータを使わずとも解析的に簡単に解くことができて，その解は

$$\left. \begin{array}{rcl} x(t) &=& v_0 t \cos \theta \\ y(t... ... \frac{1}{2} g t^2 + v_0 t \sin \theta + h \end{array}\right\}$$ (15)

であることが知られている。 しかし，$\gamma \neq 0$ の場合にはこの問題は簡単には解析的に解くことができない。 そこで，コンピュータの出番となるわけである ¹²。

この問題を数値的に解くには， 4.4.1 で説明した 微分法の導入と逆の手順を辿ればよい。 微分を定義している ()式に着目していただきたい。 ()式のように， 微分係数 $dx/dt$ は比 $\Delta x/{\Delta t}$ の ${\Delta t}\rightarrow 0$ の極限値として定義されているので， 十分に小さな ${\Delta t}$ をとれば， たとえ ${\Delta t}$ の値が有限であっても

$$\frac{dx}{dt} \simeq \frac{\Delta x}{{\Delta t}}$$ (16)

と近似でき， 原理的には， ${\Delta t}$ の値を小さくとる程近似は良くなることがわかる。 そこで，適当な小さな値として ${\Delta t}$ を定めて， $0,\; {\Delta t},\; 2 {\Delta t},\; \cdots$ といった離散的な時刻

$$t_n = n {\Delta t},\qquad ( n = 0, 1, 2, \cdots )$$ (17)

についてだけの位置や速度などを考えることにしよう。 速度が位置の時間微分であり， 加速度が速度の時間微分であることに注意して， ()式の近似を用いれば， 速度と加速度の $x$-成分 $v_x(t_n)$ と $a_x(t_n)$ は

$$\left. \begin{array}{rcl} v_x(t_n) &\simeq& \{x(t_n) - x(t_{... ...eq& \{v_x(t_n) - v_x(t_{n-1})\}/{\Delta t} \end{array}\right\}$$ (18)

と近似できる。 $y$-成分 $a_y(t_n)$ と $v_y(t_n)$ についても同様である。 () 式を用いて，()式の運動方程式を整理すると，結局， 解くべき近似的な方程式は

$$v_x(t_{n+1}) = v_x(t_n) + \frac{1}{m}f_x(t_n){\Delta t}$$ (19)

$$v_y(t_{n+1}) = v_y(t_n) + \frac{1}{m}f_y(t_n){\Delta t}$$ (20)

$$x(t_{n+1}) = x(t_n) + v_x(t_{n+1}){\Delta t}$$ (21)

$$y(t_{n+1}) = y(t_n) + v_y(t_{n+1}){\Delta t}$$ (22)

と書ける。 このような方程式は「差分方程式」と呼ばれているが， 微分方程式を差分方程式で近似する方法は一意的でない点に注意しよう。 () $\sim$ ()で採用したのはその中でも最も簡単なもので， オイラー法と呼ばれているものである。

オイラー法によって常微分方程式の数値解を 具体的に計算する手順(アルゴリズム)についてまとめてみよう。

1. 与えられた初期条件により， $n=0$ の場合について ()式と()式の右辺を計算する。
2. 初期条件で与えられた $x(0), y(0)$ と の結果と ()，()式から定まる $v_x(t_1), v_y(t_1)$ を ()，()式に代入して， $x(t_1), y(t_1)$ を得る。
3. 初期条件の代りに， 直前の段階で得られた $x(t_n), y(t_n), v_x(t_n), v_y(t_n)$ を用いて .，.と同様の手順で， $x(t_{n+1}), y(t_{n+1}), v_x(t_{n+1}), v_y(t_{n+1})$ を得る。
4. 必要なだけ，手順 . を繰り返す。

このアルゴリズムに従ったプログラムの例を章末に付録として付けておくので， 興味のある人は参照して欲しい。

図: 空気抵抗がある場合の放物体の軌跡

実際にプログラムを実行して， グラフにプロットした結果が図に示されている。

ここで述べた方法は直ちにケプラー問題にも適用することができるし， さらに自由度が多い場合についても，基本的には同様である。 多くの粒子からなる体系にたいして， 古典力学の運動法則に従って粒子を動かして， その体系の様々な性質を調べる計算機シミューレーションの手法を 「分子動力学法」という。 1980年代以後，スーパーコンピュータの発達と連動して， 分子動力学法についての新しい方法論が次々と提案され， この方法を用いて調べることのできる現象や分野が大きくひろがってきた。 分子動力学法の適用範囲は， 「物性物理学」でミクロな原子や分子の集合体としての物質の性質を調べる場合から， 「宇宙物理学」で銀河の集合としての超マクロな銀河団の形成の過程を調べる場合に 至るまで非常に広い分野に及んでいる。