Pythonを用いた非線形方程式の数値解法　ニュートン法 - モデリングと数値計算（主に化学工学向け）

はじめに

非線形方程式の数値解を得るためのニュートン法について、考え方とPythonを用いた計算のやり方を簡単に記す。

ニュートン法

次の非線形方程式を満たす $x$ を考える

$\displaystyle{ 2x^2 + x - 5 = 0 }$

上記方程式の左辺を $f(x)$ とおき、最初の近似解を $x_0$ としたとき、

$\displaystyle{ x_1 = x_0 - \frac{f(x_0)}{f'(x_0)}, \qquad x_2 = x_1 - \frac{f(x_1)}{f'(x_1)} }$

と、順次計算していく方法。一般化すると下記。

$\displaystyle{ x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)} }$

適当なところで計算を終了する条件の設定が必要。

Pythonで計算するなら scipy.optimize.newton を使うのが楽。

import scipy.optimize

fx = lambda x: 2*x**2 + x - 5

print(scipy.optimize.newton(fx, 2))

> 1.3507810593582126

自分で関数を作るなら、たとえば $f'(x)$ を得るのに数値微分を使って次のやり方がある。

# 数値微分(中心差分)
def dfx(fx, x):
    h = 1e-8
    return (fx(x+h) - fx(x-h))/(2*h)

def newton(fx, x0):
    x1 = x0 - fx(x0)/dfx(fx, x0)
    if abs(x0 - x1) < 1e-5:
        return x1
    else:
        return newton(fx, x1)

print(newton(fx, 2))

> 1.3507810593595635