ベッセル関数
モジュールのインポート
import numpy as np
from scipy.special import jv, yv # ベッセル関数
from scipy.special import jn_zeros, jnp_zeros, jnjnp_zeros, jnyn_zeros # ベッセル関数の零点
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'])
第1種ベッセル関数
関数の値を計算すると，
x0, nx = 0.0, 151
xx = np.linspace(x0, 30.0, nx)
yy3 = np.empty(0)
for i in range(nx):
y3 = jv(0,xx[i]) # 第1種ベッセル関数
print(f"{xx[i]:8.2f} {y3:15.11f}")
yy3 = np.append(yy3, y3)
0.00 1.00000000000
0.20 0.99002497224
0.40 0.96039822666
0.60 0.91200486350
0.80 0.84628735275
1.00 0.76519768656
1.20 0.67113274426
1.40 0.56685512037
1.60 0.45540216764
1.80 0.33998641104
2.00 0.22389077914
2.20 0.11036226692
2.40 0.00250768330
2.60 -0.09680495440
2.80 -0.18503603336
3.00 -0.26005195490
3.20 -0.32018816966
3.40 -0.36429559676
3.60 -0.39176898370
3.80 -0.40255641018
4.00 -0.39714980986
4.20 -0.37655705437
4.40 -0.34225679000
4.60 -0.29613781657
4.80 -0.24042532729
5.00 -0.17759677131
5.20 -0.11029043979
5.40 -0.04121010124
5.60 0.02697088469
5.80 0.09170256757
6.00 0.15064525725
6.20 0.20174722295
6.40 0.24331060482
6.60 0.27404336062
6.80 0.29309560310
7.00 0.30007927052
7.20 0.29507069140
7.40 0.27859623266
7.60 0.25160183385
7.80 0.21540780775
8.00 0.17165080714
8.20 0.12221530178
8.40 0.06915726166
8.60 0.01462299128
8.80 -0.03923380318
9.00 -0.09033361118
9.20 -0.13674837076
9.40 -0.17677157275
9.60 -0.20897871837
9.80 -0.23227602758
10.00 -0.24593576445
10.20 -0.24961706985
10.40 -0.24337175071
10.60 -0.22763504762
10.80 -0.20320196711
11.00 -0.17119030041
11.20 -0.13299193686
11.40 -0.09021450025
11.60 -0.04461567409
11.80 0.00196717331
12.00 0.04768931080
12.20 0.09077012317
12.40 0.12956102652
12.60 0.16260727175
12.80 0.18870135478
13.00 0.20692610238
13.20 0.21668592226
13.40 0.21772517873
13.60 0.21013316137
13.80 0.19433563522
14.00 0.17107347611
14.20 0.14136938466
14.40 0.10648411849
14.60 0.06786406832
14.80 0.02708231459
15.00 -0.01422447283
15.20 -0.05442079684
15.40 -0.09193622786
15.60 -0.12532596402
15.80 -0.15332574776
16.00 -0.17489907398
16.20 -0.18927494698
16.40 -0.19597482879
16.60 -0.19482785581
16.80 -0.18597386535
17.00 -0.16985425215
17.20 -0.14719114677
17.40 -0.11895585634
17.60 -0.08632791550
17.80 -0.05064644607
18.00 -0.01335580572
18.20 0.02405229241
18.40 0.06009789204
18.60 0.09337081628
18.80 0.12258534436
19.00 0.14662943966
19.20 0.16460665908
19.40 0.17586917809
19.60 0.18004072743
19.80 0.17702864315
20.00 0.16702466434
20.20 0.15049455643
20.40 0.12815707345
20.60 0.10095318521
20.80 0.07000686745
21.00 0.03657907100
21.20 0.00201673882
21.40 -0.03230108377
21.60 -0.06501790416
21.80 -0.09485305082
22.00 -0.12065147570
22.20 -0.14142816591
22.40 -0.15640546714
22.60 -0.16504191464
22.80 -0.16705151191
23.00 -0.16241278131
23.20 -0.15136731788
23.40 -0.13440799160
23.60 -0.11225734890
23.80 -0.08583714431
24.00 -0.05623027417
24.20 -0.02463667239
24.40 0.00767504678
24.60 0.03941826119
24.80 0.06933931475
25.00 0.09626678328
25.20 0.11915710520
25.40 0.13713480805
25.60 0.14952578764
25.80 0.15588238219
26.00 0.15599931552
26.20 0.14991995044
26.40 0.13793267833
26.60 0.12055766118
26.80 0.09852452025
27.00 0.07274191801
27.20 0.04426029224
27.40 0.01422926235
27.60 -0.01614857358
27.80 -0.04566560175
28.00 -0.07315701055
28.20 -0.09754657704
28.40 -0.11788862895
28.60 -0.13340455231
28.80 -0.14351244102
29.00 -0.14784876468
29.20 -0.14628125361
29.40 -0.13891255151
29.60 -0.12607455379
29.80 -0.10831371758
30.00 -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, yy3, label="scipy.special.jv")
plt.legend(ncol=1, loc='upper right', fancybox=False, frameon=True)
fig.savefig('p1_ex06_1.pdf')
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('p1_ex06_2.pdf')
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('p1_ex06_3.pdf')
ベッセル関数の導関数
def func_dbesj(x): # 第1種ベッセル関数の導関数
f = (jv(m-1,x)-jv(m+1,x))*0.5
return f
これより，ベッセル関数とその導関数の値を求めると，
x_min, x_max, n_x = 0.0, 3.0, 11
x = np.linspace(x_min, x_max, n_x)
m = 1 #ベッセル関数の次数
y1 = jv(m,x) # 第1種ベッセル関数 for TM1m-mode
y2 = func_dbesj(x) # for TE1m-mode
出力用のユーザ関数
def table2(sx,sy1,sy2,x,y1,y2):
nx = len(x)
print(f"{sx:>7s}",f"{sy1:>7s}",f"{sy2:>7s}")
for i in range(nx):
print(f'{x[i]:7.2f}',f'{y1[i]:7.2f}',f'{y2[i]:7.2f}')
return
を用いて，
table2('x','J_m(x)','dJ_m(x)/dx', x, y1, y2)
x J_m(x) dJ_m(x)/dx
0.00 0.00 0.50
0.30 0.15 0.48
0.60 0.29 0.43
0.90 0.41 0.36
1.20 0.50 0.26
1.50 0.56 0.14
1.80 0.58 0.02
2.10 0.57 -0.10
2.40 0.52 -0.21
2.70 0.44 -0.31
3.00 0.34 -0.37
これをプロットすると，
fig = plt.figure() # グラフ領域の作成
plt.plot(x,y1)
plt.plot(x,y2)
fig.savefig('p1_ex06_4.pdf')
x_min, x_max, n_x = 0.0, 7.0, 101
x = np.linspace(x_min, x_max, n_x)
m = 1 #ベッセル関数の次数
fig = plt.figure() # グラフ領域の作成
plt.plot(x,jv(m,x))
plt.plot(x,func_dbesj(x))
fig.savefig('p1_ex06_5.pdf')
整数次 $n$ の第1種ベッセル関数の零点：$J_n(x)=0$
https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.jn_zeros.html
nt = 6
n_max = 3
for n in range(n_max+1):
print(n, jn_zeros(n ,nt))
0 [ 2.40482556 5.52007811 8.65372791 11.79153444 14.93091771 18.07106397]
1 [ 3.83170597 7.01558667 10.17346814 13.32369194 16.47063005 19.61585851]
2 [ 5.1356223 8.41724414 11.61984117 14.79595178 17.95981949 21.11699705]
3 [ 6.3801619 9.76102313 13.01520072 16.22346616 19.40941523 22.58272959]
jn_zeros(0,15)
array([ 2.40482556, 5.52007811, 8.65372791, 11.79153444, 14.93091771,
18.07106397, 21.21163663, 24.35247153, 27.49347913, 30.63460647,
33.77582021, 36.91709835, 40.05842576, 43.19979171, 46.34118837])
整数次（$n\le 1200$）までの第1種ベッセル関数および導関数の零点：$J_n(x)=0$，$J_n'(x)=0$
https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.jnjnp_zeros.html
nt = 50
zo, n_, m_, t_ = jnjnp_zeros(nt)
n_ # 次数
array([ 0, 1, 0, 2, 1, 0, 3, 2, 4, 1, 0, 3, 5, 2, 1, 0, 6,
4, 3, 2, 1, 7, 0, 5, 4, 8, 3, 6, 2, 1, 0, 5, 9, 4,
7, 3, 2, 1, 6, 10, 0, 8, 5, 4, 11, 7, 3, 2, 1, 0],
dtype=int32)
m_ # m番目の零点
array([0, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 1, 1, 2, 2, 2, 1, 1, 2, 2, 3, 1,
3, 1, 2, 1, 2, 1, 3, 3, 3, 2, 1, 2, 1, 3, 3, 4, 2, 1, 4, 1, 2, 3,
1, 2, 3, 4, 4, 4], dtype=int32)
t_ # Jn(x)=0のとき０，Jn'(x)=0のとき１
array([1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 0, 0, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1,
0, 0, 1, 1, 0, 0, 1, 0, 1, 1, 1, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1,
1, 1, 0, 1, 0, 1], dtype=int32)
zo # 零点の値の配列
array([ 0. , 1.84118378, 2.40482556, 3.05423693, 3.83170597,
3.83170597, 4.20118894, 5.1356223 , 5.31755313, 5.33144277,
5.52007811, 6.3801619 , 6.41561638, 6.70613319, 7.01558667,
7.01558667, 7.50126614, 7.58834243, 8.0152366 , 8.41724414,
8.53631637, 8.57783649, 8.65372791, 8.77148382, 9.28239629,
9.64742165, 9.76102313, 9.93610952, 9.96946782, 10.17346814,
10.17346814, 10.51986087, 10.71143397, 11.06470949, 11.08637002,
11.34592431, 11.61984117, 11.7060049 , 11.73493595, 11.77087667,
11.79153444, 12.22509226, 12.3386042 , 12.68190844, 12.82649123,
12.93238624, 13.01520072, 13.17037086, 13.32369194, 13.32369194])
整数次 $n$ の第1種ベッセル関数の零点，導関数の零点，第2種ベッセル関数の零点，導関数の零点：$J_n(x)=0$，$J_n'(x)=0$，$Y_n(x)=0$，$Y_n'(x)=0$
円形導波管の遮断特性
a = 14.3
dpi = 2*np.pi*a
wlc = dpi/zo[1:]
fc = 300/wlc
print("a=",a,'[mm]')
print(f"{'No.':>7s}",f"{'mode':>7s}",f"{'xi':>10s}",f"{'wlc[mm]':>10s}",f"{'fc[GHz]':>10s}")
for i in range(1,nt):
print(f'{i:6.0f}',end=' ')
if t_[i]==1:
print(f'{"TE":2s}',end=' ')
elif t_[i]==0:
print(f'{"TM":2s}',end=' ')
print(f'{n_[i]:2.0f}',f'{m_[i]:2.0f}',f'{zo[i]:10.6f}',f'{wlc[i-1]:10.6f}',f'{fc[i-1]:10.6f}')
a= 14.3 [mm]
No. mode xi wlc[mm] fc[GHz]
1 TE 1 1 1.841184 48.799881 6.147556
2 TM 0 1 2.404826 37.362190 8.029508
3 TE 2 1 3.054237 29.418003 10.197837
4 TM 1 1 3.831706 23.448968 12.793740
5 TE 0 1 3.831706 23.448968 12.793740
6 TE 3 1 4.201189 21.386696 14.027412
7 TM 2 1 5.135622 17.495358 17.147406
8 TE 4 1 5.317553 16.896785 17.754857
9 TE 1 2 5.331443 16.852765 17.801234
10 TM 0 2 5.520078 16.276862 18.431071
11 TM 3 1 6.380162 14.082644 21.302818
12 TE 5 1 6.415616 14.004820 21.421197
13 TE 2 2 6.706133 13.398116 22.391208
14 TM 1 2 7.015587 12.807133 23.424447
15 TE 0 2 7.015587 12.807133 23.424447
16 TE 6 1 7.501266 11.977918 25.046089
17 TM 4 1 7.588342 11.840471 25.336830
18 TE 3 2 8.015237 11.209844 26.762193
19 TM 2 2 8.417244 10.674462 28.104462
20 TE 1 3 8.536316 10.525565 28.502034
21 TE 7 1 8.577836 10.474617 28.640666
22 TM 0 3 8.653728 10.382757 28.894061
23 TM 5 1 8.771484 10.243370 29.287238
24 TE 4 2 9.282396 9.679564 30.993131
25 TE 8 1 9.647422 9.313323 32.211920
26 TM 3 2 9.761023 9.204932 32.591225
27 TM 6 1 9.936110 9.042729 33.175824
28 TE 2 3 9.969468 9.012472 33.287205
29 TM 1 3 10.173468 8.831752 33.968344
30 TE 0 3 10.173468 8.831752 33.968344
31 TE 5 2 10.519861 8.540945 35.124920
32 TE 9 1 10.711434 8.388191 35.764566
33 TM 4 2 11.064709 8.120371 36.944123
34 TM 7 1 11.086370 8.104506 37.016446
35 TE 3 3 11.345924 7.919104 37.883076
36 TM 2 3 11.619841 7.732425 38.797661
37 TE 1 4 11.706005 7.675509 39.085354
38 TE 6 2 11.734936 7.656586 39.181952
39 TE 10 1 11.770877 7.633208 39.301955
40 TM 0 4 11.791534 7.619835 39.370930
41 TM 8 1 12.225092 7.349601 40.818543
42 TM 5 2 12.338604 7.281987 41.197549
43 TE 4 3 12.681908 7.084860 42.343813
44 TE 11 1 12.826491 7.004998 42.826563
45 TE 7 2 12.932386 6.947639 43.180137
46 TM 3 3 13.015201 6.903432 43.456647
47 TE 2 4 13.170371 6.822097 43.974747
48 TM 1 4 13.323692 6.743593 44.486673
49 TE 0 4 13.323692 6.743593 44.486673
円形導波管の遮断特性を求める関数は，
def circular_waveguide_xi2(imode, m, n): # imode = 1 (TE), -1 (TM)
if imode == 1:
print('TE mode:',end=' ')
xi = jnp_zeros(m,n)
elif imode == -1:
print('TM mode:',end=' ')
xi = jn_zeros(m,n)
print(xi)
return xi[n-1]
circular_waveguide_xi2(imode=1,m=1, n=3)
TE mode: [1.84118378 5.33144277 8.53631637]
8.536316366346286
第2種ベッセル関数
n = 1; nt = 5
jnyn_zeros(n ,nt)
(array([ 3.83170597, 7.01558667, 10.17346814, 13.32369194, 16.47063005]),
array([ 1.84118378, 5.33144277, 8.53631637, 11.7060049 , 14.86358863]),
array([ 2.19714133, 5.42968104, 8.59600587, 11.74915483, 14.89744213]),
array([ 3.68302286, 6.94149995, 10.12340466, 13.28575816, 16.44005801]))
前のページに戻る
「Python」の目次ページに戻る