非線形方程式をpythonで解く

非線形方程式

ある方程式 $y=f_{(x)}$の1次導関数がxに関する1次式でないとき、その方程式は非線形である。
というわけで、皆さんご存知の
$$ y = ax^2 + bx + c $$
も非線形方程式です。この子は $ y=0 $ なら中学生で習う二次方程式の解の公式 $ \frac{-b\pm\sqrt{b^2-4ac}}{2a} $で解析解が求められますが、
例えばこんなのだともう解析解を求めるのがおっくうになってしまいます。

$$ y = x^2 - \sqrt{x+3} $$

というわけで多くの非線形方程式は近似解として数値解析を用いた数値解を使うことがほとんどです。
今回は非線形方程式をPythonで解いてみようと思います。

解き方

scipy.optimize.newton(Newton-Raphson Method)

$ x_i $ における $ y = f_{(x)} $ の接線は
$$ y - f_{(x_i)} = f_{(x_i)}' (x - x_i) $$
接線とx軸の交点である $ y = 0 $をとった
$$ x = x_i - \frac{f_{(x)}}{f_{(x)}'} $$
が許容される誤差$ \varepsilon $ より小さくなったら解となり、満たさない場合は$x$を次の$x_i$として計算を続けます。

scipy.optimizeのnewtonは関数とその導関数を与えればNewton-Raphsonで計算し、導関数を与えない場合はSecant Methodで計算します。Secant Methodは導関数のかわりに有限差分を用いたもので収束性はNewton-Raphsonより良くないです。ただ、導関数を求めるのは非常に面倒。実用上求めているのは教科書的な式をそのまま投げ込んで解けるもの。というわけでここではsecantを使います。

from scipy import optimize
optimize.newton(func, estimate)

scipy.optimize.fsolve

詳しくは触れませんが、fsolveはFortranのminpackライブラリのラッパーのようです。ヤコビ行列を使って解くみたいですね。

from scipy import optimize
optimize.fsolve(func, init)

fsolveの真価はnewtonとは違い多変数関数に対応できることなので、実際はnewtonとは競合する関係ではないです。ここでは一応比較として。

Interpolate Method

この手法はちゃんとした名前が他にあるかもしれませんが、私が知らないので適当な名前をつけて呼んでいる手法で、先に解空間を作っておいて補間で解を導くものです。
事前に解空間をscipyのinterpolateで作っておき、あとから補間値だけを取り出します。解空間を保持しておかないと遅くなってしまうのでclassとあわせることが必須です。（interpolateの戻りがイテレータとかなら変わらないかも。不勉強。）

from scipy import optimize

class MyClass:
    def __init__(self):
        x_array = np.arange(0.001, 1.0, 0.00001)
        y_array = func(x_array)
        self.y_inter = interpolate.interp1d(y_array, x_array, kind='linear', fill_value='extrapolate')

    def polate(self, y):
        return self.y_inter(y)

解いてみた

以下の3つの方程式を試した。
$$ y = x\sin{x} $$
$$ y = x^2 - \sqrt{x + 3} $$
$$ \frac{A_e}{A_t} = (\frac{2}{\gamma + 1})^{\frac{1}{\gamma - 1}} (\frac{P_c}{P_e})^{\frac{1}{\gamma}} \sqrt{\frac{\gamma + 1}{\gamma - 1} (1 - (\frac{P_e}{P_c})^{\frac{\gamma - 1}{\gamma}})} $$
3つめの式は超音速ノズルにおける開口比とノズル入口圧力と出口圧力の比を関係させた式であり、$\gamma$はノズル内の流体の比熱比。開口比 $ \frac{A_e}{A_t} $を与えたときの圧力比 $ \frac{P_r}{P_c} $を求めます。というかこれを解くためにこの記事書いてます。

import time
import numpy as np
from scipy import optimize as opt
from scipy import interpolate

y = 3.49
gamma = 1.1612

# optimize.newton
class MyClass1:
    def opt_ntrs(self):
        def eq_sin(x):
            return x * np.sin(x) - y
        def eq_root(x):
            return x**2 - np.sqrt(x + 3) - y
        def eq_eps(PePc):
            # y = 開口比            
            return ((2.0 / (gamma + 1.0))**(1.0 / (gamma - 1.0))) / (PePc**(1.0 / gamma) * np.sqrt(((gamma + 1.0) / (gamma - 1.0)) * (1.0 - PePc**((gamma - 1.0) / gamma)))) - y

        # return opt.newton(eq_root, 0.0)
        # return opt.newton(eq_sin, 0.0)
        return opt.newton(eq_eps, 0.001)

# optimize.fsolve
class MyClass2:
    def opt_fsolve(self):
        def eq_sin(x):
            return x * np.sin(x) - y
        def eq_root(x):
            return x**2 - np.sqrt(x + 3) - y
        def eq_eps(PePc):
            # y = 開口比
            return ((2.0 / (gamma + 1.0))**(1.0 / (gamma - 1.0))) / (PePc**(1.0 / gamma) * np.sqrt(((gamma + 1.0) / (gamma - 1.0)) * (1.0 - PePc**((gamma - 1.0) / gamma)))) - y

        # return opt.fsolve(eq_root, 0.0)[0]
        # return opt.fsolve(eq_sin, 0.0)[0]
        return opt.fsolve(eq_eps, 0.001)[0]

# interpolate
class MyClass3:
    def __init__(self):
        def eq_sin(x):
            return x * np.sin(x)
        def eq_root(x):
            return x**2 - np.sqrt(x + 3)
        def eq_eps(PePc):
            # y = 開口比            
            return ((2.0 / (gamma + 1.0))**(1.0 / (gamma - 1.0))) / (PePc**(1.0 / gamma) * np.sqrt(((gamma + 1.0) / (gamma - 1.0)) * (1.0 - PePc**((gamma - 1.0) / gamma))))

        x_array = np.arange(0.001, 1.0, 0.00001)
        # y_array = eq_root(x_array)
        # y_array = eq_sin(x_array)
        y_array = eq_eps(x_array)
        self.y_inter = interpolate.interp1d(y_array, x_array, kind='linear', fill_value='extrapolate')

    def polate(self):
        return self.y_inter(y)

if __name__ == '__main__':
    start = time.time()
    my1 = MyClass1()
    for i in range(10000):
        a = my1.opt_ntrs()
    print(my1.opt_ntrs())
    end = time.time()
    print('time', end - start, 'sec')

    start = time.time()
    my2 = MyClass2()
    for i in range(10000):
        a = my2.opt_fsolve()
    print(my2.opt_fsolve())
    end = time.time()
    print('time', end - start, 'sec')

    start = time.time()
    my3 = MyClass3()
    for i in range(10000):
        a = my3.polate()
    print(my3.polate())
    end = time.time()
    print('time', end - start, 'sec')

一回の計算だと短くて評価しづらいので1万回ぐらい計算。式によって解領域が違うので探索開始値と解空間用意範囲は適宜修正。
加えて $ y = xsin{x} $ だけは $ y > 1 $ だと解が求まらないので $ y = 0.6 $ を使う。

結果

計算時間 [sec]	newton	fsolve	interpolate
$y=x\sin{x}$	0.44	0.98	0.45
$y=x^2-\sqrt{x+3}$	0.43	1.38	0.43
$ \frac{A_e}{A_t} = (\frac{2}{\gamma + 1})^{\frac{1}{\gamma - 1}} (\frac{P_c}{P_e})^{\frac{1}{\gamma}} \sqrt{\frac{\gamma + 1}{\gamma - 1} (1 - (\frac{P_e}{P_c})^{\frac{\gamma - 1}{\gamma}})} $	0.89	3.05	0.44

時間は3回まわした平均値。
fsolveは圧倒的に遅いですね。newtonとinterpolateの勝負になりました。
interpolateは与えられた値に対して補間値を返すだけなので、どんな式を使おうとも計算時間は変わらないものかと思います。（つまり戻りはイテレータということだろう）
シンプルな式ならほとんど変わりませんが、式が複雑になると指定の値だけを計算するinterpolateに対してnewtonは探索が入るためやや不利なようです。

解	newton	fsolve	interpolate
$y=x\sin{x}$	2.93576651985	2.93576651985	2.93576651984
$y=x^2-\sqrt{x+3}$	2.41170199519	2.41170199519	2.41170199518
$ \frac{A_e}{A_t} = (\frac{2}{\gamma + 1})^{\frac{1}{\gamma - 1}} (\frac{P_c}{P_e})^{\frac{1}{\gamma}} \sqrt{\frac{\gamma + 1}{\gamma - 1} (1 - (\frac{P_e}{P_c})^{\frac{\gamma - 1}{\gamma}})} $	0.0566476694712	0.0566476694712	0.0566476697556

newton（というかsecant）とfsolveの解は一致していて、interpolateがどこまで近いかといった雰囲気ですが、
シンプルな式なら$e^{-10}$、ノズルの式でも$e^{-8}$ぐらいの精度が出ています。

これなら割とinterpolate優位に見えますが、newtonに与える解推定値の(収束に関する)厳しさに比べてinterpolateで用意する解範囲がかなりシビアでした。与えた解範囲に対して求める解が外挿になると合わないのは当たり前ですが、解範囲を広くとっておいて内挿にしようとしても合わなくなる時があります。

たとえば2番目のルートを持つ式を解いた時。
polateに対して-2から20までを0.01刻みで解範囲を与えてやると

のように求めるyが3以下だと値が求まっていません。newtonではうまくいってますね。しきい値の3は$\sqrt{x+3}$の3かと思いきや違いました。
与える解範囲を0から5にしてやるとご覧の通り。

数学には強くないので検証できていませんが、2次方程式の解の両方が解範囲に入っているとこうなるんじゃないかと。

3番目のノズルの式に対しても。

まとめ

よっぽど複雑な式ですごい数のループをまわすのでなければnewtonが良い。
解範囲がある程度既知で絞れているならinterpolateを使うと高速化できる、ぐらいのイメージ。

ただし、newtonは1変数関数のみの対応なので多変数関数になるとfsolveになる。
一応2変数までならpolateでもinterp2dで対応できる。