Pythonによる計算物理（３）

本稿から具体的な数値計算の手法についての話に入る。まずは、古典力学の基本方程式であるNewtonの運動方程式を数値的に解く手法から。

古典力学と常微分方程式

古典力学は、基本的には質点の運動はNewtonの運動方程式で記述される。
Newtonの運動方程は、以下のように2階の常微分方程式の形である。

ここで、 $F(t,x(t),v(t))$ は外部から質点に作用する力である。
$t$ ：時間、 $x(t)$ ：質点の位置、 $v(t)$ ：質点の速度である。

一般に、 $F$ は時刻 $t$ と位置 $x(t)$ および $v(t)$ に依存する。
たとえば、重量の場合は $F=mg$ となる。
強制振動の場合は、 $F=F_0cos(\omega t)$ 。
バネの場合は、 $F=-kx(t)$ 。
空気抵抗の場合は、 $F=-bv(t)$ 。

といったように、作用する力の種類によって、 $F$ の関数形が決まる。

2階常微分方程式の変形

再び、Newtonの運動方程式に戻る。

2階の常微分方程式を数値計算するには、新たに速度 $v(t)=\frac{dx(t)}{dt}$ を導入することで、以下のような２つの1階の常微分方程式に分けると、数値解法を適用しやすくなる。

常微分方程式は、一般にN個の成分 $y_i(t)(i=1,...,N)$ それぞれに対する1階の常微分方程式として、以下のような形で表せる。

ベクトルの形で書くと以下となる。

ここで、ベクトル表記 $\boldsymbol{y}=(y_1,...,y_N)$ を用いている。
Newtonの運動方程式を一般的なベクトル表記の形式に当てはめると、以下のようになる。
- 1次元の場合　 $\boldsymbol{y}=(x, v)$ と定義できる。
- 2次元の場合　 $\boldsymbol{y}=(x, y, v_x, v_y)$ と定義できる。
- 3次元の場合　 $\boldsymbol{y}=(x, y, z, v_x, v_y, v_z)$ と定義できる。

常微分方程式の数値解法

基本的な考え方として、位置 $\boldsymbol{y}$ 、時間 $t$ を離散化する。
時刻 $t=t_n$ での位置ベクトル $\boldsymbol{y}(t_n)=\boldsymbol{y_n}$ とおく。
$\boldsymbol{y}$ を既知として、次の時刻[tex: t{n+1}]での位置ベクトル[tex: \boldsymbol{y}{n+1}]を以下のように逐次的に求める。

ここで、 $\Delta\boldsymbol{y}$ を見積もる方法として、次のようなものがある。
- 陽解法 - 陰解法 - 予測子-修正子法

陽解法

陽解法は、もっとも直接的な方法で、代表的なのがオイラー法である。

<基本的な考え方> - $t=t_n$ における関数 $\boldsymbol{y}(t)$ の接線で $\boldsymbol{y}(t)$ を近似。 - 接線は、 $\boldsymbol{y}(t_n+h)$ を $t_n=t$ の周りでテイラー展開し、[tex: h²]以上の項を無視することで得られる。

式で書くと以下となる。

Euler法（中点法）

オイラー法の考え方をより精度を高いものにした手法として、中点法がある。
以下に中点法の考え方と模式図を示す。

傾き $k_1$ で中間点 $t=t_n+h/2$ まで進み、そこで新たに傾き $k_2$ を見積もる。
傾き[yex: k_2]を使って、元の出発点から終点 $t=t_n+h$ まで進める。

式で書くと、次のようになる。

ポイント

新たな傾き $k_2$ は $t=t_n+\frac{h}{2}$ で計算するが、 $\boldsymbol{y}_{n+1}$ の計算は、元の出発点 $t=t_n$ に戻って、 $k_2$ を使って計算する。
1回の時間発展において、関数 $\boldsymbol{f}(t,\boldsymbol{y}(t))$ を2回計算することで、精度の次数が1つ向上する。