p1_ex04

関数の零点

モジュールのインポート

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

$f_1(x) = 0$ の零点

関数 $f_1(x) = -(x+1.5)x(x-1.5) = 0$ の零点を数値的に求めよう（因数分解した形で与えているので，零点はわかっている）．

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

零点を求める前に，与えらえた関数の値を計算し，プロットしておこう．関数の値を計算し，出力するためのユーザ関数は，

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

これより，計算，出力すると，

eps = 1.0e-15
x1, x2 = 0.1, 4.0
calf(func4a, 0.0, 0.2, 21)

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, 21)
yy = func4a(xx)

fig = plt.figure() # グラフ領域の作成
plt.plot(xx,yy)
fig.savefig('p1_ex04_1.pdf')

零点を，
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

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) = 0$ の零点

関数 $f_2(x)= \frac{\sin x}{x} + \cos x = 0$ の零点を数値的に求めよう．

def func4b(x):
    if x==0.0:
        t = 1.0
    else:
        t = np.sin(x)/x
    y = t + 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

関数 $f_2(x)$ は，次のようなコードでも計算すると使いやすい．

def func4bb(x):
    t = np.where(x==0.0, 1, np.sin(x)/x)
    y = t + np.cos(x)
    return y

プロットすると，

xx = np.linspace(-2.0, 2.0, 21)
yy = func4bb(xx)

fig = plt.figure() # グラフ領域の作成
plt.plot(xx,yy)
fig.savefig('p1_ex04_2.pdf')

零点は，

xzero2 = optimize.bisect(func4b, x1, x2)
print ('optimize.bisect : ',xzero2)
xzero3 = optimize.brentq(func4b, x1, x2)
print ('optimize.brentq : ',xzero3)

optimize.bisect :  2.0287578381104594
optimize.brentq :  2.0287578381106566

【問題】$f_3(x)=\frac{\sin x}{x} + \cos x$の概形を描き，零点を求めよ．