この記事は筆者(Python初学者)がプラズマトモグラフィーにおけるPythonを用いた級数展開法についての備忘録を示したものです。補足やミスなどがあれば遠慮なくコメントしてください。
そもそも級数展開法って？
ここでの級数展開法はプラズマのトモグラフィーなどで使用される"級数展開法"を指す。ここでの説明は簡単に済ますが、級数展開法とは像が基底関数の線形結合で表されるという仮定の元(ここでの像は再現したい2次元分布などの今わかっていない情報のこと)、得られた観測値(1次元)に基底関数(2次元)の投影図(1次元)を当てはめる(フィッティング)ことによって、展開係数を決定している。この展開係数をそれぞれ対応する基底関数にかけて基底関数の像を足し合わせると、観測値(1次元)から像(2次元)が得られるという手法である。
Fourier-Bessel関数を基底関数に用いた級数展開法を式で表そう！
文章で説明してもわかりにくいと思うので一旦式で補足を交えながら説明しよう。
g(r,\theta)=\Sigma^{M_{max}}_{m=0}\Sigma^{L_{max}}_{l=0}[a^{(c)l}_{m}cosm\theta+a^{(s)l}_{m}sinm\theta]g^{l}_{m}(r)
\\
g^{l}_{m}(r)=J_{m}(\lambda^{l+1}_{m}r)
がFourier-Bessel関数を基底関数に用いた級数展開法を示す式である。第一式からわかるように基底関数系の足し合わせであることがわかる。また第二式はBessel関数を示す。$\lambda^{l+1}_{m}$はモード数mにおける(l+1)番目のゼロ点である。
$a^{(c)l}_m,a^{(s)l}_m$は最小二乗法で決まる展開係数である。
(補足)Bessel関数とは？
円筒関数とも呼ばれる。ここでは説明を割愛するので詳しくは各自で調べていただきたい。ここでのBessel関数は第一種Bessel関数である。下にBessel関数を記述するコード。そして結果を表示する。
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import jv #第一種Bessel関数
# グラフ描画領域を設定
fig = plt.figure(figsize=(8, 5))
ax = fig.add_subplot(111)
ax.set_title("Bessel functions", size=15)
ax.grid()
ax.set_xlim(-10, 10)
ax.set_ylim(-1.5, 1.5)
ax.set_xlabel("$\lambda_{m}^{l+1}r$", size=15, labelpad=10)
ax.set_ylabel("Jm($\lambda_{m}^{l+1}r$)", size=15, labelpad=8)
x = np.linspace(-10, 10, 65)
# 第1種Bessel関数をプロット
for v in range(6):
ax.plot(x, jv(v, x), label="m={}".format(v))
# 凡例を表示
ax.legend()
実際に基底関数系をPythonで作ってみよう!
ここまで冗長に話してきたが、実際にコードを載っけて基底関数系を作ってみましょう。(多分大半の人はできればいいはずなので早よ載せろと思っているはず...)
基底関数系のような2次元分布は地図のような等高線で表して高さ(z値)ごとで色を振り分けている。そのためまず、基底関数を構成する関数を作って、等高線を描くようなシーケンスを考えよう。
まずは$J_{0}(\lambda^{0+1}_{0}r)cos0\theta$(m=0,l=0)の時のコードを示す。
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import jv
fig = plt.figure(figsize=(5,5))
x = np.linspace(0, 1, 100) #等間隔で値を０～1まで100個作成
y = np.linspace(0, 360, 400) #等間隔で値を０～360まで100個作成
x, y = np.meshgrid(x, y) #ｘ,ｙからメッシュグリッドを作成
x1 = x*np.cos(np.radians(y))
y1= x*np.sin(np.radians(y))
#第一種ベッセル関数の各モード数に対応する一番目のゼロ点のリスト
B1 = [2.40483, 3.83171, 5.13562, 6.38016, 7.58834, 8.77148, 9.93611, 11.0864, 12.2251, 13.3543, 14.4755]
m = 0 #モード数
def fz(x,y): #高さデータzを作る関数fz定義
return jv(m, x * B1[m]) * np.cos(m * np.radians(y)) #ここで基底関数決める→jv(m,r)*cos or sin(mθ)
z = fz(x,y)
ax = plt.contourf(x1,y1,z,cmap='jet') #等高線レベルに応じて色を塗る
plt.xlim(-1,1)
plt.ylim(-1,1)
plt.gca().set_aspect('equal')
plt.show()
結果を以下に示す。
あとはモード数を変えれば、周方向にFourier級数、径方向にBessel関数の特徴を持った基底関数が描ける。