p2_ex01

テーラ展開

モジュールのインポート

import numpy as np
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'])

$\sin x$ のテーラ展開

$\sin x$ を $x=0$ においてテイラー(Taylor)展開すると次のようになる． $$ \sin x = \sum_{n=0}^{N-1} (-1)^n \frac{x^{2n+1}}{(2n+1)!} +R_N $$ だたし，$R_N$は剰余項を示す．剰余項が十分小さくなるように展開項数を自動的に決めて計算するようにすると，

def my_sin(x): # テーラ展開によるsinの計算
    global n_ # グローバル変数
    eps = 1.0e-6 # 収束判定条件
    er = eps
    m2 =1
    s = x
    t = s
    while er >= eps:
       m1 = m2+1
       m2 = m1+1
       t = -t*(x*x/m1/m2)
       s = s+t
       er = np.abs(t)
    n_ = m1/2
    return s

フォーマット付きで計算値を出力すると，

q = 180.0/np.pi
nx = 7 # 計算点数
print(f"{'No.':>3s} {'x [deg]':>8s} {'my_sin(x)':>10s}"\
        f"{'sin(x)':>10s} {'error':>12s} {'N':>4s}")
dx = 15.0/q
for i in range(nx):
    x = i*dx
    y = my_sin(x)
    print(f"{i+1:3.0f} {x*q:8.2f} {y:10.5f}"\
            f"{np.sin(x):10.5f} {y-np.sin(x):12.8f} {n_:4.0f}")

No.  x [deg]  my_sin(x)    sin(x)        error    N
  1     0.00    0.00000   0.00000   0.00000000    1
  2    15.00    0.25882   0.25882  -0.00000000    3
  3    30.00    0.50000   0.50000   0.00000000    4
  4    45.00    0.70711   0.70711   0.00000000    4
  5    60.00    0.86603   0.86603  -0.00000000    5
  6    75.00    0.96593   0.96593  -0.00000001    5
  7    90.00    1.00000   1.00000   0.00000000    6

プロットすると，

fig = plt.figure() # グラフ領域の作成

xx2 = np.linspace(0.0, 360.0/q, 91)
plt.plot(xx2,np.sin(xx2))

nx = 6
xx = np.linspace(0.0, 90.0/q, nx)
yy = np.empty(nx)
yy[:] = [my_sin(xx[i]) for i in range(nx)]
plt.plot(xx,yy,'o')

fig.savefig('p2_ex01_1.pdf')

展開項数を与えて計算するようにし，配列計算にも対応させると，

def my_sin2(x,n): # テーラ展開によるsinの計算（展開項数n）
    m2 =1
    s = x
    t = s
    for i in range(n):
       m1 = m2+1
       m2 = m1+1
       t = -t*(x*x/m1/m2)
       s = s+t
    return s

同様にして出力すると，

q = 180.0/np.pi
nx = 24 # 計算点数
n = 12 # 展開項数
print(f"{'No.':>3s} {'x [deg]':>8s} {'my_sin(x)':>10s}"\
        f"{'sin(x)':>10s} {'error':>12s} {'N':>4s}")
dx = 15.0/q
for i in range(nx):
    x = i*dx
    y = my_sin2(x,n)
    print(f"{i+1:3.0f} {x*q:8.2f} {y:10.5f}"\
          f"{np.sin(x):10.5f} {y-np.sin(x):12.8f} {n:4.0f}")

No.  x [deg]  my_sin(x)    sin(x)        error    N
  1     0.00    0.00000   0.00000   0.00000000   12
  2    15.00    0.25882   0.25882  -0.00000000   12
  3    30.00    0.50000   0.50000   0.00000000   12
  4    45.00    0.70711   0.70711   0.00000000   12
  5    60.00    0.86603   0.86603  -0.00000000   12
  6    75.00    0.96593   0.96593  -0.00000000   12
  7    90.00    1.00000   1.00000   0.00000000   12
  8   105.00    0.96593   0.96593  -0.00000000   12
  9   120.00    0.86603   0.86603   0.00000000   12
 10   135.00    0.70711   0.70711   0.00000000   12
 11   150.00    0.50000   0.50000   0.00000000   12
 12   165.00    0.25882   0.25882   0.00000000   12
 13   180.00    0.00000   0.00000   0.00000000   12
 14   195.00   -0.25882  -0.25882   0.00000000   12
 15   210.00   -0.50000  -0.50000   0.00000000   12
 16   225.00   -0.70711  -0.70711   0.00000000   12
 17   240.00   -0.86603  -0.86603   0.00000000   12
 18   255.00   -0.96593  -0.96593   0.00000000   12
 19   270.00   -1.00000  -1.00000   0.00000000   12
 20   285.00   -0.96593  -0.96593   0.00000000   12
 21   300.00   -0.86603  -0.86603   0.00000000   12
 22   315.00   -0.70711  -0.70711   0.00000001   12
 23   330.00   -0.50000  -0.50000   0.00000003   12
 24   345.00   -0.25882  -0.25882   0.00000010   12

あるいは，リスト内包表記より出力すると，

nx2 = 6
xx = np.linspace(0.0, 90.0/q, nx2)
print(f"{'x [deg]:':>10s}",end=' ')
[print(f"{xx[i]*q:8.5f}",end=' ') for i in range(nx2)]
yy = np.empty(nx2)
yy[:] = [my_sin(xx[i]) for i in range(nx2)]
print("")
print(f"{'my_sin(x):':>10s}",end=' ')
[print(f"{yy[i]:8.5f}",end=' ') for i in range(nx2)]
print("")
print(f"{'sin(x):':>10s}",end=' ')
[print(f"{np.sin(xx[i]):8.5f}",end=' ') for i in range(nx2)]
print("")
print(f"{'error:':>10s}",end=' ')
[print(f"{yy[i]-np.sin(xx[i]):8.5f}",end=' ') for i in range(nx2)]
print("")

  x [deg]:  0.00000 18.00000 36.00000 54.00000 72.00000 90.00000 
my_sin(x):  0.00000  0.30902  0.58779  0.80902  0.95106  1.00000 
   sin(x):  0.00000  0.30902  0.58779  0.80902  0.95106  1.00000 
    error:  0.00000 -0.00000  0.00000 -0.00000 -0.00000  0.00000

展開項数$n=6,7$を比較してプロットすると，

n = 6
xx2 = np.linspace(0.0, 360.0/q, 91)

fig = plt.figure() # グラフ領域の作成
plt.plot(xx2,np.sin(xx2))
plt.plot(xx2,my_sin2(xx2,n),'.')
plt.plot(xx2,my_sin2(xx2,n+1),'.')
plt.plot(xx,yy,'o')
fig.savefig('p2_ex01_2.pdf')

sinc 関数

u = np.linspace(-6.0, 6.0, 101)

fig = plt.figure() # グラフ領域の作成
plt.xlim(-6.0, 6.0) # x軸範囲の設定
plt.ylim(-1.0, 2.0) # y軸範囲の設定
plt.xlabel(r"$x$") # x軸のラベル設定
plt.ylabel(r"$sinc(x)$") # y軸のラベル設定
plt.plot(u, np.sinc(u))# 正規化されたsinc 関数
plt.plot(u, np.sinc(u/np.pi))
plt.plot(u, np.sin(u)/u, '.')
fig.savefig('p2_ex01_3.pdf')
plt.show()