p2_ex04

関数の零点

モジュールのインポート

import numpy as np
from scipy import optimize
import matplotlib.pyplot as plt
from scipy.signal import find_peaks # ピーク値を探す

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'])

関数の計算

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

関数の零点：$f_1(x) = -(x+1.5)x(x-1.5) =0$

関数を$y = -(x+1.5)x(x-1.5)$とすると，

def func4a(x):
    y = -(x+1.5)*x*(x-1.5) 
    return y

eps = 1.0e-15
# f1(x)=0
x1, x2 = 0.1, 4.0
print("f(x)=-(x+1.5)*x*(x-1.5)")
calf(func4a, 0.0, 0.2, 21)

f(x)=-(x+1.5)*x*(x-1.5)
       x            f(x)
   0.000   0.00000000000
   0.200   0.44200000000
   0.400   0.83600000000
   0.600   1.13400000000
   0.800   1.28800000000
   1.000   1.25000000000
   1.200   0.97200000000
   1.400   0.40600000000
   1.600  -0.49600000000
   1.800  -1.78200000000
   2.000  -3.50000000000
   2.200  -5.69800000000
   2.400  -8.42400000000
   2.600 -11.72600000000
   2.800 -15.65200000000
   3.000 -20.25000000000
   3.200 -25.56800000000
   3.400 -31.65400000000
   3.600 -38.55600000000
   3.800 -46.32200000000
   4.000 -55.00000000000

フリーフォーマットでは，

xx = np.linspace(-2.0, 2.0, 51)
yy = func4a(xx)
xx,yy

    (array([-2.  , -1.92, -1.84, -1.76, -1.68, -1.6 , -1.52, -1.44, -1.36,
    -1.28, -1.2 , -1.12, -1.04, -0.96, -0.88, -0.8 , -0.72, -0.64,
    -0.56, -0.48, -0.4 , -0.32, -0.24, -0.16, -0.08,  0.  ,  0.08,
     0.16,  0.24,  0.32,  0.4 ,  0.48,  0.56,  0.64,  0.72,  0.8 ,
     0.88,  0.96,  1.04,  1.12,  1.2 ,  1.28,  1.36,  1.44,  1.52,
     1.6 ,  1.68,  1.76,  1.84,  1.92,  2.  ]),
array([ 3.5     ,  2.757888,  2.089504,  1.491776,  0.961632,  0.496   ,
     0.091808, -0.254016, -0.544544, -0.782848, -0.972   , -1.115072,
    -1.215136, -1.275264, -1.298528, -1.288   , -1.246752, -1.177856,
    -1.084384, -0.969408, -0.836   , -0.687232, -0.526176, -0.355904,
    -0.179488,  0.      ,  0.179488,  0.355904,  0.526176,  0.687232,
     0.836   ,  0.969408,  1.084384,  1.177856,  1.246752,  1.288   ,
     1.298528,  1.275264,  1.215136,  1.115072,  0.972   ,  0.782848,
     0.544544,  0.254016, -0.091808, -0.496   , -0.961632, -1.491776,
    -2.089504, -2.757888, -3.5     ]))

ピーク値は，

locs, _ = find_peaks(yy) # ピーク値
xx[locs],yy[locs]

(array([0.88]), array([1.298528]))

関数の概形をプロットすると，

fig = plt.figure() # グラフ領域の作成
plt.plot(xx,yy)
plt.plot([-10,10],[0,0],'--',color='gray',alpha=0.5)
plt.plot(xx[locs],yy[locs],'o')
plt.xlim(-2.0, 2.0)
#plt.ylim(-4.0, 4.0)
plt.xlabel(r"$x$")
plt.ylabel(r"$f(x)$")
fig.savefig('p2_ex04_1.pdf')

二分法

def bisect(func, eps, x01, x02): # 二分法
    x1 = x01
    x2 = x02
    y = func(x1)
    while np.abs(x2-x1)>eps:
        xm = (x1+x2)*0.5
        if y*func(xm)>0.0:
            x1 = xm
        else:
            x2 = xm
    return xm

区間$[1,4]$での関数の零点を二分法より求めると，

x1, x2 = 0.1, 4.0
xzero = bisect(func4a, eps, x1, x2)
print ('zero of function: ',xzero)

zero of function:  1.4999999999999991

Scipyの二分法より求めると，
https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.bisect.html#scipy.optimize.bisect

xzero2 = optimize.bisect(func4a, x1, x2)
print ('optimize.bisect : ',xzero2)

optimize.bisect :  1.4999999999991358

Scipyのブレント法より求めると，
https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.brentq.html

xzero3 = optimize.brentq(func4a, x1, x2)
print ('optimize.brentq : ',xzero3,"\n")

optimize.brentq :  1.5

関数の零点：$f_2(x) = \frac{\sin x}{x} + \cos x = 0$

def func4b(x):
    u = x==0
    t = np.where(u, 1.0, np.sin(x)/x)
    y = t + np.cos(x)
    return y

def func4c(x):
    y = np.sinc(x/np.pi) + np.cos(x)
    return y

関数の値は，

print("f(x)=sin(x)/x + cos(x)")
calf(func4b, -2.0, 0.2, 21)

f(x)=sin(x)/x + cos(x)
       x            f(x)
  -2.000   0.03850187687
  -1.800   0.31382436691
  -1.600   0.59553397960
  -1.400   0.87385980718
  -1.200   1.13905699278
  -1.000   1.38177329068
  -0.800   1.59340182297
  -0.600   1.76640640390
  -0.400   1.89460684977
  -0.200   1.97341323182
   0.000   2.00000000000
   0.200   1.97341323182
   0.400   1.89460684977
   0.600   1.76640640390
   0.800   1.59340182297
   1.000   1.38177329068
   1.200   1.13905699278
   1.400   0.87385980718
   1.600   0.59553397960
   1.800   0.31382436691
   2.000   0.03850187687

プロットすると，

xx2 = np.linspace(-4.0, 4.0, 81)
yy2 = func4c(xx2)

fig = plt.figure() # グラフ領域の作成
plt.plot(xx2,yy2)
fig.savefig('p2_ex04_2.pdf')

関数の零点は，

x1, x2 = 0.1, 3.0
xzero = bisect(func4b, eps, x1, x2)
print ('zero of function: ',xzero)
xzero2 = optimize.bisect(func4b, x1, x2)
print ('optimize.bisect : ',xzero2)
xzero3 = optimize.brentq(func4b, x1, x2)
print ('optimize.brentq : ',xzero3)

zero of function:  2.028757838110434
optimize.bisect :  2.028757838110187
optimize.brentq :  2.0287578381103835