２００４年度卒業研究

２次元円盤領域における波動方程式の研究
明治大学理工学部数学科
中 西 謙 太指導教官 桂田 祐史

2005 年 2 月 28 日

1

目次
1

はじめに

3

2

空間２次元波動方程式（円盤領域）の初期値境界値問題
2.1 円盤領域の初期値境界値問題 . . . . . . . . . . . . . . . . . . . . .
2.2 Bessel 関数 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3 Bessel 関数を使った厳密解 . . . . . . . . . . . . . . . . . . . . . . .
2.3.1 ２次元円盤のラプラシアンの固有値、固有関数 . . . . . . .

3
3
6
7
9

3

円盤による差分法
11
3.1 陽解法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
3.2 半陰スキーム . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144差分解の安定性について
4.1 正方形領域の場合（参考） . . . . . . . . . . . . . . . . . . . . . . .
4.2 円領域 (陽解法) の場合 . . . . . . . . . . . . . . . . . . . . . . . . .
4.2.1 厳密解との比較 . . . . . . . . . . . . . . . . . . . . . . . . .
4.3 円領域 (半陰スキーム) の場合 . . . . . . . . . . . . . . . . . . . . .
4.3.1 厳密解との比較 . . . . . . . . . . . . . . . . . . . . . . . . .

16
16
17
18
19
19

5

付録（プログラム）
5.1 陽解法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2 半陰スキーム . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.3 Bessel 関数 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

21
22
27
34

6

あとがき
35
6.1 あとがき . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
6.2 参考図書 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36

2

1

はじめに

今回私は、
「２次元円盤領域における波動方程式」をテーマとして研究した。２次元円盤領域における波動方程式の初期値境界値問題を差分法を用いて解くことを目標とした。その方法として、陽解法と半陰スキームという２つの方法を用いて解いていった。また、解の安定性については 1997 年度桂田研、養田孝さんの卒業研究レポートを参考にした。
まず、円盤領域 Ω = {(x, y); x2 + y 2 < 1} における２次元波動方程式

utt = uxx + uyy
を極座標変換し、円盤領域にして解いていくことにする。

2

空間２次元波動方程式（円盤領域）の初期値境界値問題

2.1

円盤領域の初期値境界値問題

解くべき領域の境界形状が円なので、極座標変換

x = r cos θ
y = r sin θ(0 ≤ r ≤ 1, 0 ≤ θ ≤ 2π)を行なう。ここで偏微分の関係式∂u∂r ∂u ∂θ ∂u=+= rx ur + θx uθ ,∂x∂x ∂r ∂x ∂θ∂u∂r ∂u ∂θ ∂u=+= ry ur + θy uθ∂y∂y ∂r ∂y ∂θを用いる。式 (1)、(2) の u を x、y とおくと、

rx xr + θx xθ = 1,rx yr + θx yθ = 0,ry xr + θy xθ = 0,ry yr + θy yθ = 1
3

(1)
(2)

となる。以上より rx 、θx 、ry 、θy を解くと、

rx = yθ /J,

θx = −yr /J,(3)

ry = −xθ /J,θy = xr /J(4)となる。J はヤコビアンとよばれる量で、次式で定義される。

J = xr y θ − xθ y r .
(3)、(4) を (1)、(2) に代入すると r − θ 面の微係数のみで表された関係式
ux = (yθ ur − yr uθ )/J,(5)

uy = (−xθ ur + xr uθ )/J(6)が得られる。２階微係数を求めるには、(3)、(4) を 2 度使う。すなわち、(1)、(2)の u の代わりに ux 、uy とおき、整理すると、

uxx = (yθ2 urr − 2yr yθ urθ + yr2 uθθ )/J 2 + [(yθ2 yrr − 2yr yθ yrθ + yr2 yθθ )(xθ ur − xr uθ )+(yθ2 xrr − 2yr yθ xrθ + yr2 xθθ )(yr uθ − yθ ur )]/J 3 ,uyy = (x2θ urr − 2xr xθ urθ + x2r uθθ )/J 2 + [(x2θ yrr − 2xr xθ yrθ + x2r yθθ )(xθ ur − xr uθ )+(x2θ xrr − 2xr xθ xrθ + x2r xθθ )(yr uθ − yθ ur )]/J 3が得られる。さらに uxx + uyy を求めるために 2 式を加えると、

uxx + uyy = (aurr − 2burθ + cuθθ )/J 3 + [(axrr − 2bxrθ + cxθθ )(yr uθ − yθ ur )+(ayrr − 2byrθ + cyθθ )(xθ ur − xr uθ )]/J 3 ,a = x2θ + yθ2 ,b = xr xθ + y r y θ ,c = x2θ + yθ2となる。実際に計算すると、

xr = cos θ, xθ = −r sin θ, yr = sin θ, yθ = r cos θ,
a = x2θ + yθ2 = r2 , b = xr xθ + yr yθ = 0, c = x2r + yr2 = 1,
xrr = 0, xrθ = − sin θ, xθθ = −r cos θ,
yrr = 0, yrθ = cos θ, yθθ = −r sin θとなるので、uxx + uyy は、

1r2 urr + uθθ+ 3 [−r cos θ(uθ sin θ − ur r cos θ) − r sin θ(−ur r sin θ − uθ cos θ)]
2
r
r
1
1
= urr + 2 uθθ + urrruxx + uyy =4

となり、問題は以下のようになる。
未知関数 u(r, θ, t)（r, θ は空間座標。t は時刻。0 ≤ r ≤ 1, 0 ≤ θ ≤ 2π, t ≥ 0）
に関する２次元波動方程式を考える。

∂ 2u∂ 2 u 1 ∂u
1 ∂ 2u
=
+
+
∂t2
∂r2
r ∂r r2 ∂θ2(0 < r < 1, 0 ≤ θ < 2π, t > 0)(7)

境界条件 (Dirichlet 境界条件) は、

u(1, θ, t) = 0 (t ≥ 0, 0 ≤ θ < 2π)(8)となり、初期条件は、

u(r, θ, 0) = f (r, θ) (0 ≤ r ≤ 1, 0 ≤ θ ≤ 2π)∂u(r, θ, 0) = g(r, θ) (0 ≤ r ≤ 1, 0 ≤ θ ≤ 2π)
∂t
となる。

5

(9)
(10)

2.2

Bessel 関数

この節の内容は NUMERICAL RECIPES in C を参考にした。任意の実数 ν に対して、ν 次第１種 Bessel 関数 Jν (z) は∞(− 14 z 2 )k1 ν∑Jν (z) = ( z)2k!Γ(ν + k + 1)
k=0で与えられる。ただし、これは一般には多価関数である。ν が 0 以上の整数の時には、任意の z ∈ C に対して定義される。
第２種 Bessel 関数 Yν (z) は、ν が整数でないときにはYν (z) =Jν (z) cos(νπ) − Jν (z)sin(νπ)(z ∈ R)で、ν が整数 n に等しいときはYν (z) = Yn (z) = lim Yµ (z)µ→nで定義される。
また変数 x が x < ν であるとき、特に 0 < x ≪ ν に対しては次の式で表される。
11Jν (x) ∼( x)ν (ν ≥ 0),
Γ(ν + 1) 2
2
Y0 (x) ∼ln(x),(ln(x) = loge x)πΓ(ν) 1 −ν( x)(ν > 0).
Yν (x) ∼ −π 2
次に x > ν であるとき、特に x ≫ ν に対しては次の式で表される。
√

Jν (x) ∼211
cos(x − νπ − π),
πx
2
4

√

Yν (x) ∼211
sin(x − νπ − π)πx24

また x ∼ ν での Bessel 関数の典型的な振幅はJν (ν) ∼21/3
10.4473∼,2
ν 1/3
32/3 Γ( 3 ) ν 1/3

0.774821/3
1Yν (ν) ∼ − 1/6 2 1/3 ∼ − 1/3ν
3 Γ( 3 ) νのように表される。

6

2.3

Bessel 関数を使った厳密解

まず (7),(8) の変数分離解を求める。u(r, θ, t) をu(r, θ, t) = R(r)Θ(θ)T (t)(11)

とおく。式 (7) に式 (13) を代入し、R(r)Θ(θ)T (t) で割ると、

T ′′R′′ 1 R′
1 Θ′′
=
+
+ 2
T
R
rR
r Θ

(12)

となる。式 (14) が成り立つには、左辺、右辺が、定数関数の場合だけなので、その定数値を −λ2 とおくと、T ′′R′′ 1 R′
1 Θ′′
=
+
+ 2
= −λ2 .
T
R
rR
r Θ
よって

T ′′ + λ2 T = 0(13)となる。また、Θ についても同様に

Θ′′R′′ 1 R′= −r2 (λ2 ++) = −n2ΘRrR(n : 定数)とおけるので、

Θ′′ + n2 Θ = 0.

(14)

R は、λ2 、n2 を使って整理すると、
r2 R′′ + rR′ + (λ2 r2 − n2 )R = 0(15)となる。式 (13)、(14) より、T (t) = C1 cos λt + C2 sin λt,(16)

Θ(θ) = C3 cos nθ + C4 sin nθ

(17)(C1 , C2 , C3 , C4 : 定数, n : 周期 2π の解を持つために 0 以上の整数)
式 (17) は Bessel の方程式で、その解 R(r) は、R(r) = AJn (λr) + BYn (λr)
7(A, B は定数)である。Jn は r = 0 で連続、Yn は r = 0 で不連続である。本題の条件では、r = 0
で連続なので B = 0 でR(r) = AJn (λr)
となる。一方、境界条件 u(1, θ, t) = 0 より

Jn (λ) = 0が得られる。いいかえれば、R(r) を円の境界上で０にするには、分離定数 λ をJn (λ) = 0 の根の一つであるように選ばなければならない。つまり、

λ = νn,mである。ここで、νn,m は Jn (r) = 0 の小さい方から数えて m 番目の正の根である。
以上より、R(r) = AJn (νn,m r)となる。こうして変数分離解

u(r, θ, t) = Jn (νn,m r)(C1 cos νn,m t + C2 sin νn,m t)(C3 cos nθ + C4 sin nθ)が求まった。
これらを重ね合わせすることにより波動方程式の混合問題の解が求まる。すな
わち初期値を、

f (r, θ) =

∞ ∑
∞
∑

a+
n,m Jn (νn,m r) cos nθ

n=0 m=1∞ ∑∞∑+n=1 m=1g(r, θ) =

∞ ∑
∞
∑

a−
n,m Jn (νn,m r) sin nθb+n,m Jn (νn,m r) cos nθ

n=0 m=1∞ ∑∞∑+n=1 m=1b−n,m Jn (νn,m r) sin nθとおけば、解はu(t, r, θ) =∞ ∑∞∑n=0 m=1a+n,m Jn (νn,m r) cos nθ cos(νn,m t)
8

+

∞ ∑
∞
∑
n=1 m=1∞∞ ∑∑a−n,m Jn (νn,m r) sin nθ cos(νn,m t)
1

b+
n,m

Jn (νn,m r) cos nθ sin(νn,m t)
νn,m1+
Jn (νn,m r) sin nθ sin(νn,m t)
b−
n,m
νn,m
n=1 m=1+n=0 m=1
∞
∞ ∑
∑

となる。
−
+
−
ここで f (r, θ), g(r, θ) が与えられたとき a+n,m , an,m , bn,m , bn,m を求める公式を紹介
しておく。
∫ 1

∫ π

a+
n,m =

2
π{Jn+1 (νn,m r)}2

a−
n,m =

∫ π
∫ 1
2
f (r, θ) sin nθdθ}dr,Jn (νn,m r){
π{Jn+1 (νn,m r)}2 0
−π

b+
n,m =

∫ π
∫ 1
2
g(r, θ) cos nθdθ}dr,
J
(ν
r){
n n,m
π{Jn+1 (νn,m r)}2 0
−π

b−
n,m =

∫ π
∫ 1
2
g(r, θ) sin nθdθ}dr.
J
(ν
r){
n n,m
π{Jn+1 (νn,m r)}2 0−π0Jn (νn,m r){−πf (r, θ) cos nθdθ}dr,

ちなみに、実際のプログラムでは、厳密解を+
a+0,1 = 1, b1,1 = 1,すなわち

u(r, θ, t) = J0 (ν0,1 r) cos(ν0,1 t) +1
J1 (ν1,1 r) cos θ sin(ν1,1 t)
ν1,1

とした。

2.3.1

２次元円盤のラプラシアンの固有値、固有関数

この部分は、金子「偏微分方程式入門」を参考にした。

9

円盤におけるラプラシアン △ の固有値、固有関数について説明しておく。λ2 が固有値で、U が λ2 に属する固有関数であるとは、
1
1
Urr + Ur + 2 Uθθ = −λ2 U,rr
U = 0 (円周上)2をみたすことをいう。上の議論により、固有値 λ2 = νnmに属する固有関数として、

Jn (νn,m r) cos nθ,n = 0, 1, 2, · · · ,m = 1, 2, · · ·Jn (νn,m r) sin nθ,n = 1, 2, · · · ,m = 1, 2, · · ·が得られた。
小さな固有値に属する固有関数について、節 (u = 0 をみたす点の集合) の位置を見ていく。最初の固有値は J0 (s) の第 1 零点 ν0,1 = 2.40483 · · · に対応する単純固
有値で、固有関数

J0 (ν0,1 r) cos 0θ = J0 (ν0,1 r)
は円盤の内部でいたるところ正であり、節を持たない。第 2 の固有値は J1 (s) の第1 零点 ν1,1 = 3.83171 · · · に対応するもので重根であり、固有関数

J1 (ν1,1 r) cos θ,

J1 (ν1,1 r) sin θはそれぞれ y 軸、および x 軸を節に持つ（上の式が 0 になるのは θ が 0, π/2, π, 3π/2のときだからである）。3 番目の固有値は J2 (s) の第 1 零点 ν2,1 = 5.13562 · · · に対応するもので重根であり、固有関数

J2 (ν2,1 r) cos 2θ,J2 (ν2,1 r) sin 2θ

の節はそれぞれの両対角線の 2 本となる。4 番目の固有値は J0 (s) の第 2 零点 ν0,2 =5.52008 · · · に戻り、固有関数

J0 (ν0,2 r) cos 0θ = J0 (ν0,2 r)
は J0 の第 1 零点に対応した位置 r = ν0,1 /ν0,2 = 0.43565 · · · に円周状の節を一つ持つ。
混合問題の定常波は、固有関数に長方形と同様の時間関数因子をつければ得られる。すなわち、




cos
cos




Jn (νn,m r)  または  nθ  または  (νn,m t)sinsinの形となる。ここに sin と cos を縦に書いたところはどちらかを順不同で選択する
意である。

10

3

円盤による差分法

3.1

陽解法

差分方程式の計算方法としてまず、陽解法から行う。領域の格子点を (ri , θj ) として各時間ステップごとに近似値を求める。
まず、差分近似式を導く。記号については hr 、hθ 、τ をそれぞれ r、θ、t 方向の格子間隔とする。空間変数 r の区間 [0, 1]、θ の区間 [0, 2π] をそれぞれ Nr 、Nθ 等分したものを間隔とし、時間変数 t については、τ を間隔とする。ri = ihr 、θj = jhθ 、
tn = nτ での近似値を uni,j と書く。まず (7) の微分方程式の各項を中心差分近似、２階中心差分近似し、整理すると以下のようになる。
n−1
n
un+1
uni+1,j − 2uni,j + uni−1,j1 uni+1,j − uni−1,j1 uni,j+1 − 2uni,j + uni,j−1
i,j − 2ui,j + ui,j
=
+
+
τ2
h2r
ri
2hr
ri2
h2θ
(1 ≤ i ≤ Nr − 1, 0 ≤ j ≤ Nθ − 1, n ≥ 1)
τ
= λr 、 hτθ = λθ として整理すると,
hrun+1= λ2r (uni+1,j − 2uni,j + uni−1,j ) +i,j

1 τ2 n(u− uni−1,j )ri 2hr i+1,jλ2
n−1+ 2θ (uni,j+1 − 2uni,j + uni,j−1 ) + 2uni,j − ui,jri
1 τ2 nλ2(u− uni−1,j )= 2(1 − λ2r − 2θ )uni,j + λ2r (uni+1,j + uni−1,j ) +ri
ri 2hr i+1,jλ2
n−1+ 2θ (uni,j+1 + uni,j−1 ) − ui,jri
(1 ≤ i ≤ Nr − 1, 0 ≤ j ≤ Nθ − 1, n ≥ 1)
ただし、uni,Nθ = uni,0 , uni,−1 = uni,Nθ −1 とする。
初期条件より、

u0i,j = fi,j(0 ≤ i ≤ Nr , 0 ≤ j ≤ Nθ )∂u(r, θ, 0) = g(r, θ) については Taylor の公式より、∂t∂uτ 2 ∂ 2u(r, θ, 0) + O(τ 3 )(r, θ, 0) +
∂t
2 ∂t2
τ2
1
1
= f (r, θ) + τ g(r, θ) + (frr (r, θ) + fr (r, θ) + 2 fθθ (r, θ)) + O(τ 3 )2rru(r, θ, τ ) = u(r, θ, 0) + τ

11導関数をいずれも中心差分商で近似すると、

τ 2 fi−1,j − 2fi,j + fi+1,j
1 fi+1,j − fi−1,j1 fi,j−1 − 2fi,j + fi,j+1
)
+
+ 2
(
2
2
hr
ri
2hr
ri
h2θ
1
1 τ2
= fi,j + τ gi,j + (λ2r (fi−1,j − 2fi,j + fi+1,j ) +(fi+1,j − fi−1,j )2
ri 2hrλ2+ 2θ (fi,j−1 − 2fi,j + fi,j+1 ))
ri
1
1 τ2
λ2
(fi+1,j − fi−1,j )= (1 − λ2r − 2θ )fi,j + (λ2r (fi−1,j + fi+1,j ) +ri2ri 2hrλ2+ 2θ (fi,j−1 + fi,j+1 )) + τ gi,jri

u1i,j = fi,j + τ gi,j +ただし、fi,Nθ = fi,0 , fi,−1 = fi,Nθ −1 とする。
次に原点 (r = 0) での微分方程式での近似を考える。（なぜなら 1r , r12 があるun+1i,jの近似式中では r = 0 では使えないので）∂u∂ 2 u 1 ∂u
1 ∂ 2u
= 2 +
+
∂t
∂r
r ∂r r2 ∂θ2を元に戻して、∂u∂2u ∂ 2u=+.
∂t
∂x2 ∂y 2
ここで半径 hr , 中心が原点の円を描く、この円が x 軸 y 軸と交わる点における値は
u1,0 ,u1, Nθ , u1, Nθ ,u1, 3Nθ , 中心の値は u0,0 となる。２階中心差分近似をして、
424u1,0 − 2u0,0 + u1,Nθ /2∂2u(0, 0, nτ ) =+ O(hr )2 ,2∂xh2r
u1,Nθ /4 − 2u0,0 + u1,3Nθ /4
∂ 2u
(0,
0,
nτ
)
=
+ O(hr )222∂yhr
よってこれから、

u1,0 + u1,Nθ /2 + u1,Nθ /2 + u1,3Nθ /4 − 4u0,0
∂ 2u
∂ 2u
(0,
0,
nτ
)
+
(0, 0, nτ ) =+ O(hr )2
2
2
∂x
∂y
h2r

12

以下、この回転を繰り返し行うことにより、

u1,0 + u1,Nθ /4 + u1,Nθ /2 + u1,3Nθ /4 − 4u0,0+ O(hr )22
hr
u1,1 + u1,Nθ /4+1 + u1,Nθ /2+1 + u1,3Nθ /4+1 − 4u0,1=+ O(hr )22
hr..
.
1 ∑Nθ−1 n4( Nθ j=0 u1,j − un0,0 )=+ O(hr )22
hr∆u(0, 0, nτ ) =となる。そこで、n−1n4( N1θ
un+1
0,0 − 2u0,0 + u0,0
=
τ2

∑Nθ−1
j=0

un1,j − un0,0 )h2r

整理すると、
2 n
un+10,j = 2(1 − 2λr )u0,0 +θ −1
4λ2r N∑n−1
un − u0,0.
Nθ j=0 1,jという差分方程式が導かれる。
また境界条件は、

un+1
Nr ,j = 0
である。

13

3.2

半陰スキーム

後述する陽解法の安定性条件についての予想を元に考えると、安定性の条件を厳しくしているのは、∆=∂ 2 u 1 ∂u
1 ∂2u
+
+
∂r2
r ∂r r2 ∂θ2の式の r12 の存在である。
よって θ 方向の微分係数を陰的に扱うことにする。
n−1
n+1
n+1
n
un+1
uni+1,j − 2uni,j + uni−1,j1 uni+1,j − uni−1,j1 un+1i,j − 2ui,j + ui,ji,j+1 − 2ui,j + ui,j−1
=
+
+
τ2
h2r
ri
2hr
ri2
h2θ
(1 ≤ i ≤ Nr − 1, 0 ≤ j ≤ Nθ − 1, n ≥ 1)

λ2
λ2θ n+1
hr n
hr n
n−1
2 n
2
(1 + 2 2θ )un+1(ui,j+1 + un+1
)ui+1,j + (1 −)u−
] − ui,j(18)
i,ji,j−1 ) = (2 − 2λr )ui,j + λr [(1 +
2
ri
ri
2ri
2ri i−1,j(1 ≤ i ≤ Nr − 1, 1 ≤ j ≤ Nθ − 1, n ≥ 1)
ただし、uni,Nθ = uni,0 ，uni,−1 = uni,Nθ −1 , 右辺に現れる unNr ,j は境界値から計算 (=u(1, θj , nt))、
un0,j は（j が何であっても）原点での値とする。ここで
2λ2

λ2

i

i

Ai := (1 + r2θ )I − r2θ J,


0

1



J := 









n
bi = 





,

1
0
..
.

I := Nθ 次の単位行列,
1

1
..
.
..
.

1

..
.
..
.
1







,

1


0

hr
hr
2(1 − λ2r )uni,0 + λ2r [(1 + 2r
)uni+1,0 + (1 − 2r
)uni−1,0 ]
i
i
..
.
hr
hr
2(1 − λ2r )uni,j + λ2r [(1 + 2r
)uni+1,j + (1 − 2r
)uni−1,j ]
i
i
..
.

hr
hr
)uni+1,Nθ−1 + (1 − 2r
)uni−1,Nθ−1 ]2(1 − λ2r )uni,Nθ−1 + λ2r [(1 + 2r
i
i

T
n+1
Uin+1 := (un+1i,0 , . . . , ui,Nθ −1 )

14








  n−1− u  i,j
 
..
 
.
 










n−1
ui,0

.. 

. 



n−1
ui,N
θ−1

とおくと (18) はAi Uin+1 = bni − Uin−1(i = 1, 2, · · · , Nr − 1)(19)と表される。
原点においては、陽解法のときと同様に
2 n
un+1
0,0 = 2(1 − 2λr )u0,0 +θ −1
4λ2r N∑n−1
un − u0,0Nθ j=0 1,jを採用する。
また、f (r, θ) と g(r, θ) も同様に、

f (r, θ) = J0 (ν0,1 r)(0 ≤ r ≤ 1, 0 ≤ θ ≤ 2π)g(r, θ) = J1 (ν1,1 r) cos θ(0 ≤ r ≤ 1, 0 ≤ θ ≤ 2π)とした。
(19) のタイプの係数行列を持った連立１次方程式の解き方とプログラムは桂田祐史 [熱方程式に対する差分法 II — 円盤領域、円柱領域、球における熱方程式 —]
による。

15

4

差分解の安定性について

4.1

正方形領域の場合（参考）まずは、参考として正方形領域においての安定性を挙げておく。この安定性の証明は 1997 年度桂田研, 養田孝さんの卒業論文を参考にした。
λ2x + λ2y > 1 となるように Nx 、Ny 、τ を選ぶと近似解が安定しないことがある。
実は、

λ2x + λ2y ≤ 1

τ2
τ2( 2 + 2 ≤ 1)hx hy

は解の安定性の十分条件であることが知られている。この条件を Courant-FriedrichsLewy の条件（CFL 条件）という。CFL 条件の証明は以下のとおりである。
√
uki,j = g k exp[ −1(iξhx + jηhy )] (g: 増幅率, ξ, η: 任意)
として、正方形領域の差分方程式に代入すると、√√
g k+1 exp[ −1(iξhx + jηhy )] = 2(1 − λ2x + λ2y )g k exp[ −1(iξhx + jηhy )]√+λ2x (g k exp[ −1((i − 1)ξhx + jηhy )]√+g k exp[ −1((i + 1)ξhx + jηhy )])√+λ2y (g k exp[ −1(iξhx + (j − 1)ηhy )]√+g k exp[ −1(iξhx + (j + 1)ηhy )])√−g (k−1) exp[ −1(iξhx + jηhy )]整理すると、√√
g 2 = [(exp[ −1ξhx ] − 2 + exp[− −1ξhx ])λ2x√√+(exp[ −1ηhy ] − 2 + exp[− −1ηhy ])λ2y + 2]g − 1= 2[λ2x (cos ξhx − 1) + λ2y (cos ηhy − 1) + 1]g − 1ξhxηhy= 2[λ2x (−2 sin2
) + λ2y (−2 sin2
) + 1]g − 1
2
2
よって、

g 2 − 2[1 − 2(λ2x sin2ξhx
ηhy
+ λ2y sin2)] + 1 = 0
2
2
16

(20)

いま、式 (22) をくと、(A = 1 − 2λ2x sin2 ξh2x − 2λ2y sin2 ηh2y ) とお

g 2 − 2Ag + 1 = 0
g =A±

√

A2 − 1| g |≤ 1 となるための条件を調べる。
√√
A2 > 1 のとき : A > 1 ⇒ A + A2 − 1 > 1, A < −1 ⇒ A − A2 − 1 < 1 で両方とも駄目。
A2 = 1 のとき : A = 1 ⇒ g = 1, A = 1 ⇒ g = −1
で両方とも条件を満たす。
√√√√
A2 < 1 のとき : g = A ± 1 − A2 −1 ⇒ | g |= A2 + ( 1 − A2 )2 = 1 で条件を満たす。
ゆえに、| A |≤ 1 のとき g が安定する。A の定義式を代入すると、−1 ≤ 1 − 2(λ2x sin2ξhx
ηhy
+ λ2y sin2)≤122となる。左の不等式から（右の不等式は明らか）、

λ2x sin2ξhx
ηhy+ λ2y sin2
≤1
2
2

よって、

λ2x + λ2y ≤ 1
となる。

4.2

円領域 (陽解法) の場合

正方形領域の場合と同様に円領域の CFL 条件を導くことは難しい。しかし正方形の場合と同様にして安定性を考えると、安定性条件は

λ2r +

1 2
λ ≤1
ri2 θ

(i = 1, · · · , Nr − 1)

であると予想される。いま r12 が一番大きくなるのは r が一番小さくなるとき、つi

まり、ri = N1r のときであるから、安定条件は

λ2r + Nr2 λ2θ ≤ 1
17

と書き換えることができる。
実際に Nr = 20、Nθ = 100、τ = 0.0031、λ2r + Nr2 λ2θ = 0.977618 < 1 で安定である。Nr = 20、Nθ = 100、τ = 0.0032、λ2r + Nr2 λ2θ = 1.041617 > 1 で不安定である。

4.2.1

厳密解との比較

厳密解と差分解（陽解法）との比較を行う。
厳密解を、

u(r, θ, t) = J0 (ν0,1 r) cos(ν0,1 t) +1
J1 (ν1,1 r) cos θ sin(ν1,1 t)
ν1,1とすると初期条件 f (r, θ)、g(r, θ) は、

f (r, θ) = J0 (ν0,1 r)(0 ≤ r ≤ 1, 0 ≤ θ ≤ 2π)g(r, θ) = J1 (ν1,1 r) cos θ(0 ≤ r ≤ 1, 0 ≤ θ ≤ 2π)となる。
２つの最大誤差が、Nr , Nθ , τ を変えていくごとにどのように変化するかを調べる。(x, y) = (ri cos θj , ri sin θj ) とし、(x, y) = (0.30, 0.00) を固定して誤差を計る。
また、ν0,1 = 2.404825557695773, ν1,1 = 3.831705970207512 として計算した。
実験結果は以下のようになった。

Nr
5
10
20
40
80

Nθ
τ
25
0.48
50
0.12×10−1
100
0.3×10−2
200 0.75×10−3
400 0.1875×10−4

最大誤差1.40e-007
1.64e-007
1.94e-007
2.23e-007
8.64e-008

Nr , Nθ をそれぞれ 2 倍にしていき、τ を 4 分の 1 していくとき、最大誤差のそれ
ぞれの差は、ほぼ同じ値になった。

18

4.3

円領域 (半陰スキーム) の場合

陽解法の場合と同様に、このスキームが安定であるためには、

2 − 2λ2r ≥ 0と予想できる。つまり安定条件は、

λ2r ≤ 1
実際に、Nr = 20、Nθ = 100、τ = 0.048、λ2r = 0.9216 < 1 で安定である。
Nr = 20、Nθ = 100、τ = 0.055、λ2r = 1.0404 > 1 で不安定である。
また、この結果からわかるように半陰スキームの場合、陽解法とは違い Nθ の値をどんなに小さくとろうが大きくとろうが、安定性の条件は Nθ によらないことが
わかる。

4.3.1

厳密解との比較

厳密解と差分解（半陰スキーム）との比較を行う。
陽解法と同様に、厳密解を、

u(r, θ, t) = J0 (ν0,1 r) cos(ν0,1 t) +1
J1 (ν1,1 r) cos θ sin(ν1,1 t)
ν1,1とすると初期条件 f (r, θ)、g(r, θ) も、

f (r, θ) = J0 (ν0,1 r)(0 ≤ r ≤ 1, 0 ≤ θ ≤ 2π)g(r, θ) = J1 (ν1,1 r) cos θ(0 ≤ r ≤ 1, 0 ≤ θ ≤ 2π)となる。
陽解法と同様に考え、今度は (x, y) = (0.15, 0.00) での点を固定し、ν0,1 と ν1,1 も同じ値として最大誤差を計る。
実験結果は以下のようになった。

19

Nr
5
10
20
40
80

Nθ
τ
25
0.16
50 0.8 ×10−1
100 0.4×10−1
200 0.2×10−1
400 0.1×10−1

最大誤差
1.42e-009
1.42e-009
1.39e-009
1.40e-009
1.42e-009

陽解法と同様に Nr , Nθ をそれぞれ 2 倍していき、今度は τ を 2 分の 1 倍してい
くとそれぞれの最大誤差はほぼ同じ値になった。

20

5

付録（プログラム）以下、実験に用いたプログラムを紹介する。
２次元円領域 Dirichlet 境界条件（陽解法）
２次元円領域 Dirichlet 境界条件（半陰スキーム）（なお、プログラムはすべて Java を用いて書いた。）

21

5.1

陽解法

/** ２次元波動方程式（円領域・陽解法）** <applet code = Wave2d width=500 height=500></applet>**/
import java.applet.*;import java.awt.*;
import java.awt.event.*;
import corejava.*;
import Mitsui.*;
import cern.jet.math.*;public class Wave2d extends Applet{
MitsuiWorld_2 m = new MitsuiWorld_2();
double xmin,xmax,ymin,ymax,zmin,zmax;
private TextArea ta;
public void init(){Graphics g=getGraphics();
m.setGraphics(g);m.setScreenSize(getSize().width,getSize().height);//[a,b] × [c,d] × [e,f]xmin=-2.0;
xmax=2.5;
ymin=-2.0;
ymax=2.5;
zmin=-1.0;
zmax=1.0;
m.setArea2(xmin,xmax,ymin,ymax,zmin,zmax);
m.setColor(3);add(ta = new TextArea(5,30));//ta.setBounds(0,0,10,7);}
public void paint(Graphics g){
int nr = 20;
int nt = 20;
double tau = 0.005;
double Tmax = 1.0;
double dt = 0.05;
double skip = dt/tau+0.5;
double hr = 1.0/nr;
double ht = 2.0*Math.PI/nt;
double lambdar = tau / hr;double lambdat = tau / ht;

22

double lambdar2 = lambdar * lambdar;
double lambdat2 = lambdat * lambdat;
double Kmax = (Tmax+0.1*tau)/tau;
double[][] u1 = new double[nr+1][nt+1];
double[][] u2 = new double[nr+1][nt+1];
double[][] u3 = new double[nr+1][nt+1];//初期値代入
for(int i=0;i<nr;i++){for(int j=0;j<nt;j++){u1[i][j] = phi(i*hr,j*ht);}}{
double mean = 0.0;
for(int j=0;j<nt;j++)mean += u1[1][j];mean /= nt;
for(int j=0;j<nt;j++)u1[0][j]=mean;}for(int i=0;i<=nr;i++){u1[i][nt]=u1[i][0];}for(int j=0;j<nt;j++){u1[nr][j]=0.0;}
g.clearRect(0,0,getSize().width,getSize().height);
m.hideBirdView(u1,nr,nt,1);
for(int i=1;i<nr;i++){double ri=i*hr;for(int j=0;j<nt;j++){
int jp1=j+1;
if(jp1==nt){jp1=0;}
int jm1=j-1;
if(jm1==-1){jm1=nt-1;}
u2[i][j]=(1.0-lambdar2-(lambdat2/(ri*ri)))*u1[i][j]+0.5*(lambdar2*(u1[i-1][j]+u1[i+1][j])+((tau*tau)/(2.0*ri*hr))*(u1[i+1][j]-u1[i-1][j])+(lambdat2/(ri*ri))*(u1[i][jm1]+u1[i][jp1]))+tau*psi(i*hr,j*ht);}}{
double mean = 0.0;
for(int j=0;j<nt;j++)23

mean += u2[1][j];mean /= nt;
for(int j=0;j<nt;j++)u2[0][j]=mean;}for(int i=0;i<=nr;i++){u2[i][nt]=u2[i][0];}for(int j=0;j<nt;j++){u2[nr][j]=0.0;}
g.clearRect(0,0,getSize().width,getSize().height);
m.hideBirdView(u2,nr,nt,1);
for(int k=1;k<=Kmax;k++){for(int i=1;i<nr;i++){double ri=i*hr;for(int j=0;j<nt;j++){
int jp1=j+1;
if(jp1==nt){jp1=0;}
int jm1=j-1;
if(jm1==-1){jm1=nt-1;}
u3[i][j]=2.0*(1.0-lambdar2-(lambdat2/(ri*ri)))*u2[i][j]+lambdar2*(u2[i+1][j]+u2[i-1][j])+((tau*tau)/(2.0*ri*hr))*(u2[i+1][j]-u2[i-1][j])+(lambdat2/(ri*ri))*(u2[i][jp1]+u2[i][jm1])-u1[i][j];
}
}
//r=0
{
double mean=0.0;
for(int j=0;j<nt;j++)mean += u3[1][j];mean /= nt;
for(int j=0;j<nt;j++)u3[0][j]=mean;}//Dir(r=1)
for(int j=0;j<nt;j++){u3[nr][j]=0.0;}//θ=2 πのところはθ=0 のところのコピー
for(int i=0;i<=nr;i++){u3[i][nt]=u3[i][0];}24g.clearRect(0,0,getSize().width,getSize().height);
m.hideBirdView(u3,nr,nt,1);//u1 に u2 の値を u2 に u3 の値を代入
for(int i=0;i<=nr;i++){for(int j=0;j<=nt;j++){u1[i][j] = u2[i][j];u2[i][j] = u3[i][j];}}double norm = maxnorm(u3,nr,nt);
if(norm > 10.0)
break;
double t = tau*k;}}//誤差を測る
double maxe = 0.0;
double x = 0.0;
double y = 0.0;
double e = 0.0;
for(int i=0;i<=nr;i++){double ri = i*hr;for(int j=0;j<=nt;j++){double theta_j = j * ht;
e = Math.abs(exactu(ri,theta_j,t) - u3[i][j]);if(e > maxe){
maxe = e;x = ri * Math.cos(theta_j);y = ri * Math.sin(theta_j);}}}
Format fmt = new Format("%7.2f");
Format fmt_e = new Format("%7.2e");String str =ta.getText();
if (str.length() > 1000)
ta.setText("誤差"+fmt_e.form(e)+"\n(x,y)="+fmt.form(x)+","+fmt.form(y)+")\n");else
ta.setText(str+"誤差"+fmt_e.form(e)+"\n(x,y)="+fmt.form(x)+","+fmt.form(y)+")\n");if(skip == 0){t = k*tau;}

25

//初期条件 ノルム
public double maxnorm(double[][] u,int nr,int nt){
int i0 = 0;
int j0 = 0;
double absu;
double tmpmax = Math.abs(u[0][0]);for(int i=0;i<=nr;i++){for(int j=0;j<=nt;j++){if((absu = Math.abs(u[i][j]))>tmpmax){
tmpmax = absu;
i0 = i;j0 = j;}}}return tmpmax;}
double mu01 = 2.404825557695773;
double mu11 = 3.831705970207512;
double mu02 = 5.52008;//初期条件 φpublic double phi(double r, double theta){
return Bessel.j0(mu01*r)+Bessel.j0(mu02*r);}//初期条件ψpublic double psi(double r,double theta){
return 0.0;
}

}

//厳密解
public double exactu(double r,double theta,double t){
return 1.0*Bessel.j0(mu01*r)*Math.cos(mu01*t)+2.0*Bessel.j0(mu02*r)*Math.cos(mu02*t);
}

26

5.2

半陰スキーム/** ２次元波動方程式（円領域・半陰スキーム）** <applet code = Wave2d_i2 width=500 height=500></applet>**/
import java.applet.*;import java.awt.*;
import java.awt.event.*;
import corejava.*;
import Mitsui.*;
import cern.jet.math.*;public class Wave2d_i2 extends Applet{
int nfunc=0;
MitsuiWorld_2 m = new MitsuiWorld_2();
double xmin,xmax,ymin,ymax,zmin,zmax;
private TextArea ta;
public void init(){Graphics g=getGraphics();
m.setGraphics(g);m.setScreenSize(getSize().width,getSize().height);//[a,b] × [c,d] × [e,f]xmin=-2.0;
xmax=2.5;
ymin=-2.0;
ymax=2.5;
zmin=-1.0;
zmax=1.0;
m.setArea2(xmin,xmax,ymin,ymax,zmin,zmax);
m.setColor(3);add(ta = new TextArea(5,30));//ta.setBounds(0,0,10,7);}
public void paint(Graphics g){
int nr = 20;
int nt = 100;
double tau = 0.045;
double Tmax = 3.0;
double dt = 0.05;
double hr = 1.0/nr;
double ht = 2.0*Math.PI/nt;
double lambda_r = tau / hr;double lambda_t = tau / ht;
double lambdar2 = lambda_r * lambda_r;27

double lambdat2 = lambda_t * lambda_t;
double[][] u= new double[nr+1][nt+1];
double[][] u1= new double[nr+1][nt+1];
double[][] newu = new double[nr+1][nt+1];
double[][] al= new double[nr][nt];double[][] ad= new double[nr][nt];double[][] au= new double[nr][nt];double[][] ab= new double[nr][nt];double[][] ar= new double[nr][nt];double[]b= new double[nt+1];
if(dt < tau){
dt = tau;}
double skip = Math.rint(dt/tau);g.clearRect(0,0,getSize().width,getSize().height);
m.hideBirdView(u,nr,nt,1);//LU 分解for(int i=1;i<nr;i++){double ri = i * hr;double ri2= ri * ri;
double d = 1+2*lambdat2/ri2;
double od = -lambdat2 / ri2;
for(int j=0;j<nt;j++){ad[i][j] = d;al[i][j] = au[i][j] = od;ab[i][j] = ar[i][j] = 0.0;}
ab[i][0] = ab[i][nt-2] = ar[i][0] = ar[i][nt-2] = od;
ptrilu0(nt,al[i],ad[i],au[i],ab[i],ar[i]);}for(int i=0;i<=nr;i++){double ri = i*hr;for(int j=0;j<=nt;j++){double theta_j = j*ht;
u1[i][j] = phi(ri,theta_j);}}
for(int i=0;i<=nr;i++){double ri = i*hr;for(int j=0;j<=nt;j++){double theta_j = j*ht;
u[i][j] = phi(ri,theta_j) +tau*psi(ri,theta_j);
}

28

}
//時間ループ
double nMax = Math.rint(Tmax / tau);
for(int n=1;n<=nMax;n++){for(int i=1;i<nr;i++){//右辺を作る
double alpha = lambdar2 * (1.0+0.5/i);
double beta = lambdar2 * (1.0-0.5/i);for(int j=0;j<nt;j++){b[j] = 2*(1 - lambdar2)*u[i][j] + alpha * u[i+1][j]+ beta * u[i-1][j] -u1[i][j];}//連立１次方程式を解く
ptrisol0(nt,al[i],ad[i],au[i],ab[i],ar[i],b);//コピー
for(int j=0;j<nt;j++)newu[i][j] = b[j];newu[i][nt] = b[0];}//Dirfor(int j=0;j<=nt;j++)newu[nr][j] = 0.0;//原点
double s = 0.0;
for(int j=0;j<nt;j++)s += u[1][j];newu[0][0] = 4*lambdar2 * (s/nt-u[0][0]) + 2*u[0][0] - u1[0][0];for(int j=1;j<=nt;j++)newu[0][j] = newu[0][0];// newu, u, u1 更新
for(int i=0;i<=nr;i++)for(int j=0;j<=nt;j++)u1[i][j] = u[i][j];for(int i=0;i<=nr;i++)for(int j=0;j<=nt;j++)u[i][j] = newu[i][j];g.clearRect(0,0,getSize().width,getSize().height);
m.hideBirdView(u,nr,nt,1);
double t = n * tau;//誤差を測る29

}}

double maxe,e,x=0,y=0;
maxe = 0;e = 0;
for(int i=0;i<=nr;i++){double ri = i*hr;for(int j=0;j<=nt;j++){double theta_j = j * ht;
e = Math.abs(exactu(ri,theta_j,t) - u[i][j]);if (e > maxe) {maxe = e;x = ri * Math.cos(theta_j);y = ri * Math.sin(theta_j);}}}
Format fmt = new Format("%7.2f");
Format fmt_e = new Format("%7.2e");String str =ta.getText();
if (str.length() > 1000)
ta.setText("誤差"+fmt_e.form(e)+"\n(x,y)="+fmt.form(x)+","+fmt.form(y)+")\n");else
ta.setText(str+"誤差"+fmt_e.form(e)+"\n(x,y)="+fmt.form(x)+","+fmt.form(y)+")\n");if(skip == 0){t = n*tau;}

double mu01 = 2.404825557695773;
double mu11 = 3.831705970207512;//初期値
public double phi(double r, double theta){
return Bessel.j0(mu01*r);}
public double psi(double r,double theta){
return Bessel.j1(mu11*r)*Math.cos(theta);}//厳密解public double exactu(double r,double theta,double t){
return Bessel.j0(mu01*r)*Math.cos(mu01*t)+1/mu11*Bessel.j1(mu11*r)*
Math.cos(theta)*Math.sin(mu11*t);}
public double maxnorm(int m,int n,double[][] u){30

}

int i0 = 0;
int j0 = 0;
double absu;
double tmpmax = Math.abs(u[0][0]);for(int i=0;i<m;i++){for(int j=0;j<n;j++){if((absu = Math.abs(u[i][j]))>tmpmax){
tmpmax = absu;
i0 = i;j0 = j;}}}return tmpmax;

public void ptrilu0(int n,double[] al,double[] ad,double[] au,
double[] ab,double[] ar){int nm1 = n - 1;
int nm2 = n - 2;
for(int k=0;k<nm2;k++){int kp1 = k + 1;
al[kp1] /= ad[k];ad[kp1] -= al[kp1]*au[k];if(kp1 < nm2)
ar[kp1] -= al[kp1]*ar[k];else
au[kp1] -= al[kp1]*ar[k];ab[k] = ab[k]/ad[k];if(kp1 < nm2)
ab[kp1] -= ab[k]*au[k];else
al[nm1] -= ab[k]*ar[k];ad[nm1] -= ab[k]*ar[k];}}al[nm1] /= ad[nm2];ad[nm1] -= al[nm1]*au[nm2];

public void ptrisol0(int n,double[] al,double[] ad,double[] au,
double[] ab,double[] ar,double[] b){
int nm1 = n - 1;
for(int k=0;k<nm1;k++){b[k+1] -= al[k+1]*b[k];if(k+1 < nm1 && k+1 > nm1)b[nm1] -= ab[k]*b[k];31}
b[nm1] = b[nm1]/ad[nm1];
for(int k = n-2;k>=0;k--){b[k] -= au[k]*b[k+1];
if(k+1 < nm1 && k+1 > nm1)
b[k] -= ar[k]*b[nm1];
b[k] /= ad[k];}}public void ptrilu1(int n,double[] al,double[] ad,double[] au,
double[] ab,double[] ar){int nm1 = n - 1;
for(int k=1;k<nm1;k++){int kp1 = k + 1;
al[kp1] /= ad[k];ad[kp1] -= al[kp1] * au[k];if(kp1<nm1)
ar[kp1] -= al[kp1]*ar[k];else
au[kp1] -= al[kp1]*ar[k];ab[k] /= ad[k];if(kp1<nm1)
ab[kp1] -= ab[k]*au[k];else
al[n] -= ab[k]*au[k];ad[n] -= ab[k]*ar[k];}}al[n] /= ad[nm1];
ad[n] -= al[n]*au[nm1];

public void ptrisol1(int n,double[] al,double[] ad,double[] au,
double[] ab,double[] ar,double[] b){
for(int k=1;k<n;k++){b[k+1] -= al[k+1]*b[k];if(k+1<n && k+1>n)b[n] -= ab[k]*b[k];}b[n] = b[n]/ad[n];for(int k=n-1;k>=1;k--){b[k] -= au[k]*b[k+1];
if(k+1 < n && k+1 > n)b[k] = b[k] - ar[k] * b[n];b[k] /= ad[k];
}

32

}

}

33

5.3

Bessel 関数

Bessel 関数を求めるコードついて、最初は NUMERICAL RECIPES in C を参考にしたが、後からColt Project (http://dsd.lbl.gov/~hoschek/colt/)のライブラリィを利用することにした。プログラム先頭部分で¶import cern.jet.math.*;µ³´としているのが、colt ライブラリィを使うための宣言である。

34

6

あとがき

6.1

あとがき

私は今回２次元円盤領域における波動方程式について研究をしたのだが，もう少しいろいろな差分近似の方法を試してみたかったと思う。安定性の条件についてもいろいろな方法を研究してみたかったと思う。予想ではなく式によって証明をしたかった。さらにこれから円盤領域だけでなく楕円領域についても研究してくれる桂田研の学生がいてくれることを願っている。
Java についてももう少し勉強しておくべきだったと思う。
熱心に研究を手伝って下さり、申し訳ないと思うくらいに親身に見て下さった
桂田先生に本当に深く感謝したい。
そして、４年間数学科でお世話になった先生方、そして友達に本当に「ありがとうございました。」と言いたい。この４年間で得たことを今後の自分の将来に役
立てたい。

35

6.2

参考図書

偏微分方程式：スタンリー · ファーロウ著、伊理正夫、伊理由美訳：ワイリー ·ジャパン
偏微分方程式入門：金子晃著：東京大学出版会

数値解析入門 [増訂版]：山本哲郎著：サイエンス社
応用偏微分方程式：河村哲也著：共立出版株式会社

NUMERICAL RECIPES in C 〜 C 言語による数値計算のレシピ〜：William
H.Press,Saul A.Teakolsky,William T.Vetterling,Brian P.Flannery 著、丹慶勝市、奥村晴彦、佐藤俊郎、小林誠訳 :技術評論社
養田孝、1997 年度卒業研究レポート：「２次元領域における波動方程式の研究」
桂田祐史, 熱方程式に対する差分法 II — 円盤領域、円柱領域、球における熱方
程式 —, http://www.math.meiji.ac.jp/~mk/labo/text/heat-fdm-2.pdf

36