ニュートン法で方程式の近似解を求める - 薬剤師のプログラミング学習日記

前回「トラフ値からkeを推定する」で定常状態のCss_minを求める式からkeを求めました。このときscipy.optimizeモジュールのnewton()を使いましたが、これはニュートン法（ニュートン・ラフソン法）というアルゴリズムを実装した関数です（引数に導関数を与えた場合はニュートン法、与えない場合は割線法で処理されます）。今回はこのアルゴリズムをPythonで書いてみようと思います。

ニュートン法とは
ニュートン法のアルゴリズム
数値微分
ニュートン法の実装
数値計算における注意点
参考文献

ニュートン法とは

方程式の解の一つを求める方法で、反復計算によってf(ｘ)＝０となるxを近似的に求めることができます。

ニュートン法のアルゴリズム

下図青のグラフのような $y=f(x)$ で表せる曲線があるとします。方程式の解は $f(x)=0$ となるxの値なので、曲線とｘ軸との交点ということになり、次の手順で求めることができます。

$f(x)=0$ の解に近い適当な初期値 $x_n$ を定める。
$(x_n, f(x_n))$ における $f(x)$ の接線を引く。
接線とx軸との交点を得る。このときのｘ座標 $x_{n+1}$ を新たな予測値とする。
2、3を一定の回数繰り返す。

f:id:enokisaute:20200317231632p:plain

$x_{n+1}$ と $x_n$ の関係についてみてみます。

$y=f(x)$ 上の点 $(x_n, f(x_n))$ における接線の方程式は次のように表せました。

$y - f(x_n) = f'(x_n)(x - x_n)$

この直線がx軸上の点 $(x_{n+1}, 0)$ を通るので、 $x=x_{n+1}, \ y=0$ を代入すると、

$-f(x_n) = f'(x_n)(x_{n+1}-x_n)$

これを整理すると次のような漸化式となり、上の1～4はこの式の $x_{n+1}$ を繰り返し求めることに相当します。

$\large{x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}}$

数値微分

上のアルゴリズムで $x_n$ を求めるには $f(x)$ の微分が必要となります。簡単に微分が求められるような関数の場合は良いのですが、複雑な関数の場合は数式から解析的に微分を求めるのは面倒です。
そこで、数値微分を使うとある程度の誤差を許容することにはなりますが、計算によって微分の近似値を求めることができます。

微分（導関数）は数式で書くと次のように定義されます。

$f'(x) =\displaystyle \lim_{h \to 0} \frac{f(x+h) - f(x)}{h}$

ここでは、この式をそのまま実装するのではなく、誤差をより小さくするために、ｘの前後(x+h)と(x-h)での関数fの差分（中心差分）から求めた傾きを微分の近似値とします。また、差分で用いる微小量のhには $10^{-4}$ を用いています。

h = 1e-4    # 0.0001
df = (f(x + h) - f(x - h)) / (2 * h)

これを使って $f(x)=x^2-2$ の微分を求めてみます。解析的には $f'(x)=2x$ となり、x=2での $f(x)$ の微分は4となります。

# 微分を求めたい関数
def func(x):
    return x ** 2 - 2
 
# 中心差分で微分を求める
def numerical_diff(f, x0, h=1e-4):
    return (f(x0 + h) - f(x0 - h)) / (2 * h)
 
print(numerical_diff(func, 2))  # 4.000000000004

解析的な値と一致はしませんが、とても小さい誤差で微分を求めることができます。ちなみに、ここでは深入りしませんがhの値は小さすぎても駄目です。丸め誤差というコンピュータによる計算で起こる誤差の影響により、かえって結果の誤差が大きくなってしまうという現象が起きます。

ニュートン法の実装

 
# (解を求める方程式, 初期値, 微分で用いる微小量, 許容(絶対)誤差, 反復の最大数)
def newton_method(f, x0, h=1e-4, tol=1.48e-8, maxiter=50):
    for cnt in range(1, maxiter+1):
        # 中心差分で微分を求める
        df = (f(x0 + h) - f(x0 - h)) / (2 * h)
        next_x = x0 - f(x0) / df
 
        # 反復回数, 解の近似値x, f(x)を表示
        print('{:2d}: x={:.15f}, f(x)={:.15f}'.format(cnt, x0, f(x0)))
 
        # 収束判定条件. 1ステップの差が許容(絶対)誤差内なら終了
        if abs(next_x - x0) < tol:
            break
        x0 = next_x
    return x0

ニュートン法自体シンプルなアルゴリズムであるため、簡単に書くことができます。
引数の微分の微小量については前の項で書いた通りです。許容(絶対)誤差tol、反復の最大数maxiterについてはScipy.optimizeのニュートン法のデフォルト値と同じにしました。結果の精度をもっと良くしたいという場合はtolをより小さな値（1e-12とか）にするといいと思います。
また、収束していく様子がわかるように途中print()でf(x)（0に近づく）も表示するようにしています。
では次のように初期値を与えて $x^2-2=0$ を解いてみます。

print(newton_method(func, 2))

結果は、次のようになりました。#は小数点以下16桁までの $\sqrt{2}$ の値です。

 1: x=2.000000000000000, f(x)=2.000000000000000
 2: x=1.500000000000500, f(x)=0.250000000001500
 3: x=1.416666666666719, f(x)=0.006944444444593
 4: x=1.414215686274509, f(x)=0.000006007304879
 5: x=1.414213562374690, f(x)=0.000000000004511
1.4142135623746899    # 1.4142135623730950

次は前回記事で書いたCss_minの式からトラフ値9μg/mLとしてkeを求めてみます。

def func(ke, trough=9):
    F = 1             # バイオアベイラビリティ
    D = 200           # 投与量(mg)
    vd = 0.46 * 63.6  # 分布容積(L/kg) * 理想体重(kg)
    ka = 0.36         # 吸収速度定数(hr^-1)
    tau = 12          # 投与間隔(hr)
 
    A = F * D * ka / vd / (ka - ke)
 
    return A * np.exp(-ke * tau) / (1 - np.exp(-ke * tau)) \
            - A * np.exp(-ka * tau) / (1 - np.exp(-ka * tau)) - trough
 
 print(newton_method(func, 0.1)) # 0.052556754839818746

scipy.optimizeのnewton()の値が0.05255675488882937だったので、このケースでは実用的には充分と言える程度には解の近似値を求めることができました。

数値計算における注意点

今回は解もだいたいわかっている簡単な例を扱い、実装も単純なものにしましたが、数値計算ではコンピュータによる計算特有の誤差を考えたり、どのような収束判定条件を設定するかなどの計算アルゴリズムとはまた別に考えなくてはいけない問題があります。
これらのことを踏まえた上で、関数の形を確かめたり、何度も条件を変えながら実行したりして、結果の信頼性を評価することが必要です。

最後に今回のコードです。

# 方程式. =0となる形で書く
def func(x):
    return x ** 2 - 2
 
 
# (解を求める方程式, 初期値, 差分での微小量, 許容(絶対)誤差, 反復の最大数)
def newton_method(f, x0, h=1e-4, tol=1.48e-8, maxiter=50):
    for cnt in range(1, maxiter+1):
        # 中心差分で微分を求める
        df = (f(x0 + h) - f(x0 - h)) / (2 * h)
        next_x = x0 - f(x0) / df
 
        # 反復回数, 解の近似値x, f(x)を表示
        print('{:2d}: x={:.15f}, f(x)={:.15f}'.format(cnt, x0, f(x0)))
 
        # 収束判定条件. 1ステップの差が許容(絶対)誤差内なら終了
        if abs(next_x - x0) < tol:
            break
 
        x0 = next_x
    return x0
 
 
if __name__ == '__main__':
    print(newton_method(func, 1.5))