mono.log(e)

旅とバイクとプログラミング

非線形連立方程式の解法

勉強用 数値計算 python

非線形連立方程式の数値計算

非線形連立方程式の数値計算を行うには一般的にニュートン法が用いられるそう。 具体的には、



\begin{align}
f_1(x_1, x_2, \cdots, x_n) &= 0\\\
f_1(x_1, x_2, \cdots, x_n) &= 0\\\
&\cdots\\\
f_n(x_1, x_2, \cdots, x_n) &= 0
\end{align}

のような非線形連立方程式の解x = (x_1, x_2, \cdots, x_n)を求める。

このとき方程式を



\begin{align}
f(x) = 0
\end{align}

とまとめ、この方程式に対するニュートン法の漸化式は以下のようになる。



\begin{align}
x^{(n+1)} &= x^{(n)}+\delta x\\\
&= x^{(n)}-J(x^{(n)})^{-1}f(x^{(n)})
\end{align}

ここでJ(x)はヤコビ行列である。 したがって



\begin{align}
\delta x = -J^{-1}f(x^{(n)})
\end{align}

としてヤコビ行列の逆行列から更新後のx^{(n+1)}が求まるが、数値計算上は



\begin{align}
J\delta x = -f(x^{(n)})
\end{align}

としてガウスの消去法を用いて求解するのが通常のようである。このときの注意としては、ニュートン法の特性により、初期値は求めたい解の近傍に取る必要がある。(多分)

以上を実際にpythonを用いて実行してみるが、scipy にライブラリがあるのでこれを利用してみる。

Scipy.optimize.root

以下使ってみた結果。

gist28808d39e6a1649f14db1d28489c9e6c

参考

org-technology.com

qiita.com