ベッセル関数
モジュールのインポート
import numpy as np
from scipy.special import jv, yv # ベッセル関数
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の設定
第1種ベッセル関数の級数展開
# x: 変数
# n: ベッセル関数の次数
# m: 展開項数
def my_besselj(x, n): # 第1種ベッセル関数の級数展開
eps=1.0e-12 # 収束判定条件
a = 1.0
z = x/2.0
for i in range(n):
a = a*z/(i+1)
y = a
i = 0
while np.abs(a)>=eps:
i = i+1
a = -a*z*z/i/(i+n)
y = y+a
m = i
return y, m
第1種ベッセル関数の漸近展開
# x: 変数
# n: ベッセル関数の次数
# m: 展開項数
def besselj(x, n): # 第1種ベッセル関数の級数展開と漸近展開
eps=1.0e-12# 収束判定条件
a = 1.0
if np.abs(x)<14.0: # 級数展開
z = x/2.0
for i in range(n):
a = a*z/(i+1)
y = a
i = 0
while np.abs(a)>=eps:
i = i+1
a = -a*z*z/i/(i+n)
y = y+a
m = i
else: # 漸近展開
za = 1.0
zb = (4.0*n*n-1.0)/8.0/x
an = za
bn = zb
i = 0
i2 = 1
while np.abs(zb)>=eps:
i = i+1
i1 = i2+1
i2 = i1+1
za = -zb*(4.0*n*n-(2.0*i1-1.0)*(2.0*i1-1.0))/i1/8.0/x
zb = za*(4.0*n*n-(2.0*i2-1.0)*(2.0*i2-1.0))/i2/8.0/x
an = an+za
bn = bn+zb
m = i
th = x-(2.0*n+1.0)*np.pi/4.0
y = np.sqrt(2.0/np.pi/x)*(an*np.cos(th)-bn*np.sin(th))
return y, m
級数展開，漸近展開，および Scipy の関数 jv より計算すると，
x0, nx = 0.0, 151
xx = np.linspace(x0, 30.0, nx)
yy1 = np.empty(0)
yy2 = np.empty(0)
yy3 = np.empty(0)
print(f"{'x':>8s} {'m1':>4s} {'my_besselj':>15s} {'m2':>4s}"\
f"{'besselj':>15s} {'scipy.special.jv':>15s}")
for i in range(nx):
y1, m1 = my_besselj(xx[i], 0)
y2, m2 = besselj(xx[i], 0)
y3 = jv(0,xx[i]) # 第1種ベッセル関数
print(f"{xx[i]:8.2f} {m1:4.0f} {y1:15.11f} {m2:4.0f}"\
f"{y2:15.11f} {y3:15.11f}")
yy1 = np.append(yy1, y1)
yy2 = np.append(yy2, y2)
yy3 = np.append(yy3, y3)
print("")
x m1 my_besselj m2 besselj scipy.special.jv
0.00 1 1.00000000000 1 1.00000000000 1.00000000000
0.20 5 0.99002497224 5 0.99002497224 0.99002497224
0.40 6 0.96039822666 6 0.96039822666 0.96039822666
0.60 7 0.91200486350 7 0.91200486350 0.91200486350
0.80 7 0.84628735275 7 0.84628735275 0.84628735275
1.00 8 0.76519768656 8 0.76519768656 0.76519768656
1.20 8 0.67113274426 8 0.67113274426 0.67113274426
1.40 9 0.56685512037 9 0.56685512037 0.56685512037
1.60 9 0.45540216764 9 0.45540216764 0.45540216764
1.80 10 0.33998641104 10 0.33998641104 0.33998641104
2.00 10 0.22389077914 10 0.22389077914 0.22389077914
2.20 10 0.11036226692 10 0.11036226692 0.11036226692
2.40 11 0.00250768330 11 0.00250768330 0.00250768330
2.60 11 -0.09680495440 11 -0.09680495440 -0.09680495440
2.80 12 -0.18503603336 12 -0.18503603336 -0.18503603336
3.00 12 -0.26005195490 12 -0.26005195490 -0.26005195490
3.20 12 -0.32018816966 12 -0.32018816966 -0.32018816966
3.40 13 -0.36429559676 13 -0.36429559676 -0.36429559676
3.60 13 -0.39176898370 13 -0.39176898370 -0.39176898370
3.80 13 -0.40255641018 13 -0.40255641018 -0.40255641018
4.00 14 -0.39714980986 14 -0.39714980986 -0.39714980986
4.20 14 -0.37655705437 14 -0.37655705437 -0.37655705437
4.40 14 -0.34225679000 14 -0.34225679000 -0.34225679000
4.60 15 -0.29613781657 15 -0.29613781657 -0.29613781657
4.80 15 -0.24042532729 15 -0.24042532729 -0.24042532729
5.00 15 -0.17759677131 15 -0.17759677131 -0.17759677131
5.20 16 -0.11029043979 16 -0.11029043979 -0.11029043979
5.40 16 -0.04121010124 16 -0.04121010124 -0.04121010124
5.60 16 0.02697088469 16 0.02697088469 0.02697088469
5.80 17 0.09170256757 17 0.09170256757 0.09170256757
6.00 17 0.15064525725 17 0.15064525725 0.15064525725
6.20 17 0.20174722295 17 0.20174722295 0.20174722295
6.40 18 0.24331060482 18 0.24331060482 0.24331060482
6.60 18 0.27404336062 18 0.27404336062 0.27404336062
6.80 18 0.29309560310 18 0.29309560310 0.29309560310
7.00 18 0.30007927052 18 0.30007927052 0.30007927052
7.20 19 0.29507069140 19 0.29507069140 0.29507069140
7.40 19 0.27859623266 19 0.27859623266 0.27859623266
7.60 19 0.25160183385 19 0.25160183385 0.25160183385
7.80 20 0.21540780775 20 0.21540780775 0.21540780775
8.00 20 0.17165080714 20 0.17165080714 0.17165080714
8.20 20 0.12221530178 20 0.12221530178 0.12221530178
8.40 21 0.06915726166 21 0.06915726166 0.06915726166
8.60 21 0.01462299128 21 0.01462299128 0.01462299128
8.80 21 -0.03923380318 21 -0.03923380318 -0.03923380318
9.00 22 -0.09033361118 22 -0.09033361118 -0.09033361118
9.20 22 -0.13674837076 22 -0.13674837076 -0.13674837076
9.40 22 -0.17677157275 22 -0.17677157275 -0.17677157275
9.60 22 -0.20897871837 22 -0.20897871837 -0.20897871837
9.80 23 -0.23227602758 23 -0.23227602758 -0.23227602758
10.00 23 -0.24593576445 23 -0.24593576445 -0.24593576445
10.20 23 -0.24961706985 23 -0.24961706985 -0.24961706985
10.40 24 -0.24337175071 24 -0.24337175071 -0.24337175071
10.60 24 -0.22763504762 24 -0.22763504762 -0.22763504762
10.80 24 -0.20320196711 24 -0.20320196711 -0.20320196711
11.00 24 -0.17119030041 24 -0.17119030041 -0.17119030041
11.20 25 -0.13299193686 25 -0.13299193686 -0.13299193686
11.40 25 -0.09021450025 25 -0.09021450025 -0.09021450025
11.60 25 -0.04461567409 25 -0.04461567409 -0.04461567409
11.80 26 0.00196717331 26 0.00196717331 0.00196717331
12.00 26 0.04768931080 26 0.04768931080 0.04768931080
12.20 26 0.09077012317 26 0.09077012317 0.09077012317
12.40 26 0.12956102652 26 0.12956102652 0.12956102652
12.60 27 0.16260727175 27 0.16260727175 0.16260727175
12.80 27 0.18870135478 27 0.18870135478 0.18870135478
13.00 27 0.20692610238 27 0.20692610238 0.20692610238
13.20 28 0.21668592226 28 0.21668592226 0.21668592226
13.40 28 0.21772517873 28 0.21772517873 0.21772517873
13.60 28 0.21013316137 28 0.21013316137 0.21013316137
13.80 29 0.19433563522 29 0.19433563522 0.19433563522
14.00 29 0.17107347611 9 0.17107347611 0.17107347611
14.20 29 0.14136938466 9 0.14136938466 0.14136938466
14.40 29 0.10648411849 9 0.10648411849 0.10648411849
14.60 30 0.06786406832 8 0.06786406832 0.06786406832
14.80 30 0.02708231458 8 0.02708231459 0.02708231459
15.00 30 -0.01422447282 8 -0.01422447283 -0.01422447283
15.20 31 -0.05442079685 8 -0.05442079684 -0.05442079684
15.40 31 -0.09193622786 8 -0.09193622786 -0.09193622786
15.60 31 -0.12532596403 8 -0.12532596402 -0.12532596402
15.80 31 -0.15332574775 7 -0.15332574776 -0.15332574776
16.00 32 -0.17489907397 7 -0.17489907398 -0.17489907398
16.20 32 -0.18927494696 7 -0.18927494698 -0.18927494698
16.40 32 -0.19597482880 7 -0.19597482879 -0.19597482879
16.60 32 -0.19482785581 7 -0.19482785581 -0.19482785581
16.80 33 -0.18597386530 7 -0.18597386535 -0.18597386535
17.00 33 -0.16985425216 7 -0.16985425215 -0.16985425215
17.20 33 -0.14719114667 7 -0.14719114677 -0.14719114677
17.40 34 -0.11895585622 7 -0.11895585634 -0.11895585634
17.60 34 -0.08632791551 7 -0.08632791550 -0.08632791550
17.80 34 -0.05064644612 7 -0.05064644607 -0.05064644607
18.00 34 -0.01335580590 6 -0.01335580572 -0.01335580572
18.20 35 0.02405229239 6 0.02405229241 0.02405229241
18.40 35 0.06009789197 6 0.06009789204 0.06009789204
18.60 35 0.09337081654 6 0.09337081628 0.09337081628
18.80 36 0.12258534424 6 0.12258534436 0.12258534436
19.00 36 0.14662943987 6 0.14662943966 0.14662943966
19.20 36 0.16460665912 6 0.16460665908 0.16460665908
19.40 36 0.17586917782 6 0.17586917809 0.17586917809
19.60 37 0.18004072735 6 0.18004072743 0.18004072743
19.80 37 0.17702864303 6 0.17702864315 0.17702864315
20.00 37 0.16702466460 6 0.16702466434 0.16702466434
20.20 38 0.15049455633 6 0.15049455643 0.15049455643
20.40 38 0.12815707287 6 0.12815707345 0.12815707345
20.60 38 0.10095318643 6 0.10095318521 0.10095318521
20.80 38 0.07000686787 6 0.07000686745 0.07000686745
21.00 39 0.03657906881 6 0.03657907100 0.03657907100
21.20 39 0.00201673865 6 0.00201673882 0.00201673882
21.40 39 -0.03230108312 6 -0.03230108377 -0.03230108377
21.60 39 -0.06501790346 6 -0.06501790416 -0.06501790416
21.80 40 -0.09485305425 6 -0.09485305082 -0.09485305082
22.00 40 -0.12065146429 5 -0.12065147570 -0.12065147570
22.20 40 -0.14142816201 5 -0.14142816591 -0.14142816591
22.40 41 -0.15640546436 5 -0.15640546714 -0.15640546714
22.60 41 -0.16504192246 5 -0.16504191464 -0.16504191464
22.80 41 -0.16705150060 5 -0.16705151191 -0.16705151191
23.00 41 -0.16241277065 5 -0.16241278131 -0.16241278131
23.20 42 -0.15136732526 5 -0.15136731788 -0.15136731788
23.40 42 -0.13440799715 5 -0.13440799160 -0.13440799160
23.60 42 -0.11225733834 5 -0.11225734890 -0.11225734890
23.80 43 -0.08583717239 5 -0.08583714431 -0.08583714431
24.00 43 -0.05623024426 5 -0.05623027417 -0.05623027417
24.20 43 -0.02463662996 5 -0.02463667239 -0.02463667239
24.40 43 0.00767509218 5 0.00767504678 0.00767504678
24.60 44 0.03941817704 5 0.03941826119 0.03941826119
24.80 44 0.06933932507 5 0.06933931475 0.06933931475
25.00 44 0.09626663715 5 0.09626678328 0.09626678328
25.20 44 0.11915710342 5 0.11915710520 0.11915710520
25.40 45 0.13713488057 5 0.13713480805 0.13713480805
25.60 45 0.14952557947 5 0.14952578764 0.14952578764
25.80 45 0.15588249714 5 0.15588238219 0.15588238219
26.00 46 0.15599937257 5 0.15599931552 0.15599931552
26.20 46 0.14992010900 5 0.14991995044 0.14991995044
26.40 46 0.13793277681 5 0.13793267833 0.13793267833
26.60 46 0.12055775542 5 0.12055766118 0.12055766118
26.80 47 0.09852459735 5 0.09852452025 0.09852452025
27.00 47 0.07274168414 5 0.07274191801 0.07274191801
27.20 47 0.04426043719 5 0.04426029224 0.04426029224
27.40 48 0.01422857645 5 0.01422926235 0.01422926235
27.60 48 -0.01614587597 5 -0.01614857358 -0.01614857358
27.80 48 -0.04566615084 5 -0.04566560175 -0.04566560175
28.00 48 -0.07315711454 5 -0.07315701055 -0.07315701055
28.20 49 -0.09754530343 5 -0.09754657704 -0.09754657704
28.40 49 -0.11789004546 5 -0.11788862895 -0.11788862895
28.60 49 -0.13340153204 5 -0.13340455231 -0.13340455231
28.80 49 -0.14351250512 5 -0.14351244102 -0.14351244102
29.00 50 -0.14784115220 5 -0.14784876468 -0.14784876468
29.20 50 -0.14628263202 5 -0.14628125361 -0.14628125361
29.40 50 -0.13890540622 5 -0.13891255151 -0.13891255151
29.60 51 -0.12606301096 5 -0.12607455379 -0.12607455379
29.80 51 -0.10831830045 5 -0.10831371758 -0.10831371758
30.00 51 -0.08636246307 5 -0.08636798358 -0.08636798358
これらをプロットすると，
fig = plt.figure() # グラフ領域の作成
plt.grid(color = "gray", linestyle="--")
plt.minorticks_on() # #補助目盛りをつける
plt.xlabel(r"$x$") # x軸のラベル設定
plt.ylabel(r"Bessel function of the first kind $J_0(x)$") # y軸のラベル設定
plt.plot(xx, yy1, 'o', label="series expansion")
plt.plot(xx, yy2, '.', label="series/asymptotic expansions")
plt.plot(xx, yy3, label="scipy.special.jv")
plt.legend(ncol=1, loc='upper right', fancybox=False, frameon = True)
fig.savefig('p2_ex06_1.pdf')
fig = plt.figure() # グラフ領域の作成
plt.grid(color = "gray", linestyle="--")
plt.minorticks_on() # #補助目盛りをつける
plt.xlabel(r"$x$") # x軸のラベル設定
plt.ylabel(r"Bessel function of the first kind $J_0(x)$") # y軸のラベル設定
plt.plot(xx[0:nx//3-1], yy1[0:nx//3-1], label="series expansion")
plt.plot(xx[nx//3:2*nx//3-1], yy2[nx//3:2*nx//3-1], label="series/asymptotic expansions")
plt.plot(xx[2*nx//3:nx], yy3[2*nx//3:nx], label="scipy.special.jv")
plt.legend(ncol=1, loc='upper right', fancybox=False, frameon = True)
fig.savefig('p2_ex06_2.pdf')
0,1,2次の第1種ベッセル関数
fig = plt.figure() # グラフ領域の作成
plt.grid(color = "gray", linestyle="--")
plt.minorticks_on() # #補助目盛りをつける
plt.xlabel(r"$x$") # x軸のラベル設定
plt.ylabel(r"Bessel functions") # y軸のラベル設定
plt.ylim(-1.0, 2.0) # y軸範囲の設定
plt.plot(xx, jv(0,xx), label=r"$J_0(x)$")
plt.plot(xx, jv(1,xx), label=r"$J_1(x)$")
plt.plot(xx, jv(2,xx), label=r"$J_2(x)$")
plt.legend(ncol=1, loc='upper right', fancybox=False, frameon = True)
fig.savefig('p2_ex06_3.pdf')
0,1,2次の第2種ベッセル関数
fig = plt.figure() # グラフ領域の作成
plt.grid(color = "gray", linestyle="--")
plt.minorticks_on() # #補助目盛りをつける
plt.xlabel(r"$x$") # x軸のラベル設定
plt.ylabel(r"Bessel functions") # y軸のラベル設定
plt.ylim(-2.0, 1.0) # y軸範囲の設定
plt.plot(xx, yv(0,xx), label=r"$Y_0(x)$")
plt.plot(xx, yv(1,xx), label=r"$Y_1(x)$")
plt.plot(xx, yv(2,xx), label=r"$Y_2(x)$")
plt.legend(ncol=1, loc='lower right', fancybox=False, frameon = True)
fig.savefig('p2_ex06_4.pdf')
前のページに戻る
「目次」のページに戻る