円盤での波動方程式の差分法

中西謙太 2005年2月21日

– Typeset by FoilTEX –

２次元波動方程式

未知関数u(r, θ, t) （r, θ は空間座標,tは時刻。 0 ≤ r ≤ 1, 0 ≤ θ ≤ 2π, t ≥ 0）に関する２次元波動方程式。

\frac{∂^{2}u}{∂t^{2}} = \frac{∂^{2}u}{∂r^{2}} + \frac{1}{r}\frac{∂u}{∂r} + \frac{1}{r^{2}} \frac{∂^{2}u}{∂θ^{2}} (0 ≤ r ≤ 1, 0 ≤ θ ≤ 2π, t ≥ 0)  (1)

境界条件(Dirichlet境界条件)

u(1, θ, t) = 0 (t ≥ 0, 0 ≤ θ ≤ 2π)

初期条件 u(r, θ, 0) = f(r, θ) (0 ≤ r ≤ 1, 0 ≤ θ ≤ 2π)

\frac{∂u}{∂t}(r, θ, 0) = g(r, θ) (0 ≤ r ≤ 1, 0 ≤ θ ≤ 2π)

– Typeset by FoilTEX – 1

差分法（陽解法）

h_{r}、 h_{θ}、 τ をそれぞれr、 θ、 t 方向の格子間隔、 r_{i} = ih_{r}、 θ_{j} = jh_{θ}、 t_{n} = nτ での近似値をu_{i,j}^{n}、 (1)の微分方程式の各項を中心差分近似、２階中心差分近似すると、

\frac{u_{i,j}^{n+1} − 2u_{i,j}^{n} + u_{i,j}^{n−1}}{τ^{2}} = \frac{u_{i+1,j}^{n} − 2u_{i,j}^{n} + u_{i−1,j}^{n}}{h_{r}^{2}}

+\frac{1}{r_{i}} \frac{u_{i+1,j}^{n} − u_{i−1,j}^{n}}{2h_{r}} + \frac{1}{r_{i}^{2}} \frac{u_{i,j+1}^{n} − 2u_{i,j}^{n} + u_{i,j−1}^{n}}{h_{θ}^{2}}

(1 ≤ i ≤ N_{r} − 1, 1 ≤ j ≤ N_{θ} − 1, n ≥ 1)

– Typeset by FoilTEX – 2

\frac{τ}{h_{r}} = λ_{r}、\frac{τ}{h_{θ}} = λ_{θ} とおくと、

u_{i,j}^{n+1} = 2(1 − λ_{r}^{2} − \frac{λ_{θ}^{2}}{r_{i}^{2}})u_{i,j}^{n} + λ_{r}^{2}(u_{i+1,j}^{n} + u_{i−1,j}^{n}) + \frac{1}{r_{i}} \frac{τ^{2}}{2h_{r}}(u_{i+1,j}^{n} − u_{i−1,j}^{n})

+\frac{λ_{θ}^{2}}{r_{i}^{2}}(u_{i,j+1}^{n} + u_{i,j−1}^{n}) − u_{i,j}^{n−1}

(1 ≤ i ≤ N_{r} − 1, 1 ≤ j ≤ N_{θ} − 1, n ≥ 1)

初期条件より、

u_{i,j}^{0} = f_{i,j} (0 ≤ i ≤ N_{r}, 0 ≤ j ≤ N_{θ})

差分法（半陰スキーム）

θ 方向の微分係数を陰的に扱う。

\frac{u_{i,j}^{n+1} − 2u_{i,j}^{n} + u_{i,j}^{n−1}}{τ^{2}} = \frac{u_{i+1,j}^{n} − 2u_{i,j}^{n} + u_{i−1,j}^{n}}{h_{r}^{2}} + \frac{1}{r_{i}} \frac{u_{i+1,j}^{n} − u_{i−1,j}^{n}}{2h_{r}}

+\frac{1}{r_{i}^{2}} \frac{u_{i,j+1}^{n+1} − 2u_{i,j}^{n+1} + u_{i,j−1}^{n+1}}{h_{θ}^{2}}

(1 ≤ i ≤ N_{r} − 1, 1 ≤ j ≤ N_{θ} − 1, n ≥ 1)

ただし、 u_{i,N_{θ}}^{n} = u_{i,0}^{n}， u_{i,−1}^{n} = u_{i,N_{θ}−1}^{n}

– Typeset by FoilTEX – 4

A_{i} := (1 + \frac{2λ_{θ}^{2}}{r_{i}^{2}})I − \frac{λ_{θ}^{2}}{r_{i}^{2}}J J := ( \begin{matrix}0 & 1 & & & 1 \\ 1 & 0 & 1 \\ & ... & ... & ... \\ & & ... & ... & 1 \\ 1 & & & 1 & 0 \end{matrix} )

b^{n} = ( 2(1 − λ_{r}^{2})u_{i,0}^{n} + λ_{r}^{2}[(1 + \frac{h_{r}}{2r_{i}})u_{i+1,0}^{n} + (1 − \frac{h_{r}}{2r_{i}})u_{i−1,0}^{n}] ; ⋮ ; 2(1 − λ_{r}^{2})u_{i,j}^{n} + λ_{r}^{2}[(1 + \frac{h_{r}}{2r_{i}})u_{i+1,j}^{n} + (1 − \frac{h_{r}}{2r_{i}})u_{i−1,j}^{n}] ; ⋮ ; 2(1 − λ_{r}^{2})u_{i,N_{θ−1}}^{n} + λ_{r}^{2}[(1 + \frac{h_{r}}{2r_{i}})u_{i+1,N_{θ−1}}^{n} + (1 − \frac{h_{r}}{2r_{i}})u_{i−1,N_{θ−1}}^{n}] )-( u_{i,0}^{n−1} ; ⋮ ; u_{i,j}^{n−1} ; ⋮ ; u_{i,N_{θ−1}}^{n−1} )

U_{i}^{n+1} := (u_{i,0}^{n+1}, ⋯, u_{i,Nθ−1}^{n+1})^{T}

– Typeset by FoilTEX – 5

とすると、

A_{i}U_{i}^{n+1} = b_{i}^{n} − U_{i}^{n−1}(i = 1, 2, ⋯, N_{r} − 1)

という連立１次方程式が導かれる。また原点は、

u_{0,0}^{n+1} = 2(1 − 2λ_{r}^{2})u_{0,0}^{n} + \frac{4λ_{r}^{2}}{N_{θ}} ∑_{j=0}^{N_{θ}−1} u_{1,j}^{n} − u_{0,0}^{n−1}

実験（厳密解との比較）

厳密解uを

u(r, θ, t) = J_{0}(ν_{0,1}r) cos(ν_{0,1}t) + \frac{1}{ν_{1,1}}J_{1}(ν_{1,1}r) cos θ sin(ν_{1,1}t)

とした。初期条件ϕ(r, θ)、 ψ(r, θ)は、

ϕ(r, θ) = J_{0}(ν_{0,1}r)  (0 ≤ r ≤ 1, 0 ≤ θ ≤ 2π)
ψ(r, θ) = J_{1}(ν_{1,1}r) cos θ (0 ≤ r ≤ 1, 0 ≤ θ ≤ 2π)

(ただし、 ν_{n,m} はJ_{n}(r) = 0のm番目の根、 J_{n} はn次第一種Bessel関数) 波動方程式（半陰スキーム）グラフと厳密解との誤差を表示

– Typeset by FoilTEX – 7

安定性条件（陽解法）

安定性条件は

λ_{r}^{2} + \frac{1}{r_{i}^{2}}λ_{θ}^{2} ≤ 1

であると予想される。つまり安定条件は

λ_{r}^{2} + N_{r}^{2}λ_{θ}^{2} ≤ 1

と書き換えることができる。実際にN_{r} = 20、 N_{θ} = 100、 τ = 0.0031、 λ_{r}^{2} + N_{r}^{2}λ_{θ}^{2} = 0.977618 < 1で安定。 N_{r} = 20、 N_{θ} = 100、 τ = 0.0032、 λ_{r}^{2} + N_{r}^{2}λ_{θ}^{2} = 1.041617 > 1で不安定。

– Typeset by FoilTEX – 8

安定性条件（半陰スキーム）

このスキームが安定性の条件として、

2 − 2λ_{r}^{2} ≥ 0

であると予想される。つまり、

λ_{r}^{2} ≤ 1

実際に、 N_{r} = 20、 N_{θ} = 100、 τ = 0.045、 λ_{r}^{2} = 0.81 < 1で安定。 N_{r} = 20、 N_{θ} = 100、 τ = 0.055、 λ_{r}^{2} = 1.21 > 1で不安定。

原点（付録）

\frac{∂u}{∂t} = \frac{∂^{2}u}{∂r^{2}} + \frac{1}{r}\frac{∂u}{∂r} + \frac{1}{r^{2}} \frac{∂^{2}u}{∂θ^{2}}

を元に戻して、

\frac{∂u}{∂t} = \frac{∂^{2}u}{∂x^{2}} + \frac{∂^{2}u}{∂y^{2}}

ここで半径h_{r},中心が原点の円を描くとすると、この円がx軸y 軸と交わる点における値はu_{1,0},u_{1,\frac{N_{θ}}{4}}, u_{1,\frac{N_{θ}}{2}},u_{1,\frac{3N_{θ}}{4}}, 中心の値はu_{0,0} となる。２階中心差分近似をして、

\frac{∂^{2}u}{∂x^{2}}(0, 0, nτ) = \frac{u_{1,0} − 2u_{0,0} + u1, \frac{N_{θ}}{2}}{h_{r}^{2}} + O(h_{r})^{2}

\frac{∂^{2}u}{∂y^{2}}(0, 0, nτ) = \frac{u_{1,\frac{N_{θ}}{4}} − 2u_{0,0} + u_{1,\frac{3N_{θ}}{4}}}{h_{r}^{2}} + O(h_{r})^{2}

よってこれから、

\frac{∂^{2}u}{∂x^{2}}(0, 0, nτ) + \frac{∂^{2}u}{∂y^{2}}(0, 0, nτ) = \frac{u_{1,0} + u_{1,\frac{N_{θ}}{2}} + u_{1,\frac{N_{θ}}{2}} + u_{1,\frac{3N_{θ}}{4}} − 4u_{0,0}}{h_{r}^{2}} + O(h_{r})^{2}

と導くことができる。そこで x 軸,y 軸を少しだけ回転させる。この回転を繰り返し行うことにより、

– Typeset by FoilTEX – 11

Δu(0, 0, kτ) = \frac{u_{1,0} + u_{1,\frac{N_{θ}}{4}} + u_{1,\frac{N_{θ}}{2}}u_{1,\frac{3N_{θ}}{4}} − 4u_{0,0}}{h_{r}^{2}} + O(h_{r})^{2}

= \frac{u_{1,1} + u_{1,\frac{N_{θ}}{4} +1} + u_{1,\frac{N_{θ}}{2} +1} + u_{1,\frac{3N_{θ}}{4} +1} − 4u_{0,0}}{h_{r}^{2}} + O(h_{r})^{2}

= \frac{4(\frac{1}{N_{θ}} ∑_{j=0}^{N_{θ−1}} u_{1,j}^{n} − u_{0,0}^{n})}{h_{r}^{2}} + O(h_{r})^{2}

となる。そこで、

\frac{u_{0,0}^{n+1} − 2u_{0,0}^{n} + u_{0,0}^{n−1}}{τ^{2}} = \frac{4(\frac{1}{N_{θ}} ∑_{j=0}^{N_{θ}−1} u_{1,j}^{n} − u_{0,0}^{n})}{h_{r}^{2}}

– Typeset by FoilTEX – 12

整理すると、

u_{0,j}^{n+1} = 2(1 − 2λ_{r}^{2})u_{0,0}^{n} + \frac{4λ_{r}^{2}}{N_{θ}} ∑_{j=0}^{N_{θ}−1} u_{1,j}^{n} − u_{0,0}^{n−1}

という差分方程式が導かれる。

– Typeset by FoilTEX – 13