計算物理学: No.4

物理学を含む自然現象一般はしばしば微分方程式の形で表される。ここではその中でも基本的な常微分方程式の数値解法を扱う。

1 常微分方程式の数値解法

1.1 常微分方程式の正規形

常微分方程式の問題は、常に 1 階の微分方程式の組に置き換えることができる。例えば、 2 階の微分方程式

\frac{d^{2}y}{dx^{2}} + q(x)\frac{dy}{dx} = r(x)  (1)

は新しい変数 z を補助的に用いて、連立の 1 階微分方程式

\frac{dy(x)}{dx} = z(x)
\frac{dz(x)}{dx} = r(x) − q(x)z(x)  (2)

に書き直すことができる。一般に（多変数、高階の）常微分方程式の問題は、関数 y^{i}, i = 1, 2, …, m に対して、

\frac{dy^{i}(x)}{dx} = f^{i}(x, y^{1}, y^{2}, …, y^{m}),  i = 1, …, m  (3)

の正規形をした m 元連立 1 階微分方程式に帰着する。初期条件 y^{i}(0) の組を与えると、これを数値的に解くことができる。

1.2 Euler 法

刻み幅 h に対し、独立変数の離散的な値を

x_{n} = a + nh (n = 0, 1, 2, ⋯)  (4)

とする。常微分方程式を差分に直すと

\frac{dy^{i}}{dx} (= f^{i}(x_{n})) ≈ \frac{y_{n+1}^{i} − y_{n}^{i}}{x_{n+1} − x_{n}} = \frac{y_{n+1}^{i} − y_{n}^{i}}{h}  (5)

となる。これから x_{n} における y^{i}(x_{n}) の近似値を y_{n}^{i} とした時、次のステップ x_{n+1} での未知関数 y_{n+1}^{i} を

y_{n+1}^{i} = y_{n}^{i} + hf_{n}^{i},  f_{n}^{i} ≡ f^{i}(x_{n}, y_{n}^{1}, ⋯, y_{n}^{m})  (6)

によって次々と決めるのがオイラー法である。

1.3 Euler 法の誤差

y(x) のテイラー展開式を考えると、

y(x + h) = y(x) + hy′(x) + \frac{1}{2}h^{2}y′′(x) + ⋯  (7)

従って常微分方程式を差分に直す誤差は、

\frac{dy^{i}}{dx} = \frac{y_{n+1}^{i} − y_{n}^{i}}{h} ＋ O （h^{2})  (8)

計算する区間を a から b までとし、それを n 等分して計算を進めると h = (b − a)/n なので、 x = a から x = b までの n ステップの誤差は、

nO(h^{2}) = \frac{(b − a)}{h}O(h^{2}) = O(h)  (9)

になる。つまり、刻み幅を半分にすれば誤差は半分になる。

2 課題

課題 1.

\frac{dy}{dx} = 1 − y^{2}, y(0) = 0  (10)

の数値解をオイラー法によって求めよ。この問題の解析解 tanh x と比較せよ。刻み幅 h に対する誤差の振舞いを確認せよ。

課題 2.

単振動の微分方程式

\frac{d^{2}y}{dx^{2}} + y = 0  (11)

は 2 元 1 階常微分方程式系

\frac{dy}{dx} = z,
\frac{dz}{dx} = −y  (12)

に置き換えられる。

この数値解をオイラー法によって求めよ（初期条件は適当に選べ）。結果を解析解と比較せよ。刻み幅 h に対する誤差の振舞いを確認せよ。