Zenn
🍵

方程式の解法をPythonを使って解いてみる【ニュートン法】

2020/12/11に公開
Python
ニュートン法
球根アルゴリズム
tech

大学で物理学を専攻しています。今回（これが初投稿だけど・・・）はニュートン法と二分法をPythonを使って解いてみようと思います。

ニュートン法とは

ニュートン法、またはニュートン・ラフソン法は、数値解析の分野いおいて、方程式を数値計算によって解くための反復法による球根アルゴリズムの一つです。

手順としては、関数f(x)について、

  1. あるx_1の数値を決定する
  2. f(x_1)を求め、そのx_1の接戦の式よりy=0と交わるxの値をx_2とする
  3. f(x_2)を求め、そのx_2の接線の式よりy=0と交わるxの値をx_3とする
  4. これを繰り返すことによって、x_1, x_2はf(x)=0となるx=aの値に近づいていく
    といった流れになっています。
    2020-03a.png

この一連の手順は、定式化されています。
接線の式はf(x_n)=f'(x_n)(x_n - x_{n+1})と表すことができるので、式変形すると、

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

といった形になります。

Pythonを用いて表現してみよう

具体例として、f(x)=e^x-4x に対して f(x)=0を解いていきます。
微分についてはf'(x)=e^x -4であることを用いります。f(x), f'(x)をそれぞれf(x), df(x)と定義し、初期値をx=4とします。

実際のコード

import numpy as np
import matplotlib.pylab as plt

def f(x):
  return np.exp(x)-4*x

def df(x):
  return np.exp(x)-4

#初期値
x=4

#epsilon1は要求精度。これよりも数値精度がよくなったら二分法の繰り返しを終了する。
epsilon1=0.00001

# n, c, f(c) を値を出力する。\t はタブという空白を適当にあける命令です。
print("n\t x\t f(x)")
n=1
# while文でf(x)=0になるまで繰り返す
while True:
  x = x - f(x)/df(x)
  print("{}\t{:.5f}\t{:.5f}".format(n, x, f(x)))
  n += 1
  if abs(f(x)) < epsilon1:
    break 

print("x= %f" % x)

すると結果はx= 2.153292という値になります。

確認してみよう

実際にPythonを用いてグラフを作成して、確認してみたいと思います。
使用するコードは以下の通りです。

import numpy as np
import matplotlib.pylab as plt

def f(x):
  return np.exp(x)-4*x

#0<x<4の間で描画
x=np.arange(0.0,4.0,0.05)
#zは、f(x)=0を描画
z=np.zeros(len(x))
plt.plot(x,f(x),x,ｚ)

これより、グラフは次のように描かれます。
ダウンロード.png
このグラフから、実際に、x= 2.153292に答えがあるとわかります。

余談

ニュートン法は、scipy.optimizeに入っているので

import scipy.optimize as spopt

def f(x):
  return np.exp(x)-4*x

spopt.newton(f,4)

とすれば、簡単に求めることもできます。scipy.optimizeについてはこちらよりご参照ください

ご高覧ありがとうございました。

参考文献

[ニュートン法]（https://mathworld.wolfram.com/NewtonsMethod.html）

Discussion