2 - 明治大学

楕円体の酒場
明治大学理工学部数学科
嵯峨野美希
指導教員桂田祐史准教授
2014 年 3 月 3 日
39
目次
1
概要
41
2
楕円体領域での Shortley-Weller 近似
2.1 Shortley-Weller 近似とは . . . . . . . .
2.2 求める領域を差分格子で覆う . . . . .
2.3 格子点の領域内外の判定 . . . . . . . .
2.4 境界上の不等間隔格子点 . . . . . . . .
2.5 境界の座標の計算 . . . . . . . . . . . .
2.6 S-W 近似の差分方程式 . . . . . . . . .
2.7 差分スキームの安定性条件、CFL 条件
2.8 丸め誤差への対処 . . . . . . . . . . . .
2.9 第 n + 1 段の値の求め方 . . . . . . . .
.
.
.
.
.
.
.
.
.
42
42
43
44
44
45
45
46
47
48
数値実験
3.1 楕円体領域の焦点 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 初期値の設定 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3 実行結果 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
48
48
49
49
A 付録
A.1 daentai.c . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
57
57
3
40
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
1
概要
インターネット上からのコピーだが以下の様な話がある。
「シャーロックホームズはこうやって泥棒をつかまえた。」[1]
コナンドイルの名作の主人公、シャーロック・ホームズが、ある日酒場に入りました。よく
みると泥棒が何人かでテーブルを囲んでひそひそなにやら相談をしていました。ところが、声
が小さくてよく聞き取れません。よくみるとこの酒場は、天井が楕円形をしていてその一方の
焦点に泥棒たちのテーブルがあるのに気づきました。そこでホームズはこの焦点と反対の位置
（焦点）にテーブル席を取りました。すると、泥棒たちの話が手にとるように聞こえてくるで
はありませんか。こうしてなんなく泥棒の現場をおさえることが出来たといいます。このよう
に、音には、光のように反射していくというような性質も備えています。
この現象が波動方程式を解くことによって再現することに挑戦した。
2011 年度、桂田研卒業の浜勇樹先輩が 2 次元の楕円領域で検証しているが、焦点から焦点
へ波が伝わる様子を、はっきりと見ることが出来なかった。これは 2 次元の波動方程式の解で
は、ホイヘンスの原理が成り立たないためかと考えられる。そこで今回は浜先輩のプログラム
を改良し、3 次元の波動方程式を解くことによって再現できるかを試した。
x2
a2
y2
b2
< 1 を x 軸の回りに回転して出来る回転楕円体領域を Ω とする：
}
{
x2 y 2 z 2
Ω = (x, y, z); 2 + 2 + 2 < 1 .
a
b
b
√
√
この境界は F ( a2 − b2 , 0, 0), F ′ (− a2 − b2 , 0, 0) からの距離の和が 2a である点の軌跡になっ
ている。
このとき Ω において波動方程式の初期値境界値問題を考える。

1

u = uxx + uyy + uzz
((x, y, z) ∈ Ω, t > 0)

c2 tt

 u(x, y, z, 0) = φ(x, y, z)
((x, y, z) ∈ Ω)

ut (x, y, z, 0) = ψ(x, y, z)
((x, y, z) ∈ Ω)



u(x, y, z, t) = 0
((x, y, z) ∈ ∂Ω, t > 0).
楕円領域
+
c は正定数であり、初期値 φ(x, y, z), ψ(x, y, z) は既知関数である。
2001 年度桂田研卒業、関口洋正先輩の「3 次元波動方程式の差分法による数値シミュレー
ション」[2] の 3 次元波動方程式と、2011 年度桂田研卒業、浜勇樹先輩の「S-W 近似による楕
円領域での波動方程式のシミュレーション」[3] を参考にし、C 言語を用いてプログラミング
を行い、計算結果を GLSC を用いて可視化した。
41
楕円体領域での Shortley-Weller 近似
2
Shortley-Weller 近似とは
2.1
通常、差分法は等間隔格子点を用いて行われる。これは領域が直方体の場合には問題ないが、
円盤領域などには合わない。これらの領域で問題を解くには、変数変換によって領域を長方形
に変換することで差分法を適用するのが普通である。しかし、以下説明する Shortley-Weller
近似ではそのままの領域に対して差分法を行うことが出来る。簡単のため 1 変数関数 f の場合
で説明する。
1
1
1
f (x + h) = f (x) + f ′ (x)h + f ′′ (x)h2 + f (3) (x)h3 + f (4) (x + θh)h4 ,
2
3!
4!
1
1
1
f (x − h′ ) = f (x) − f ′ (x)h′ + f ′′ (x)h′2 − f (3) (x)h′3 + f (4) (x − θ′ h′ )h′4 .
2
3!
4!
(1)
(2)
(1)/h + (2)/h′ を整理すると
f (x + h) − f (x) f (x) − f (x − h′ ) 1 ′′
−
= f (x)(h + h′ )
h
h′
2
1
1
+ f (3) (x)(h2 − h′2 ) + (f (4) (x + θh)h3 + f (4) (x − θ′ h′ )h′3 ).
3!
4!
ゆえに
2
f (x) =
h + h′
′′
(
f (x + h) − f (x) f (x) − f (x − h′ )
−
h
h′
)
(
′
− O(h − h) + O
h3 + h′3
h + h′
)
.
となり、この近似を f ′′ (x) の Shortley-Weller 近似と呼ぶ。以下では、Shortley-Weller 近似の
ことを S-W 近似と略記する。
42
2.2
求める領域を差分格子で覆う
楕円体領域 Ω を覆う直方体領域の左辺の x 座標を xmin , 右辺の x 座標を xmax , 下辺の y 座標
を ymin , 上辺の y 座標を ymax 、z 方向も同様に zmin , zmax とおく。x 方向の分割数を Nx , y 方
向の分割数を Ny , z 方向の分割数を Nz , x 方向の格子点の間隔を hx , y 方向の格子点の間隔を
hy ,z 方向の分割数を hz , 時間の刻み幅 τ > 0、とおく。
このとき
xmax − xmin
hx =
Nx
ymax − ymin
hy =
Ny
hz =
zmax − zmin
.
Nz
が成り立つ。さらに、xi , yj , zk を
xi = xmin + ihx
(0 ≤ i ≤ Nx )
yj = ymax + jhy
(0 ≤ j ≤ Ny )
zk = zmax + khz .
(0 ≤ k ≤ Nk )
とする。
43
2.3
格子点の領域内外の判定
楕円体領域 Ω に対して、関数 F を F (x, y, z) = 1 −
P = (x, y, z) に対して
(
x2
a2
+
y2
b2
+
z2
b2
)
で定めると、任意の



F (x, y, z) > 0 のとき領域内部
F (x, y, z) = 0 のとき境界上


F (x, y, z) < 0 のとき境界外部
が成り立つ。この F により任意の点が領域の内部、外部、境界のどこにあるか判断できる。
2.4
境界上の不等間隔格子点
等間隔の格子点によって分けられた領域内のある格子点 P = (x, y, z) に対して、その東西
北南上下の等間隔格子点が領域外に出てしまう場合、境界上に新しく不等間隔格子点を取る。
東側、西側、南側、北側、上側、下側の点をそれぞれ
Pe = (xe , y, z), Pw = (xw , y, z), Ps = (x, ys , z), Pn = (x, yn , z), Pt = (x, y, zt ), Pu = (x, y, zu )
と置く。新たに出来た東西南北上下の格子点の間隔をそれぞれ
he = xe − x = εw hx ,
hw = x − xw = εw hx
hs = y − ys = εs hy ,
hn = yn − y = εn hy
ht = zt − z = εt hz ,
hu = z − zu = εu hz .
(0 < εe , εw , εn , εs , εt , εu ≤ 1)
と置く。
44
2.5
境界の座標の計算
格子点の間に存在する境界の座標を計算する。楕円体領域 Ω の境界上では
x2 y 2 z 2
+ 2 + 2 = 1.
a2
b
b
が満たされ、x, y, z, それぞれについて解くと
√
√
y2 z2
x2 z 2
x = ±a 1 − 2 − 2 ,
y = ±b 1 − 2 − 2 ,
b
b
a
b
√
z = ±b
1−
x2 y 2
− 2.
a2
b
となる。これによって境界上の「pe , pw の x 座標 xe , xw 」と「pn , ps の y 座標 yn , ys 」と「pt , pu
の z 座標 zt , zu 」が分かるので、格子点の東西南北上下のどの等間隔格子点が境界外にでるか
に合わせて、Pe , Pw , Ps , pn , Pt , Pu を決める。
2.6
S-W 近似の差分方程式
通常の差分法では、ラプラシアンを
△u ∼
=
ue − 2u + uw un − 2u + us ut − 2u + uu
+
+
h2x
h2y
h2z
と近似する。S-W 近似では、
(
)
(
)
(
)
2
ue − u u − uw
2
un − u u − us
2
ut − u u − uu
∼
△u =
−
−
−
+
+
he + hw
he
hw
hn + hs
hn
hs
ht + hu
ht
hu
と近似する。波動方程式
1
utt = △u
c2
の左辺を２階の中心差分で近似すると、差分方程式は
(
)
n+1
n−1
ue − uni,j,k uni,j,k − uw
1 ui,j,k − 2uni,j,k + ui,j,k
2
−
=
c2
τ2
he + hw
h
hw
( e n
)
(
)
n
u
u
−
u
ut − uni,j,k uni,j,k − ut
2
2
n
i,j,k − us
i,j,k
+
−
−
+
hn + hs
hn
hs
ht + hu
ht
hu
という形になる。un+1 について整理すると
(
)
ue − uni,j,k uni,j,k − uw
2τ 2 c2
n−1
n+1
n
ui,j,k = 2ui,j,k − ui,j,k +
−
he + hw
he
hw
(
)
(
)
n
n
2 2
un − ui,j,k ui,j,k − us
ut − uni,j,k uni,j,k − uu
2τ c
2τ 2 c2
−
−
+
+
hn + hs
hn
hs
ht + hu
ht
hu
となる。これは he = hw = hx hn = hs = hy , ht = hu = hz のとき、通常の差分方程式に帰着
する。
45
2.7
差分スキームの安定性条件、CFL 条件
Wikipedia[4] には
In mathematics, the Courant-Friedrichs-Lewy condition (CFL condition) is a necessary condition for stability while solving certain partial diﬀerential equations
(usually hyperbolic PDEs) numerically by the method of finite diﬀerences. It arises
in the numerical analysis of explicit time integration schemes, when these are used
for the numerical solution. As a consequence, the time step must be less than a
certain time in many explicit time-marching computer simulations, otherwise the
simulation will produce incorrect results.
とある。
Courant-Friedrichs-Lewy 条件 (CFL 条件) とは、ある種の偏微分方程式 (通常双曲線型偏微
分方程式) を差分法により数値的に解く際に現れる。差分解で「情報が伝播する速さ」を「実
際の現象で波や物理量が伝播する速さ」よりも速くしなければならない、ということである。
今回の波動方程式の場合に当てはめて考えると、等間隔格子を使った差分法の場合
( )2 ( )2 ( )2
cτ
cτ
cτ
+
+
≤1
hx
hy
hz
となる。これを τ について解くと
√
τ≤
h2x h2y h2z
c2 (h2y h2z + h2x h2z + h2x h2y )
という τ の条件が得られる。
次に S-W 近似における安定性条件について考える。まず
(
)
(
)
(
)
2
ue − u u − uw
2
un − u u − us
2
ut − u u − un
∼
△u =
−
−
−
+
+
he + hw
he
hw
hn + hs
hn
hs
ht + hu
ht
hn
に he = εe hx , hw = εw hx , hn = εn hy , hs = εs hy , ht = εt hz , hu = εu hz を代入すると
△u ∼
=
2
εe hx + εw hx
(
ue − u u − uw
−
εe hx
εw hx
)
(
)
2
un − u u − us
−
+
εn hy + εs hy
εn hy
εs hy
(
)
2
ut − u u − uu
−
+
εt hz + εt hz
εu hz
εu hz
これを変形すると、
△u ∼
=
2(εe uw + εw ue )
2(εn us + εs un )
2(εt uu + εu ut )
+
+
2
2
εe εw hx (εe + εw ) εn εs hy (εn + εs ) εt εu h2z (εt + εu )
2(εs εn h2y εn εt h2z + εe εw h2x εt εu h2z + εe εw h2x εs εn h2y )
u
−
εe εw εs εn εt εu h2x h2y h2z
46
となる。ここで
utt ∼
=
n−1
n
un+1
i,j,k − 2ui,j,k + ui,j,l
τ2
であるから、差分方程式は、
n−1
n+1
1 ui,j,k − 2uni,j,k + ui,j,k
2(εn us + εs un )
2(εt uu + εu ut )
2(εw ue + εe uw )
+ 2
+ 2
= 2
2
2
c
τ
hx εe εw (εe + εw ) hy εs εn (εs + εn ) hz εu εt (εu + εt )
2(εs εn hy 2 εu εt hz 2 + εe εw hx 2 εu εt hz 2 + εe εw hx 2 εn εs hy 2 ) n
−
ui,j,k
εe εw εs εn εu εt hx 2 hy 2 hz 2
となる。un+1 について解くと、
{
}
2(εw ue + εe uw )
2(εn us + εs un )
2(εt uu + εu ut )
n+1
n−1
ui,j,k = −ui,j,k +
+ 2
+ 2
c2 τ 2
2
hx εe εw (εe + εw ) hy εs εn (εs + εn ) hz εu εt (εu + εt )
{
}
2c2 τ 2 (εs εn hy 2 εu εt hz 2 + εe εw hx 2 εu εt hz 2 + εe εw hx 2 εn εs hy 2 )
+ 2−
uni,j,k
εe εw εs εn εu εt hx 2 hy 2 hz 2
となる。第三項の係数が正になるようにすると、
2≥
2c2 τ 2 (εs εn hy 2 εu εt hz 2 + εe εw hx 2 εu εt hz 2 + εe εw hx 2 εn εs hy 2 )
εe εw εs εn εu εt hx 2 hy 2 hz 2
となる。これを τ について整理すると
√
τ≤
εe εw εs εn εu εt hx 2 hy 2 hz 2
c2 (εs εn hy 2 εu εt hz 2 + εe εw hx 2 εu εt hz 2 + εe εw hx 2 εn εs hy 2 )
となる。εe = εw = εn = εs = εt = εu = 1 のとき、この式は上に述べた等間隔格子の CFL 条
件に一致している。２次元波動方程式を S-W 近似によって解く場合の CFL 条件であると考え
られる。
また、
ε = min{εe , εw , εn , εs , εt , εu }
とおく。
2.8
丸め誤差への対処
格子点が境界上にあるか、または外部であるが境界に非常に近い場合に、プログラムが丸め
誤差のせいで領域の内部だと判断することがある。よってその対策をプログラムに組み込む必
要がある。何の対策もしなければ、その等間隔格子点上を領域の内部だと勘違いをして ε を極
端に小さい値として返す。double 型では、だいたい 10−14 ∼10−15 ぐらいである。この対策と
して εr = 10−14 ∼10−15 というものを置き、これを用いて点 (x, y, z) が領域の内部、外部、境
界のどこに位置するのか判断する。
47



F (x, y, z) > εr
のとき領域内部
|F (x, y, z)| ≤ εr


F (x, y, z) < −ε
r
のとき境界上
のとき領域外部
と判断することにする。また、この誤差は、領域の形や Nx , Ny , Nz の大小によって多少変化す
る。なので、もし ε が εr を下回るとき、そのことを出力するようにプログラムを書いておく。
2.9
第 n + 1 段の値の求め方
差分方程式を見ると第 n + 1 段での値を求めるために、一段前の n での値のみならず、もう
一段前の n − 1 での値が必要になる。これは、もとの方程式が時刻 t に関して２階であること
に対応している。そのため計算結果を出発させるためには、n = 0 での値だけでなく、n = 1
での値も必要になるため、計算を行う。x, y, z を固定した１変数関数 t 7→ u(x, y, z, t) の t = 0
におけるテーラー展開を行う。
∂u
τ 2 ∂ 2u
u(x, y, z, t) ∼
(x, y, z, 0).
= u(x, y, z, 0) + τ (x, y, z, 0) +
∂t
2 ∂t2
の右辺第 3 項に
1 ∂ 2u
(x, y, z, 0) = △u(x, y, z, 0) = △ϕ.
c2 ∂t2
を代入する。ゆえに
τ 2 c2
∂u
u(x, y, z, τ ) ∼
(x, y, z, 0) +
△u(x, y, z, 0)
= u(x, y, z, 0) + τ
∂t
2
∼
= ϕ(x, y, z) + τ ψ(x, y, z)+
(
)
(
)
(
))
(
2
2
ue − u
u − uw
un − u
u − us
ut − u
u − uu
τ 2 c2
2
+
+
.
−
−
−
2
hw + he
he
hw
hn + hs
hn
hs
ht + hu
ht
hu
これを整理すると
u1i,j,k =ϕ(xi , yj , zk ) + τ ψ(xi , yj , zk )
(
)
(
)
(
)}
{
ue − u0i,j,k
u0i,j,k − uw
un − u0i,j,k
u0i,j,k − us
ut − u0i,j,k
u0i,j,k − uu
τ 2 c2
τ 2 c2
τ 2 c2
+
+
.
+
−
−
−
hw + he
he
hw
hn + hs
hn
hs
hu + ht
ht
hu
という式が得られる。これで数値計算の準備ができた。
3
3.1
数値実験
楕円体領域の焦点
領域 Ω で波動方程式を解くにあたって、(x0 , 0, 0) を a > b の場合の x 軸上にある正の値の
焦点の座標とする。以下の実験では a = 5, b = 3 とする。すなわち
{
}
x2 y 2 z 2
Ω = (x, y, z);
+
+
<1 .
25
9
9
√
√
√
x0 = a2 − b2 = 7 であり、( 7, 0, 0) は楕円の焦点となる。
48
3.2
初期値の設定
検証したいのは、シャーロックホームズの話が再現できるか、つまり、片方の焦点から発せ
られた波がもう片方の焦点に向かって反射し収束するのかということである。よって片方の焦
点上にピークをもち領域内のいたるところで０になるような初期値をあたえる。この実験では
このような初期値を２種類用意した。
(nfunc=1)
a < b の場合、H(x, y, z) = K exp(−M {(x − x0 )2 + y 2 + z 2 }) を用意し、
φ(x, y, z) =
{
H(x, y, z) − 0.1 (H(x, y, z) ≥ 0.1)
0
(H(x, y, z) < 0.1).
とし、ψ については
ψ(x, y, z) = 0.
とした。
(nfunc=2)
f (r) =
{
(r2 − 1)4
(r < 1)
(r ≥ 1),
)
(
∥x∥
.
g(x, y, z) = Kf
0.5
0
を用意する。g は x = 0 で最大値 K 、∥x∥ > 0.5 で f (x) = 0 となる。
φ(x, y, z) = g(x − 4, y, z).
とし、
ψ(x, y, z) = 0.
とした。K, M で初期値の波のサイズ（幅、高さ）を変更できる。今回は K = 3, M = 10 と
した。
3.3
実行結果
楕円は平面上の二つの焦点からの距離の和が一定であるという性質を持っている。長軸 a = 5、
短軸 b = 3、波の速さ c = 1 という条件から片方の焦点から出た波が境界で反射し、もう片方
の焦点へ行くのにかかる時間は 10 と予想されるため、t = 10 と t = 20 の場合とその前後の結
果を選び載せた。
また、格子の分割数と時間の刻み幅を
49
(1)Nx = 150, Ny = Nz = 90, τ = 0.002 のとき
(2)Nx = Ny = Nz = 162, τ = 0.002 のとき
(3)Nx = 150, Ny = Nz = 90, τ = 0.004 のとき
と変えて実行した。
(1),(2) の結果は、焦点から出た波が t = 10 あたりでもう片方の焦点に行く様子が何となく
見られた。しかし初期値の波とは形がちがい、波の大きさも少し小さかった。そして t = 20
で元の焦点に戻ってくると予想したが、実際に戻ってきたのは、t = 22.5 あたりである。(3)
は CFL 条件を守らなかったため、境界に近いところで崩壊した。
以下に (2) と (3) の様子を載せた。
50
The numerical simulation of a 3-dimensional wave wquation
Nx=162, Ny=162, Nz=162, tau=0.0004, Tmax=25, dt=0.01
t = 0.500000
t=0
t = 1.000000
t = 1.500000
t = 2.000000
t = 2.500000
t = 3.000000
t = 3.500000
t = 4.000000
t = 4.500000
51
t = 5.000000
t = 5.500000
t = 6.000000
t = 6.500000
t = 7.000000
t = 7.500000
t = 8.000000
t = 8.500000
t = 9.000000
t = 9.500000
52
t = 10.000000
t = 10.500000
t = 11.000000
t = 11.500000
t = 12.000000
t = 12.500000
t = 13.000000
t = 13.500000
t = 14.000000
t = 14.500000
53
t = 15.000000
t = 15.500000
t = 16.000000
t = 16.500000
t = 17.000000
t = 17.500000
t = 18.000000
t = 18.500000
t = 19.000000
t = 19.500000
54
t = 20.000000
t = 20.500000
t = 21.000000
t = 21.500000
t = 22.000000
t = 22.500000
t = 23.000000
t = 23.500000
55
The numerical simulation of a 3-dimensional wave wquation
Nx=150, Ny=90, Nz=90, tau=0.004, Tmax=1.5m, dt=0.01
t = 0.400000
t=0
t = 1.200000
t = 0.800000
56
A
A.1
付録
daentai.c
/*
* daentai.c --- 楕円体での波動方程式を SW 近似で解く
*/
#include <stdio.h>
#include <math.h>
#define M (10)
#define K (3)
double a, b, c, xmin, xmax, ymin, ymax, zmin, zmax, distance, x0, y00;
double pi;
/* to use matrix, new_matrix() */
#include <matrix.h>
#include <unistd.h>
/* to use GLSC */
#define G_DOUBLE
#include <glsc.h>
/* #include <glsc2.h> */
/* 2 乗の関数 */
double d(double x)
{
return x * x;
}
/* 領域の内部・境界・外部の判定用 */
double F(double x, double y, double z)
{
return 1 - (d(x / a) + d(y / b) + d(z / c));
}
/* 東側の境界 */
double E(double x, double y, double z)
{
return a * sqrt(1 - d(y / b) - d(z / c));
57
}
/* 西側の境界 */
double W(double x, double y, double z)
{
return - a * sqrt(1 - d(y / b) - d(z / c));
}
/* 北側の境界 */
double N(double x, double y, double z)
{
return b * sqrt(1 - d(x / a) - d(z / c));
}
/* 南側の境界 */
double S(double x, double y, double z)
{
return - b * sqrt(1 - d(x / a) - d(z / c));
}
/* 上側の境界 */
double T(double x, double y, double z)
{
return c * sqrt(1 - d(x / a) - d(y / b));
}
/* 下側の境界 */
double U(double x, double y, double z)
{
return - c * sqrt(1 - d(x / a) - d(y / b));
}
/* 大小比較の関数 */
double min(double x, double y)
{
if (x < y)
return x;
else
return y;
}
58
/* 初期条件： φ、ψ */
double H(double x, double y, double z)
{
return K*exp(-M*(d(x-x0) + d(y-y00) + d(z)));
}
double f(double r)
{
if (r < 1)
return pow(r*r-1, 4.0);
else
return 0;
}
double g(double x, double y, double z)
{
return f(sqrt(x*x+y*y+z*z)/0.1);
}
double phi(double x, double y, double z, int nfunc)
{
double HH=H(x,y,z);
if (nfunc == 1) {
if( HH >= 0.1 )
return HH - 0.1;
else
return 0.0;
}
else if (nfunc == 2) {
return K*g(x-x0,y-y00,z);
}
else
return 0;
}
double psi(double x, double y, double z, int nfunc)
{
return 0.0;
}
59
#define MAXN (400)
#define MAXNX MAXN
#define MAXNY MAXN
#define MAXNZ MAXN
double u1[MAXNX + 1][MAXNY + 1][MAXNZ + 1];
double u2[MAXNX + 1][MAXNY + 1][MAXNZ + 1];
double u3[MAXNX + 1][MAXNY + 1][MAXNZ + 1];
int main()
{
int Nx, Ny, Nz, i, j, k, n, skip, LMax;
int nfunc;
double hx, hy, hz, lambdax, lambday, lambdaz, lambdax2, lambday2,
lambdaz2, p;
double tau, Tmax, t, C, dt, ew, es, ee, en, et, eu, er, mine, maxtau;
double x, y, z, hw, he, hn, hs, ht, hu;
double xi, xiw, xie, yj, yjs, yjn, zk, zkt, zku;
double uw, ue, un, us, ut, uu;
char str[100];
double sheta;
matrix u;
char message[100];
pi = 4 * atan(1.0);
/* 波の速さ */
C = 1.0;
/* a,b,c の値 */
a = 5.0;
b = 3.0;
c = b;
/* 直方体領域の端の座標 */
xmin = - a * 1.1;
xmax = a * 1.1;
ymin = - b * 1.1;
ymax = b * 1.1;
zmin = - c * 1.1;
zmax = c * 1.1;
60
/* 分割数 */
printf("Nx= ");
scanf("%d", &Nx);
printf("Ny= ");
scanf("%d", &Ny);
printf("Nz= ");
scanf("%d", &Nz);
u = new_matrix(Nx + 1, Ny + 1);
/* 初期条件のピークの位置 */
x0 = sqrt(d(a) - d(b));
y00 = 0.0;
/* 距離 */
distance = 60;
/* 格子間の長さ */
hx = (xmax - xmin) / Nx;
hy = (ymax - ymin) / Ny;
hz = (zmax - zmin) / Nz;
printf("hx=%g, hy=%g, hz=%g\n", hx, hy, hz);
/* 誤差 */
er = 1.0e-13;
/****************最小εと最大τを求める *****************/
mine = 1.0;
maxtau = 1.0;
for (i = 0; i <= Nx; i++) {
xi = xmin + i * hx;
xiw = xmin + (i - 1) * hx;
xie = xmin + (i + 1) * hx;
for (j = 0; j <= Ny; j++) {
yj = ymin + j * hy;
yjs = ymin + (j - 1) * hy;
yjn = ymin + (j + 1) * hy;
for (k = 0; k <= Nz; k++) {
61
zk = zmin + k * hz;
zku = zmin + (k - 1) * hz;
zkt = zmin + (k + 1) * hz;
if (F(xi, yj, zk) >= er) {
/* WEST */
if (F(xiw, yj, zk) >= er) {
hw = hx;
ew = 1.0;
} else if (fabs(F(xiw, yj, zk)) < er) {
ew = 1.0;
hw = hx;
} else {
hw = xi - W(xi, yj, zk);
ew = hw / hx;
}
/* SOUTH */
if (F(xi, yjs, zk) >= er) {
hs = hy;
es = 1.0;
} else if (fabs(F(xi, yjs, zk)) < er) {
hs = hy;
es = 1.0;
} else {
hs = yj - S(xi, yj, zk);
es = hs / hy;
}
/* EAST */
if (F(xie, yj, zk) >= er) {
he = hx;
ee = 1.0;
} else if (fabs(F(xie, yj, zk)) < er) {
ee = 1.0;
he = hx;
} else {
he = E(xi, yj, zk) - xi;
ee = he / hx;
}
/* NORTH */
if (F(xi, yjn, zk) >= er) {
62
en = 1.0;
hn = hy;
} else if (fabs(F(xi, yjn, zk)) < er) {
en = 1.0;
hn = hy;
} else {
hn = N(xi, yj, zk) - yj;
en = hn / hy;
}
/* TOP */
if (F(xi, yj, zkt) >= er) {
et = 1.0;
ht = hz;
} else if (fabs(F(xi, yj, zkt)) < er) {
et = 1.0;
ht = hz;
} else {
ht = T(xi, yj, zk) - zk;
et = ht / hz;
}
/* UNDER */
if (F(xi, yj, zku) >= er) {
eu = 1.0;
hu = hz;
} else if (fabs(F(xi, yj, zku)) < er) {
eu = 1.0;
hu = hz;
}else {
hu = zk - U(xi, yj, zk);
eu = hu / hz;
}
} else {
ew = ee = en = es = et = eu = 1.0;
hw = he = hx;
hn = hs = hy;
ht = hu = hz;
}
mine =
min(mine,
min((min(ew, ee), min(en, es)), min(et, eu)));
maxtau = min(maxtau,
63
sqrt((ew * es * ee * en * eu
d(hy) * d(hz)) / (es *
et *
ee *
et *
ee *
es *
* et * d(hx)
en * d(hy) *
d(hz) +
ew * d(hx) *
d(hz) +
ew * d(hx) *
d(hy))));
*
eu *
eu *
en *
}
}
}
printf("εの最小値= %g\n", mine);
printf("τの最大値=%g\n", maxtau);
if (mine < er) {
printf("mine=%g: *******εが小さすぎます********\n", mine);
exit(-1);
}
printf("τ (限界は%g) : ", maxtau);
scanf("%lf", &tau);
lambdax = tau / hx;
lambday = tau / hy;
lambdaz = tau / hz;
printf("λ x = %g になりました。\n", lambdax);
printf("λ y = %g になりました。\n", lambday);
printf("λ z = %g になりました。\n", lambdaz);
lambdax2 = lambdax * lambdax;
lambday2 = lambday * lambday;
lambdaz2 = lambdaz * lambdaz;
p = lambdax2 + lambday2 + lambdaz2;
printf("λ x * λ x + λ y * λ y + λ z * λ z = %g になりました。\n", p);
printf("Tmax = : ");
scanf("%lf", &Tmax);
printf("Δ t= : ");
scanf("%lf", &dt);
printf("nfunc: ");
scanf("%d", &nfunc);
skip = rint(dt / tau);
/* GLSC の開始 */
g_init("1", 250.0, 180.0);
64
g_device(G_BOTH);
g_def_scale(0,
- a, a, - a * (140.0 / 80.0), a * (140.0 / 80.0),
150.0, 35.0, 80.0, 140.0);
g_def_line(0, G_BLACK, 0, G_LINE_SOLID);
g_def_line(1, G_BLACK, 0, G_LINE_DOTS);
g_def_text(0, G_BLACK, 2);
g_sel_scale(0);
skip = rint(dt / tau);
if (skip == 0) {
printf("Δ t が小さすぎるので、Δ t=
skip = 1;
dt = skip * tau;
}
τ
とします。\n");
/* 初期値の入力 */
t = 0.0;
for (i = 0; i <= Nx; i++) {
xi = xmin + i * hx;
for (j = 0; j <= Ny; j++) {
yj = ymin + j * hy;
for (k = 0; k <= Nz; k++) {
zk = zmin + k * hz;
u1[i][j][k] = phi(xi, yj, zk, nfunc);
}
}
}
/* t=0 のグラフ */
g_cls();
g_text(10.0, 10.0,
"The numerical simulation of a 3-dimensional wave wquation");
sprintf(message, "Nx=%d, Ny=%d, Nz=%d, tau=%g, Tmax=%gm, dt=%g",
Nx, Ny, Nz, tau, Tmax, dt);
g_sel_text(0);
g_text(10.0, 20.0, message);
sprintf(message, "t=%g", 0 * tau);
g_sel_text(0);
g_text(10.0, 30.0, message);
for (i = 0; i <= Nx; i++)
65
for (j = 0; j <= Ny; j++)
u[i][j] = u1[i][j][Nz / 2];
/* 座標軸を描く */
g_sel_line(1);
g_move(-0.1, 0.0);
g_plot(1.1, 0.0);
g_move(0.0, -0.1);
g_plot(0.0, 2.0);
/* 断面図を描く */
g_sel_line(0);
g_move(xmin, u[0][0]);
for (i = 1; i <= Nx; i++)
g_plot(xmin + i * hx, u[i][Ny / 2]);
/* 鳥瞰図を描く */
g_hidden2(5.0, 5.0, 1.0, -1.0, 1.0, 15.0, 25.0, 30.0, 20.0, 60.0,
100.0, 80.0, u, Nx + 1, Ny + 1, 1, G_SIDE_NONE, 4, 4);
printf("u1 まで描けた \n");
/* u2 の計算 */
for (i = 0; i <= Nx; i++) {
xi = xmin + i * hx;
xiw = xmin + (i - 1) * hx;
xie = xmin + (i + 1) * hx;
for (j = 0; j <= Ny; j++) {
yj = ymin + j * hy;
yjn = ymin + (j + 1) * hy;
yjs = ymin + (j - 1) * hy;
for (k = 0; k <= Nz; k++) {
zk = zmin + k * hz;
zkt = zmin + (k + 1) * hz;
zku = zmin + (k - 1) * hz;
if (F(xi, yj, zk) >= er) {
/* WEST */
if (F(xiw, yj, zk) >= er) {
hw = hx;
uw = u1[i - 1][j][k];
} else if (fabs(F(xiw, yj, zk)) < er) {
hw = hx;
uw = u1[i - 1][j][k];
} else {
66
hw = xi - W(xi, yj, zk);
uw = 0;
}
/* EAST */
if (F(xie, yj, zk) >= er) {
he = hx;
ue = u1[i + 1][j][k];
} else if (fabs(F(xie, yj, zk))
he = hx;
ue = u1[i + 1][j][k];
} else {
he = E(xi, yj, zk) - xi;
ue = 0;
}
/* NORTH */
if (F(xi, yjn, zk) >= er) {
hn = hy;
un = u1[i][j + 1][k];
} else if (fabs(F(xi, yjn, zk))
hn = hy;
un = u1[i][j + 1][k];
} else {
hn = N(xi, yj, zk) - yj;
un = 0;
}
/* SOUTH */
if (F(xi, yjs, zk) >= er) {
hs = hy;
us = u1[i][j - 1][k];
} else if (fabs(F(xi, yjs, zk))
hs = hy;
us = u1[i][j - 1][k];
} else {
hs = yj - S(xi, yj, zk);
us = 0;
}
/* TOP */
if (F(xi, yj, zkt) >= er) {
ht = hz;
ut = u1[i][j][k + 1];
} else if (fabs(F(xi, yj, zkt))
67
< er) {
< er) {
< er) {
< er) {
ht = hz;
ut = u1[i][j][k + 1];
} else {
ht = T(xi, yj, zk) - zk;
ut = 0;
}
/* UNDER */
if (F(xi, yj, zku) >= er) {
hu = hz;
uu = u1[i][j][k - 1];
} else if (fabs(F(xi, yj, zk)) < er) {
hu = hz;
uu = u1[i][j][k - 1];
} else {
hu = zk - U(xi, yj, zk);
uu = 0;
}
u2[i][j][k] =
u1[i][j][k] + tau * psi(xi, yj, zk, nfunc)
+ d(tau) * d(C) * ((ue - u1[i][j][k]) / he (u1[i][j][k] - uw) / hw) / (he +
hw)
+ d(tau) * d(C) * ((un - u1[i][j][k]) / hn (u1[i][j][k] - us) / hs) / (hn +
hs)
+ d(tau) * d(C) * ((ut - u1[i][j][k]) / ht (u1[i][j][k] - uu) / hu) / (ht +
hu);
} else {
u2[i][j][k] = 0;
}
}
}
}
/* n = 1(t = τ) のグラフを描く */
g_cls();
sprintf(message, "t = %g", tau);
g_sel_text(0);
g_text(50.0, 30.0, message);
for (i = 0; i <= Nx; i++)
for (j = 0; j <= Ny; j++)
68
u[i][j] = u2[i][j][Nz / 2];
/* 座標軸を描く */
g_sel_line(1);
g_move(-0.1, 0.0);
g_plot(1.1, 0.0);
g_move(0.0, -0.1);
g_plot(0.0, 2.0);
/* 断面図を描く */
g_sel_line(0);
g_move(xmin, u[0][0]);
for (i = 1; i <= Nx; i++)
g_plot(xmin + i * hx, u[i][Ny / 2]);
/* 鳥瞰図を描く */
g_hidden2(5.0, 5.0, 1.0, -1.0, 1.0, 15.0, 25.0, 30.0, 20.0, 60.0,
100.0, 80.0, u, Nx + 1, Ny + 1, 1, G_SIDE_NONE, 4, 4);
printf("u2 まで描けた \n");
LMax = rint(Tmax / tau);
for (n = 2; n <= LMax; n++) {
/**********************領域***********************/
for (i = 0; i <= Nx; i++) {
xi = xmin + i * hx;
xiw = xmin + (i - 1) * hx;
xie = xmin + (i + 1) * hx;
for (j = 0; j <= Ny; j++) {
yj = ymin + j * hy;
yjn = ymin + (j + 1) * hy;
yjs = ymin + (j - 1) * hy;
for (k = 0; k <= Nz; k++) {
zk = zmin + k * hz;
zkt = zmin + (k + 1) * hz;
zku = zmin + (k - 1) * hz;
if (F(xi, yj, zk) >= er) {
/* WEST */
if (F(xiw, yj, zk) >= er) {
hw = hx;
uw = u2[i - 1][j][k];
} else if (fabs(F(xiw, yj, zk)) < er) {
69
hw = hx;
uw = u2[i - 1][j][k];
} else {
hw = xi - W(xi, yj, zk);
uw = 0;
}
/* EAST */
if (F(xie, yj, zk) >= er) {
he = hx;
ue = u2[i + 1][j][k];
} else if (fabs(F(xie, yj, zk)) < er) {
he = hx;
ue = u2[i + 1][j][k];
} else {
he = E(xi, yj, zk) - xi;
ue = 0;
}
/* NORTH */
if (F(xi, yjn, zk) >= er) {
hn = hy;
un = u2[i][j + 1][k];
} else if (fabs(F(xi, yjn, zk)) < er) {
hn = hy;
un = u2[i][j + 1][k];
} else {
hn = N(xi, yj, zk) - yj;
un = 0;
}
/* SOUTH */
if (F(xi, yjs, zk) >= er) {
hs = hy;
us = u2[i][j - 1][k];
} else if (fabs(F(xi, yjs, zk)) < er) {
hs = hy;
us = u2[i][j - 1][k];
} else {
hs = yj - S(xi, yj, zk);
us = 0;
}
/* TOP */
if (F(xi, yj, zkt) >= er) {
70
ht = hz;
ut = u2[i][j][k + 1];
} else if (fabs(F(xi, yj, zkt)) < er) {
ht = hz;
ut = u2[i][j][k + 1];
} else {
ht = T(xi, yj, zk) - zk;
ut = 0;
}
/* UNDER */
if (F(xi, yj, zku) >= er) {
hu = hz;
uu = u2[i][j][k - 1];
} else if (fabs(F(xi, yj, zku)) < er) {
hu = hz;
uu = u2[i][j][k - 1];
} else {
hu = zk - U(xi, yj, zk);
uu = 0;
}
u3[i][j][k] = 2.0 * u2[i][j][k] - u1[i][j][k]
+
2.0 * d(tau) * d(C) * ((ue - u2[i][j][k]) /
he - (u2[i][j][k] uw) / hw) / (he +
hw)
+
2.0 * d(tau) * d(C) * ((un - u2[i][j][k]) /
hn - (u2[i][j][k] us) / hs) / (hn +
hs)
+
2.0 * d(tau) * d(C) * ((ut - u2[i][j][k]) /
ht - (u2[i][j][k] uu) / hu) / (ht +
hu);
} else {
u3[i][j][k] = 0;
}
}
71
}
}
/* 配列交換 */
for (i = 0; i <= Nx; i++)
for (j = 0; j <= Ny; j++)
for (k = 0; k <= Nz; k++) {
u1[i][j][k] = u2[i][j][k];
u2[i][j][k] = u3[i][j][k];
}
t = tau * n;
if (n % skip == 0) {
for (i = 0; i <= Nx; i++)
for (j = 0; j <= Ny; j++)
u[i][j] = u2[i][j][Nz / 2];
g_cls();
g_line_color(G_BLACK);
sprintf(str, "t = %lf", t);
g_text(50, 30, str);
g_sel_line(0);
g_move(xmin, u[0][0]);
for (i = 1; i <= Nx; i++)
g_plot(xmin + i * hx, u[i][Ny / 2]);
g_line_color(G_BLUE);
for (k = -10; k <= 0; k++)
g_contln2(xmin, xmax, ymin, ymax, u, Nx + 1, Ny + 1, k * 0.2);
g_line_color(G_RED);
for (k = 0; k <= 10; k++)
g_contln2(xmin, xmax, ymin, ymax, u, Nx + 1, Ny + 1, k * 0.2);
g_hidden2(5.0, 5.0, 1.0, -1.0, 1.0, 15.0, 25.0, 30.0, 20.0,
60.0, 100.0, 80.0, u, Nx + 1, Ny + 1, 1, G_SIDE_NONE,
4, 4);
}
}
/* マウスでクリックされるのを待つ */
g_sleep(-1.0);
/* ウィンドウを消す */
g_term();
return 0;
}
72
さいごに
本論文を作成するにあたり、丁寧かつ熱心なご指導をして下さった桂田先生に深く感謝いた
します。今回の論文は C 言語の作成がなかなか進まず、自分自身の知識不足を実感しました
が、多くの助言やご指導のおかげで完成させることが出来ました。また卒業研究発表会では発
表練習、話し方のアドバイスやスライドの添削など前日まで付き合って頂きました。そして発
表で使用したムービーに関しては、数値を変えて良い結果が出るまで何度も試行錯誤して下
さったおかげで、発表当日は自信を持って話すことが出来ました。11 月から本格的に研究を
始め、約 4ヶ月間熱心にご指導して頂き、本当にありがとうございました。
73
参考文献
[1]「音のエピソード」, 音響システムエンジニアリング社の WWW サイト
http://www13.plala.or.jp/sound-planner/episo-do.html
[2] 関口洋正,「3 次元波動方程式の差分法による数値シミュレーション」
2001 年度桂田研卒業研究レポート
http://nalab.mind.meiji.ac.jp/~mk/labo/report/pdf/2001-sekiguchi.pdf
[3] 浜勇樹、「S-W 近似による楕円領域での波動方程式のシミュレーション」
2011 年度桂田研卒業研究レポート
http://nalab.mind.meiji.ac.jp/~mk/labo/report/pdf/2011-hama.pdf
[4]wikipedia, Courant-F riedrichs-Lewy condition(2014 年 3 月 3 日) http://en.wikipedia.
org/wiki/Courant-Friedrichs-Lewy_condition
[5] 川端康司、ヒュームにおける想像と想起のちがいについて.
人文論究,38(4),pp 17-31,1989 年
kgur.kwansei.ac.jp/dspace/bitstream/10236/5539/1/384-02.pdf
74

Download Report