ニュートン・ラフソン法¶

In [1]:
import numpy as np # 数値計算ライブラリ
from scipy import optimize # 非線形方程式解法など
import matplotlib.pyplot as plt # グラフ描画
import pandas as pd # 表形式のデータを扱うライブラリ
In [2]:
import scienceplots # 論文・科学技術系のプロットスタイル
plt.style.use(['science', 'notebook'])
#plt.rcParams['font.family'] = 'Times New Roman' # font familyの設定
plt.rcParams['font.family'] = 'serif' # font familyの設定
plt.rcParams['mathtext.fontset'] = 'cm' # math fontの設定
In [3]:
np.set_printoptions(precision=3) # 出力桁数の設定
In [4]:
# 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
In [5]:
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
In [6]:
def func5(x):
    y = x*x-4.0
    return y
In [7]:
def dfunc5(x):
    y = 2.0*x
    return y
In [8]:
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
In [9]:
xx = np.arange(0.0, 4.2, 0.20111)
df = pd.DataFrame({
    '$x$':xx,
    'func5':func5(xx),
    'dfunc5':dfunc5(xx),
})
#df.style.format({"$x$":"{:.3f}"})
df.style.format()
Out[9]:
  $x$ func5 dfunc5
0 0.000000 -4.000000 0.000000
1 0.201110 -3.959555 0.402220
2 0.402220 -3.838219 0.804440
3 0.603330 -3.635993 1.206660
4 0.804440 -3.352876 1.608880
5 1.005550 -2.988869 2.011100
6 1.206660 -2.543972 2.413320
7 1.407770 -2.018184 2.815540
8 1.608880 -1.411505 3.217760
9 1.809990 -0.723936 3.619980
10 2.011100 0.044523 4.022200
11 2.212210 0.893873 4.424420
12 2.413320 1.824113 4.826640
13 2.614430 2.835244 5.228860
14 2.815540 3.927265 5.631080
15 3.016650 5.100177 6.033300
16 3.217760 6.353979 6.435520
17 3.418870 7.688672 6.837740
18 3.619980 9.104255 7.239960
19 3.821090 10.600729 7.642180
20 4.022200 12.178093 8.044400
In [10]:
def func5a(x,a):
    y = a*x*x-4.0+x
    return y
In [11]:
def dfunc5a(x,a):
    y = 2.0*a*x
    return y
In [12]:
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.00000000000   2.00000000000
   1.000  -1.00000000000   4.00000000000
   1.500   2.00000000000   6.00000000000
   2.000   6.00000000000   8.00000000000
   2.500  11.00000000000  10.00000000000
   3.000  17.00000000000  12.00000000000
   3.500  24.00000000000  14.00000000000
   4.000  32.00000000000  16.00000000000
In [13]:
xx = np.linspace(0.0, 4.0, 9)
df = pd.DataFrame({
    '$x$':xx,
    'func5a':func5a(xx,a),
    'dfunc5a':dfunc5a(xx,a),
})
df.style.format({"x":"{:.3f}"})
Out[13]:
  $x$ func5a dfunc5a
0 0.000000 -4.000000 0.000000
1 0.500000 -3.000000 2.000000
2 1.000000 -1.000000 4.000000
3 1.500000 2.000000 6.000000
4 2.000000 6.000000 8.000000
5 2.500000 11.000000 10.000000
6 3.000000 17.000000 12.000000
7 3.500000 24.000000 14.000000
8 4.000000 32.000000 16.000000
In [14]:
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')
No description has been provided for this image
In [15]:
eps = 1.0e-9
x_ini = 0.000001
xzero = newton(func5, dfunc5, eps, x_ini, 50)
print ('zero of function: ',xzero)
print ('f(xzero): ',func5(xzero))

zero of function:  2.0
f(xzero):  0.0
In [16]:
optimize.newton(func5, x_ini, dfunc5)
Out[16]:
2.0
In [17]:
a=1.0
optimize.newton(func5a, x_ini, dfunc5a, args=(a,)) # 可変長引数
Out[17]:
1.561552814394198
In [18]:
optimize.newton(func5, x_ini)
Out[18]:
2.0000000000000275
In [19]:
optimize.newton(func5a, x_ini, args=(a,)) # 可変長引数
Out[19]:
1.5615528128088303
In [20]:
def func5ab(x,a,b):
    y = x*(x-a)*(x-b) # x**3-(a+b)*x*2+a*b*x
    return y
In [21]:
def dfunc5ab(x,a,b):
    y = 3.0*x**2-2*(a+b)*x+a*b
    return y
In [22]:
a, b = 1.0, 2.0
calf2_(func5ab, dfunc5ab, 0.0, 0.2, 21, a,b)

       x            f(x)           f2(x)
   0.000   0.00000000000   2.00000000000
   0.200   0.28800000000   0.92000000000
   0.400   0.38400000000   0.08000000000
   0.600   0.33600000000  -0.52000000000
   0.800   0.19200000000  -0.88000000000
   1.000  -0.00000000000  -1.00000000000
   1.200  -0.19200000000  -0.88000000000
   1.400  -0.33600000000  -0.52000000000
   1.600  -0.38400000000   0.08000000000
   1.800  -0.28800000000   0.92000000000
   2.000   0.00000000000   2.00000000000
   2.200   0.52800000000   3.32000000000
   2.400   1.34400000000   4.88000000000
   2.600   2.49600000000   6.68000000000
   2.800   4.03200000000   8.72000000000
   3.000   6.00000000000  11.00000000000
   3.200   8.44800000000  13.52000000000
   3.400  11.42400000000  16.28000000000
   3.600  14.97600000000  19.28000000000
   3.800  19.15200000000  22.52000000000
   4.000  24.00000000000  26.00000000000
In [23]:
xx = np.linspace(0.0, 4.0, 21)
df = pd.DataFrame({
    '$x$':xx,
    'func5ab':func5ab(xx,a,b),
    'dfunc5ab':dfunc5ab(xx,a,b),
})
df.style.format({"$x$":"{:.3f}"})
Out[23]:
  $x$ func5ab dfunc5ab
0 0.000 0.000000 2.000000
1 0.200 0.288000 0.920000
2 0.400 0.384000 0.080000
3 0.600 0.336000 -0.520000
4 0.800 0.192000 -0.880000
5 1.000 -0.000000 -1.000000
6 1.200 -0.192000 -0.880000
7 1.400 -0.336000 -0.520000
8 1.600 -0.384000 0.080000
9 1.800 -0.288000 0.920000
10 2.000 0.000000 2.000000
11 2.200 0.528000 3.320000
12 2.400 1.344000 4.880000
13 2.600 2.496000 6.680000
14 2.800 4.032000 8.720000
15 3.000 6.000000 11.000000
16 3.200 8.448000 13.520000
17 3.400 11.424000 16.280000
18 3.600 14.976000 19.280000
19 3.800 19.152000 22.520000
20 4.000 24.000000 26.000000
In [24]:
a_ , b_ = 1.0, 2.0
x_ini = (-0.5, 0.9, 2.5)
optimize.newton(func5ab, x_ini, args=(a_, b_)) # 可変長引数
Out[24]:
array([-3.42e-19,  1.00e+00,  2.00e+00])
In [25]:
optimize.newton(func5ab, x_ini, dfunc5ab, args=(a_, b_)) # 可変長引数
Out[25]:
array([-3.403e-22,  1.000e+00,  2.000e+00])
In [26]:
xx = np.linspace(-1.0, 2.5, 101)
fig = plt.figure() # グラフ領域の作成
plt.grid(color = "gray", linestyle="--")
plt.plot(xx,func5ab(xx,a,b), label='$f(x)$')
plt.plot(xx,dfunc5ab(xx,a,b), label=r'$\frac{df(x)}{dx}$')
plt.xlabel(r"$x$") # x軸のラベル
plt.ylabel(r"$f_1(x), f_2(x)$") # y軸のラベル
plt.legend(ncol=1, loc='best', fancybox=False, frameon=True)
fig.savefig('p2_ex05_2.pdf')
No description has been provided for this image