逐次代入法、Wegstein法 - モデリングと数値計算（主に化学工学向け）

はじめに

実際に使う機会は少ないと思われるが、基礎として知っておいた方がよいと思われる逐次代入法、Wegstein法について触れ、 Pythonで計算した例を示す。

逐次代入法

方程式を次の形にして、

$x = f(x)$

解の最初の予測値を右辺に代入して計算し、得られた値を第2の予測値として再び右辺に代入して計算、を繰り返し近似解を得る数値解法。新しい解と1つ前の解の差が予め決めた許容値を下回るまで計算を繰り返す。

例えば次の方程式の解を求める場合は（解は $x=1.25$ )、

$\displaystyle{ \frac{1}{x} = 0.8 }$

まず次の形に式変形し、

$\displaystyle{ x = \frac{1}{x} - 0.8 + x }$

例えば次の手順で解が求まる。

x1 = 2  # 最初の予測値
err = 1e10  # 計算終了判定用

while err > 1e-5:  # 計算終了判定
    x2 = 1.0/x1 - 0.8 + x1  # x = f(x)
    err = abs(x1 - x2)
    x1 = x2

print(x1)

Pythonコード実行すると次のとおり解が得られる。

1.250005085698445

Wegstein法

逐次代入法は $x=f(x)$ の中の $f(x)$ が次の条件でないと数値が収束しない。

$|f(x)'|<1$

これを解決する方法がWegstein法である。

例えば次の問題を考える。

$x^2 = 2$

逐次代入法で解ける形に式変形すると

$x = x^2 + x - 2$

$f(x)=x^2+x-2$ であり、その微分は $f(x)'=2x+1$ となり $|f(x)'|<1$ を多くの場合で満たさない。逐次代入法を試してみると、

x1 = 2

for i in range(5):
    x2 = x1**2 + x1 - 2
    x1 = x2
    print(x1)

結果は次の通り発散していく。

4
18
340
115938
13441735780

Wegstein法は最初の予測値から第2の予測値を求めるところまでは同じだが、第3以降の予測値は次式により算出する。

$x_{k+1} = t\times f(x_k) + (1-t)\times x_k \qquad (k=2, 3, \cdots)$

ここで、

$\displaystyle{ t = \frac{1}{1-s},\qquad s=\frac{f(x_k) - f(x_{k-1})}{x_k - x_{k-1}} }$

先程の問題 ( $x^2 = 2$ ) は次のように解ける。

def f(x):
    return x**2 + x - 2  # f(x)

x1 = 2  # 最初の予測値
x2 = f(x1)  # 第2の予測値
err = 1e10

while err > 1e-5:
    s = (f(x2) - f(x1))/(x2 - x1)
    t = 1.0/(1.0 - s)
    x3 = t * f(x2) + (1.0 - t) * x2
    err = abs(x3 - x2)
    x1 = x2
    x2 = x3

print(x2)