ニュートン法

モジュールのインポート

import numpy as np
from scipy import optimize
import matplotlib.pyplot as plt

import scienceplots
#Warning : As of version 2.0.0, you need to add import scienceplots before setting the style (plt.style.use('science')).
plt.style.use(['science', 'notebook'])

# rcPramsを使ってTeXのフォントをデフォルトにする
#plt.rcParams['font.family'] = 'Times New Roman' # font familyの設定
plt.rcParams['font.family'] = 'serif' # font familyの設定
plt.rcParams['mathtext.fontset'] = 'cm' # math fontの設定

ニュートン法の関数

# eps: 収束判定条件
# x00: 解の初期値
# max: 繰り返し回数の上限
def newton(func, dfunc, eps, x00, nmax): # ニュートン・ラフソン法
    global nx
    icon = -1
    x0 = x00
    for i in range(nmax):
       x = x0-func(x0)/dfunc(x0)
       if np.abs(x-x0)>eps:
           x0 = x
       else:
           icon = 0
           break
    xz = x
    nx = i
    if icon==-1:# Error in newton
       raise Exception('Error in newton')
    return xz

関数の計算

def calf2(func, func2, x0, dx, nx):# 与えられた関数計算，出力
    print(f"{'x':>8s} {'f(x)':>15s} {'f2(x)':>15s}")
    for i in range(nx):
        x = x0 + dx*i
        y = func(x)
        y2 = func2(x)
        print(f"{x:>8.3f} {y:>15.11f} {y2:>15.11f}")
    return

関数と導関数

def func5(x):
    y = x*x-4.0
    return y

def dfunc5(x):
    y = 2.0*x
    return y

calf2(func5, dfunc5, 0.0, 0.2, 21)

x            f(x)           f2(x)
0.000  -4.00000000000   0.00000000000
0.200  -3.96000000000   0.40000000000
0.400  -3.84000000000   0.80000000000
0.600  -3.64000000000   1.20000000000
0.800  -3.36000000000   1.60000000000
1.000  -3.00000000000   2.00000000000
1.200  -2.56000000000   2.40000000000
1.400  -2.04000000000   2.80000000000
1.600  -1.44000000000   3.20000000000
1.800  -0.76000000000   3.60000000000
2.000   0.00000000000   4.00000000000
2.200   0.84000000000   4.40000000000
2.400   1.76000000000   4.80000000000
2.600   2.76000000000   5.20000000000
2.800   3.84000000000   5.60000000000
3.000   5.00000000000   6.00000000000
3.200   6.24000000000   6.40000000000
3.400   7.56000000000   6.80000000000
3.600   8.96000000000   7.20000000000
3.800  10.44000000000   7.60000000000
4.000  12.00000000000   8.00000000000

def calf2_(func, func2, x0, dx, nx, *args):# 与えられた関数計算，出力
    print(f"{'x':>8s} {'f(x)':>15s} {'f2(x)':>15s}")
    for i in range(nx):
        x = x0 + dx*i
        y = func(x,*args)
        y2 = func2(x,*args)
        print(f"{x:>8.3f} {y:>15.11f} {y2:>15.11f}")
    return

def func5a(x,a):
    y = a*x*x-4.0
    return y

def dfunc5a(x,a):
    y = 2.0*a*x
    return y

a=2.0
calf2_(func5a, dfunc5a, 0.0, 0.5, 9, a)

       x            f(x)           f2(x)
   0.000  -4.00000000000   0.00000000000
   0.500  -3.50000000000   2.00000000000
   1.000  -2.00000000000   4.00000000000
   1.500   0.50000000000   6.00000000000
   2.000   4.00000000000   8.00000000000
   2.500   8.50000000000  10.00000000000
   3.000  14.00000000000  12.00000000000
   3.500  20.50000000000  14.00000000000
   4.000  28.00000000000  16.00000000000

xx = np.linspace(-4.0, 4.0, 101)
fig = plt.figure() # グラフ領域の作成
plt.grid(color = "gray", linestyle="--")
plt.plot(xx,func5(xx), label='$f(x)$')
plt.plot(xx,dfunc5(xx), label=r'$\frac{df(x)}{dx}$')
plt.xlabel(r"$x$") # x軸のラベル
plt.ylabel(r"$f(x), f'(x)$") # y軸のラベル
plt.legend(ncol=1, loc='lower right', fancybox=False, frameon=True)
fig.savefig('p2_ex05_1.pdf')

eps = 1.0e-9
x_ini = 1.0
xzero = newton(func5, dfunc5, eps, x_ini, 50)
print ('zero of function: ',xzero,nx)
print ('f(xzero): ',func5(xzero))

zero of function:  2.0 5
f(xzero):  0.0

optimize.newton(func5, x_ini, dfunc5)

2.0

a=1.0
optimize.newton(func5a, x_ini, dfunc5a, args=(a,)) # 可変長引数

2.0

optimize.newton(func5, x_ini)

2.0

optimize.newton(func5a, x_ini, args=(a,)) # 可変長引数

2.0

xzero2 = optimize.newton(func5, x_ini, dfunc5)
print ('optimize.newton: ',xzero2)

optimize.newton:  2.0

xzero3 = optimize.newton(func5, x_ini)
print ('optimize.newton: ',xzero3)

optimize.newton:  2.0000000000000004