非線形方程式を解く¶
本稿では SciPy を用いて非線形方程式を数値計算で解く方法について記す。
まず、方程式を解くためのツールはモジュール scipy.optimize に色々とある。試しに ドキュメントのリファレンス から Root
finding という文字列を検索してみるとよい。問題となっている方程式の特徴を活かして、どのソルバーを採用するのかを検討するのが一般的には重要なのだが、ここではあまり細かいところには立ち寄らず、とにかくそれらしい根が求まればよいという態度で通す。
単純な多項式については、専用の計算法がある。本稿で例示する。
なお、線形連立方程式のそれについてはページを改めて述べる。
Newton-Raphson 法¶
関数 scipy.optimize.newton は Newton-Raphson 法により、与えられた方程式の根の計算をする。この計算法については有名かつ基本的なものなので、忘れない自信があるのでここには記さない。
コード的な手順は次のとおりとなる。
根を求めたい式を Python の関数の形式で定義する。本説明では
funcとする。ここが面倒なのだが、別途
funcの解の近似値root_guessを用意しておく。関数
newtonに引数として上述のものを渡して呼び出す。
方程式 \(x - 2 \sqrt{x - 1} = 0\) を解くコード例を示す。暗算で求まる例で申し訳ない。
#!/usr/bin/env python
"""newton.py: Demonstrate equation solvers of SciPy.
"""
from scipy.optimize import newton
import numpy as np
# pylint: disable=invalid-name
def func(x):
"""A function."""
return x - 2 * np.sqrt(x - 1)
# Solve func(x) == 0.
root_guess = 1.5
print("Root:", newton(func, root_guess))
実行結果は私の環境では次のようになった。もちろん推定根の与え方で戻り値が揺れるが、だいたい 10 のマイナス 7, 8 乗の精度で問題ない。
Root: 1.99999999195
ここでは採り上げないが、関数 newton にはキーワード引数 fprime や
fprime2 で func の一次、二次の導関数を別途指定できる。導関数を指示する場合、解計算のアルゴリズムと収束率が良いほうに変わるようだ。
多項式方程式¶
SciPy というよりは NumPy だけで多項式の定義と根の計算が実現できるのだが本稿に記す。
コード的な手順は次のとおりとなる：
一変数多項式をクラス
numpy.polynomial.polynomial.Polynomialで定義する。コンストラクターでは多項式の係数を昇べきの順で指示する。メンバーメソッド
rootsを呼び出し、すべての根を得る。
コード例を示す。
#!/usr/bin/env python
"""polynomial.py: Demonstrate equation solvers of SciPy.
"""
from numpy.polynomial.polynomial import Polynomial
# pylint: disable=invalid-name
# Define a polynomial x^3 - 2 x^2 + x - 2.
f = Polynomial([-2, 1, -2, 1])
print("polynomial: ", f)
# Solve f(x) == 0.
print("roots: ", f.roots())
私の環境での実行結果を示す。解 \(\pm i, 2\) がすべて求まっていることがわかる。
polynomial: poly([-2. 1. -2. 1.])
roots: [ -1.83880688e-16-1.j -1.83880688e-16+1.j 2.00000000e+00+0.j]
Note
本節で例示した方程式は SymPy を用いれば代数的処理で解を得られる。 SymPy 利用ノート 参照。
In [1]: solve(x - 2 * sqrt(x - 1))
Out[1]: [2]
In [2]: solve(x**3 - 2 * x**2 + x - 2)
Out[2]: [2, -I, I]
