非線形連立方程式ってPythonで簡単に解けるんですね

数値計算をやっていると、非線形連立方程式を解く機会ってまぁまぁありますよね。もう何年も前ですが大学院生で研究をやっていた頃は、Fortranでコトコト実装していました。最近Pythonでプログラムを書く機会が増えてきて、再び非線形連立方程式を解く場面に遭遇しました。非線形連立方程式の解法として有名なものと言えばNewton法ですが、ゼロから実装するの面倒なのでいいライブラリがないか探しました。

「Scipy.optimize.root」で解けちゃう

半端じゃなく簡単でした。今回は例題として以下のような非線形連立方程式を解きたいと思います。

$$
\begin{align}
x^{2}+y^{2}-1=0 \\
x=0
\end{align}
$$

ソースコードは以下の通りです。あとで簡単な説明をします。

import numpy as np
from scipy import optimize

# 解きたい関数をリストで戻す
def func(x):
    return [x[0]**2 + x[1]**2 -1.0,
            x[0]]


result = optimize.root( func, [ 1.0, 0.0], method="broyden1")
print(result)

実行結果は以下の通りです。

     fun: array([-1.98898568e-07, -5.14009858e-06])
 message: 'A solution was found at the specified tolerance.'
     nit: 9
  status: 1
 success: True
       x: array([-5.14009858e-06,  9.99999901e-01])

ソースコード/実行結果の説明

scipy.optimize.rootの利用方法の詳細はscipyドキュメントを参照してください。ここではざっくりの説明をします。optimize.rootの第一引数は、解きたい関数を定義したfunc関数です。第二引数は、問題を解き始めるときに使用する初期値です。第三引数は、解き方を指定する所となります。詳細は、scipyドキュメントを参照してもらいたいのですが、ここで1つ注意点。今回はmethodでbroyden1を指定してますが、引数によってはヤコビ行列が別に定義してやる必要があります。僕は面倒だったので、ヤコビ行列が不要なものを選びました。

また、実行結果について大事なところはx：array()の部分でこれが実際の解を表しており、$x=0,y=1$っぽいのであってそう。解の詳細については、scipy.optimize.OptimizeResultを参照。

最後に

Pythonって色々簡単に実装できて便利ですね。大学院生の時に出会いたかった。いや、出会っていたんだけど乗り換えが面倒だったので見て見ぬふりをしてました。ごめんなさい。