8.1. 非線型方程式の解法

8.1.1. 概観

この章では非線型方程式

f(x)=0

の解を求める計算法を紹介します。

非線型方程式を解くにはSciPyの scipy.optimize モジュールを使います。このモジュールは、最適化（最小化・最大化）や最小自乗フィットなども含む大きなモジュールです。方程式の解を求める関数は Root Finding にまとまっています。

公式ドキュメントを見ると、多くのソルバーがあることがわかります。ソルバーの特徴は Root Finding に載っている表によくまとまっています。ソルバーは大きく分けると2つに分類されます：

  1. 区間を指定する方法（bracketing methods）

  2. 初期値を指定する方法

1.の方法は、指定した区間に解があれば必ず答えにたどり着くという特長があります。アルゴリズムによってたどり着くまでに必要なイテレーションの数（f(x) を実行する回数）が変わります。ただし、この方法は x が実数の場合にしか使えません。

x が複素数の場合、または、変数が2つ以上ある場合は、2.の方法しか使えません。初期値が解に近ければ収束は早いですが、必ずしも初期値に近い解にたどり着くとは限りません。また、収束しない場合もあります。 f'(x) や f''(x) を利用できるかで選ぶアルゴリズムが変わります。

1.の方法が使える場合には1.を使い、それ以外の場合には2.を使うという考えで基本的によいと思います。ただし、1.が使える場合でも、初期値を適切に与えられる、かつ f'(x) が利用できる場合には2.の選択肢もあります。1.の場合にはBrent法（method='brentq'）を使っておけば間違いないです。

8.1.2. 区間を指定する法

このカテゴリーでもっとも単純なものが bisection法（二分法） です。根を挟み込む区間を半分にする操作を繰り返します。始めに区間内に解が存在すれば必ず解にたどり着く安全な方法です。反復回数 n に対し、精度が 2^{-n} に比例して小さくなっていきます（10回で3桁）。

二分法の安定性を残しつつ収束性を改善したものが Brent法 です。二分法に加えて、割線法（後述）による区間の見積りも同時に行い、複数の条件からそのどちらかを採用します。

8.1.3. 初期値を指定する方法

こちらのカテゴリーの方法は、f(x) の微分が使えるかどうかで選ぶアルゴリズムが変わります。

f(x) のみ → secant法
f'(x) まで → Newton法（Newton-Raphson法とも呼ぶ）
f''(x) まで → Halley法

Newton法 は候補点 x_i において、関数を直線（f(x_i) の接線）で近似し、その直線がゼロと交わる点を次の候補とする方法です。 f'(x_i) が小さい場合（x_i が f(x) の極大・極小に近い場合）にこの方法は不安定となります。初期値の選び方がとても重要です。

secant法（割線法） は f(x_i) と f(x_{i-1}) を結ぶ直線がゼロと交わる点を次の候補とする方法です。Newton法で使用する接線の代わりに、2点を通る直線で代用している、とみることもできます。

8.1.4. scipy.optimize.root_scalarの使い方

1変数の場合は scipy.optimize.root_scalar 関数を使います。 公式ドキュメントによると、引数は以下のようになっています：

scipy.optimize.root_scalar(f, args=(), method=None, bracket=None,
                           fprime=None, fprime2=None, x0=None, x1=None,
                           xtol=None, rtol=None, maxiter=None, options=None)

第一引数fに f(x) の値を返す関数を指定します。もし、その関数が x 以外のパラメータを引数として持つ場合、例えば func(x, param1, param2) という関数の場合、param1 と param2 の値を args=(param1, param2) としてroot_scalar関数に与えます。これについては、scipy.integrateの使い方 にも説明があります。

利用可能なソルバー（method）は公式ドキュメント Root Finding にまとまっています。 以下に、表を引用して示します：

../_images/solvers.png

Domain of fは関数の値が実数か複素数か、Bracket?は挟み込む方法（区間を指定する方法）かどうか、Derivatives?は微分を使うかどうか、Solversは解法の名前、Convergence Guranteed?は収束が保証されているか、Rate(s)は収束の速さ（大きいほど速い）です。 詳細は公式ドキュメント Root Finding を参照してください。

Brent法を使う場合は method='brentq' として、区間を bracket=(xmin, xmax) のように与えます。Newton法の場合は、method='newton' として、導関数 fprime と初期値 x0 を与える必要があります。それ以外の引数や別のmethodの使い方については公式ドキュメント scipy.optimize.root_scalar を参照してください。

8.1.5. 多変数の場合

変数が2つ以上ある連立非線形方程式の場合を考えます。

\bm{f}(\bm{x}) = 0

ここで、関数 \bm{f} および 変数 \bm{x} は N 次元ベクトル、 \bm{f}=(f_1, \cdots, f_N) および \bm{x}=(x_1, \cdots, x_N) です。

この場合は区間を指定して挟み込む方法は使えません。初期ベクトル \bm{x}_0 を指定して、勾配を使って更新する方法になります。一般に、初期値に近い解に収束するかどうか、あるいはいずれかの解に収束するかどうか、は保証されません。

SciPyでは scipy.optimize.root 関数で多次元連立方程式を解きます。公式ドキュメントによると、この関数の引数は以下のようになっています：

scipy.optimize.root(func, x0, args=(), method='hybr', jac=None, tol=None,
                    callback=None, options=None

1変数の場合と違い method にデフォルト値が設定されていて、method='hybr' はPowellのハイブリッド法と呼ばれる方法です。第5引数の jac はヤコビアンです。ヤコビアンの表式が分かっている場合にはその関数を jac に与えます。この引数を与えない場合には、ヤコビアンを \bm{f}(\bm{x}) から数値的に見積もるので、その分計算量が多くなります。詳細は公式ドキュメント scipy.optimize.root を参照してください。