Shortley-Weller 近似による Poisson 方程式のシミュレーション
明治大学総合数理学部現象数理学科
船見 俊哉
2018
年2
月15
日目 次
第1章 はじめに 2
1.1 Shortley-Weller
近似とは. . . . 2
1.2
研究の狙い. . . . 2
第2章 S-W近似 3
2.1
言葉の定義. . . . 3
2.2
差分格子. . . . 3
2.3
不等間隔格子点. . . . 3
2.4
等間隔格子点の領域内外の判断. . . . 4
2.5
不等間隔格子点の座標と格子点間の距離. . . . 4
2.6 Laplacian
のS-W
近似. . . . 5
第3章 円盤領域におけるPoisson方程式のS-W近似 6
3.1
領域と方程式. . . . 6
3.2
差分方程式. . . . 7
3.3
未知ベクトルU
の構成法. . . . 7
3.4
係数行列A
の構成法. . . . 7
3.5
ベクトルF
の構成法. . . . 8
3.6
シミュレーション結果. . . . 9
3.7
まとめ. . . . 10
付 録A 使用したプログラムのソースコード 11
参考文献 22
第 1 章 はじめに
1.1 Shortley-Weller 近似とは
Shortley-Weller
近似(
以下S-W
近似と呼ぶ)
は、G. H. Shortley
とR. Weller
が考えだした差分法 の一種である[1]
。あとで定義する「不等間隔格子点」を用いて行う差分法で、通常の差分法に比べ ていくらか複雑になっている。その一方、局所誤差がO(h)
であっても差分解の誤差はO(h
2)
で、通 常の差分法と同じになるという利点がある[2][3]
。1.2 研究の狙い
この研究では、
S-W
近似によるPoisson
方程式の数値計算プログラムの作成を目標とした。先行研 究ではS-W
近似による熱方程式や波動方程式の陽解法プログラムは存在したが、陰解法プログラム は見つからなかった。先程述べたような利点がありながらも研究が進展しないのは、あくまでも想像 だが、プログラムの作り方がわかりにくいせいかもしれない。ひとまずプログラムを作ることで、次 の研究の一助になれば良いと考えている。またプログラムが完成すれば、山本
[2]
の主張の確認ができるだろう。さらに熱方程式や波動方程 式の陰解法プログラムによる数値実験も行うことができると考えている。第 2 章 S-W 近似
2.1 言葉の定義
定義 2.1
2次元領域を縦横に平行に区切る線分を格子軸と呼び、縦の格子軸と横の格子軸の交点を格子点 と呼ぶ。特に格子軸同士の幅が等間隔なものを等間隔格子と呼び、同様に等間隔格子点と呼ぶ。
2.2 差分格子
R
2の有界領域Ω
上で問題を解くことを考える。このとき、Ω
の閉包をΩ
とし、Ω
を含む長方形領 域領域D
を以下で定義する。
定義 2.2
D := (x
min, x
max)
×(y
min, y
max)
⊃Ω
このようにして定義した
D
を等間隔格子で分割する。N
x, N
y をそれぞれx, y
方向のD
の分割数、h
x, h
yをそれぞれx, y
方向のD
の分割幅とするとき、h
x= x
max−x
minN
x, h
y= y
max−y
minN
y が成り立つ。さらにx
i, y
jをx
i:= x
min+ i
×h
x(0
≤i
≤N
x, i
∈N) y
j:= y
min+ j
×h
y(0
≤j
≤N
y, j
∈N)
と定義すると、
(x
i, y
j)
は等間隔格子点の座標であり、これを(i, j)
の格子点と呼ぶことにする。2.3 不等間隔格子点
定義 2.3
Ω
の境界∂Ω
と格子軸との交点を不等間隔格子点という。
Ω
内の等間隔格子点P = (x
i, y
j)
に対して、その左右上下の等間隔格子点をP
W, P
E, P
N, P
Sと置 く(West, East, North, South
のつもり)
。このとき、∂Ω
付近の点P
に関して、P
W, P
E, P
N, P
Sのい ずれかまたは複数がΩ
に含まれない点P
が存在する。このような場合、図2.1
のように∂Ω
上に不等 間隔格子点を取り、P
W, P
E, P
N, P
Sと置き換える。P
WP
P
SP
EP
N図
2.1:
不等間隔格子点の取り方2.4 等間隔格子点の領域内外の判断
2.3
節で説明した手法を達成するためには、等間隔格子点がΩ
の内部、外部、または境界上のいず れに属するかを判断する必要がある。特に(i, j)
の格子点に対して、次のような関数M
を定義する。
定義 2.4
M(i, j) =
1 ((x
i, y
j)
∈Ω) 0 ((x
i, y
j)
∈∂Ω)
−
1 ((x
i, y
j)
∈D
\Ω)
2.5 不等間隔格子点の座標と格子点間の距離
不等間隔格子点における未知関数の値は境界条件を使うことにする。境界条件を
g(x, y)
とすると、不等間隔格子点の座標が必要になる。
定義 2.5
•
E(x
i, y
j)
は、(E(x
i, y
j), y
j)
が(x
i, y
j)
の右側にある最も近い不等間隔格子点のx
座標で ある。•
W (x
i, y
j)
は、(W (x
i, y
j), y
j)
が(x
i, y
j)
の左側にある最も近い不等間隔格子点のx
座標で ある。•
N (x
i, y
j)
は、(x
i, N (x
i, y
j))
が(x
i, y
j)
の上側にある最も近い不等間隔格子点のy
座標で ある。•
S(x
i, y
j)
は、(x
i, S(x
i, y
j))
が(x
i, y
j)
の下側にある最も近い不等間隔格子点のy
座標である。
以上を用いて、
P = (x
i, y
j)
と周りの4
点との距離h
E, h
W, h
N, h
Sを求める。h
E=
{h
x(P
E ∈Ω)
|
E(x
i, y
j)
−x
i|(else) , h
W=
{h
x(P
W ∈Ω)
|
W (x
i, y
j)
−x
i|(else) h
N=
{
h
y(P
N ∈Ω)
|
N (x
i, y
j)
−y
j|(else) , h
S=
{h
y(P
S ∈Ω)
|
S(x
i, y
j)
−y
j|(else)
2.6 Laplacian の S-W 近似
例として
Laplacian : ∆u(P ) = u
xx(P ) + u
yy(P)
に対してS-W
近似を適用すると、以下のように なる。∆u(P)
≃2
h
E+ h
W(
u(P
E)
−u(P )
h
E −u(P )
−u(P
W) h
W)
+ 2
h
N+ h
S(
u(P
N)
−u(P )
h
N −u(P )
−u(P
S) h
S)
第 3 章 円盤領域における Poisson 方程式の S-W 近似
3.1 領域と方程式
Ω =
{(x, y)
∈R2 |x
2+ y
2< 1
}とし、 {
−
∆u(x, y) = f(x, y) ((x, y)
∈Ω)
u(x, y) = g(x, y) ((x, y)
∈∂Ω) (3.1)
を考える。
2.4
節で定義したM
は、この問題では次のようになる。M(i, j) =
1 (x
2i+ y
j2< 1) 0 (x
2i+ y
j2= 1)
−
1 (x
2i+ y
j2> 1)
数学的にはこれで良いが、丸め誤差のため、数値計算において浮動小数点数の完全な一致を判定す るプログラムは書くことができない。そこで
∂Ω
にε
程度の幅を持たせて、∂Ω
上ではないが極めて 近い等間隔格子点は∂Ω
上にあるという判断をする。例えば、図3.1
では、点P
は、半径1
の円周= ∂Ω(
赤線)
上にはない。しかし、半径1
±ε
の円周=
青線の間に存在しているため、このようなとき は∂Ω
上にあるという判断をする。double
型の丸め誤差は10
−14ほどであるため、特にε = 10
−13と し、先程のM
を次のように書き換える。M (i, j) =
1 (1
−(x
2i+ y
j2) > ε) 0 (|1
−(x
2i+ y
2j)| ≤ ε)
−
1 (1
−(x
2i+ y
j2) <
−ε)
P
図
3.1:
境界線の幅また、
2.5
節で定義したE, W, N, S
は、円盤領域において次のようになる。E(x
i, y
j) =
√
1
−y
j2, W (x
i, y
j) =
−√1
−y
j2, N (x
i, y
j) =
√
1
−x
2i, S(x
i, y
j) =
−√1
−x
2i3.2 差分方程式
Poisson
方程式のS-W
近似による差分化は次のようになる。{ −(
2 hE+hW
(u(PE)−u(P)
hE −u(P)−hWu(PW))
+
h 2N+hS
(u(PN)−u(P)
hN − u(P)−hSu(PS)))
= f (P ) (P
∈Ω)
u(P ) = g(P ) (P
∈∂Ω)
(3.2)
このとき、u(P ) (
∀P
∈Ω)
が未知数となる。続いて、この差分方程式を領域内および境界上の全ての等間隔格子点
P
について並べた連立方程 式を構成する。そのような連立方程式を行列の書き方に沿ってAU = F
とするが、ここでA
は係数 行列、U
は未知ベクトル、F
は既知ベクトルであることを断っておく。3.3 未知ベクトル U の構成法
まずは、領域内部および境界上の等間隔格子点に
0
番目から番号をつけていく。例えば、N
x= N
y= 4
のときは図3.2
のように番号をつけていく。すなわち最も下の格子軸のうち左にある等間隔格子点か ら順番に番号をつけ、その格子軸上の領域内部または境界上の等間隔格子点が無くなったら次の格子 軸に移るということを繰り返して行う。0 12
1 2 3
4 5 6 7 8
9 10 11
図
3.2:
番号のつけ方後で番号
k
からその点の(i, j)
を求める必要があるため、次の定義を行う。
定義 3.1
k
番目の点をP
kと表すことにする。P
kが(i, j)
の格子点であるとき、p(k) = i, q(k) = j, V (i, j) = k
となるp(k), q(k), V (i, j)
を定義する。
このようにして数え上げた等間隔格子点の個数を
n
とする。ここでu(P
k) = u
kと表すことにする(0
≤k
≤n
−1)
。3.2
節でも述べたようにu
kは未知数であるから、U =
t(u
0, u
1,
· · ·, u
n−1)
として未知ベクトルU
を 構成する。3.4 係数行列 A の構成法
係数行列
A
∈M (n,
R)
は次のようにして構成する。u
kの係数、すなわち対角成分A
kkは、A
kk=
2(h
Eh
W+ h
Nh
S) h
Eh
Wh
Nh
S(P
k ∈Ω)
1 (P
k ∈∂Ω)
である。
P
E, P
W, P
N, P
Sの番号をk
E, k
W, k
N, k
Sとすると、A
kkE=
−
2
h
E(h
E+ h
W) (P
E ∈Ω)
0 (else)
, A
kkW=
−
2
h
W(h
E+ h
W) (P
W ∈Ω)
0 (else)
A
kkN=
−
2
h
N(h
N+ h
S) (P
N ∈Ω)
0 (else)
, A
kkS=
−
2
h
S(h
N+ h
S) (P
S ∈Ω)
0 (else)
それ以外の
A
の成分は全て0
である。3.5 ベクトル F の構成法
天下り的だが、まず
n
次元縦ベクトルF
′=
t(F
0′, F
1′,
· · ·, F
n′−1)
を以下のようにして定義する。F
k′=
{f (P
k) (P
k ∈Ω) g(P
k) (P
k ∈∂Ω)
S-W
近似では領域外の等間隔格子点の代わりに境界上の不等間隔格子点を使うが、そのときの未知 関数の値は境界条件を使う。そのため実は既知関数である。そこで境界条件を使うようなときは、不 等間隔格子点における未知関数の値=
境界条件の値とその係数を、右辺に移項してから計算を行う。そのために次のような
n
次元縦ベクトルF
E(
右補正ベクトルと呼ぶことにする)
を定義する。F
Ek=
2
h
E(h
E+ h
W) (M (p(k) + 1, q(k)) =
−1)0 (else)
また左、上、下補正ベクトル
F
W, F
N, F
Sを同様に定義する。F
W k=
2
h
W(h
E+ h
W) (M (p(k)
−1, q(k)) =
−1)0 (else)
F
N k=
2
h
N(h
N+ h
S) (M (p(k), q(k) + 1) =
−1)0 (else)
F
Sk=
2
h
S(h
N+ h
S) (M(p(k), q(k)
−1) =
−1)
0 (else)
以上の
4
つの補正ベクトルとF
′を用いてF = F
′+ F
E+ F
W+ F
N+ F
Sと定義する。
3.6 シミュレーション結果
以上のようにして構成した連立方程式
AU = F
をU
について解き、その結果を可視化した。f
とg
はf (x, y) =
−4(x
2+ y
2−1)e
−x2−y2g(x, y) = e
−1としている。
図
3.3: N
x= N
y= 10
図3.4: N
x= N
y= 20
図
3.5: N
x= N
y= 40
図3.6: N
x= N
y= 80
図
3.7: N
x= N
y= 160
図
3.8:
誤差の推移次に厳密解
u
と数値解U
の誤差を調べた。ここでいう誤差とは、各等間隔格子点P
における|u(P)
−U (P)|
の最大値のことを指す。先程のf
とg
に対して厳密解u
はu(x, y) = e
−x2−y2 となる。図
3.8
が分割数を増やしたときの誤差の推移である。横軸が分割数、縦軸が誤差max
P∈Ω{
u(P )
−U (P)
} である。3.7 まとめ
図
3.8
からわかるように、山本[2]
の主張を確認することができた。また今回作成したプログラムをたたき台にして、熱方程式や波動方程式など様々な偏微分方程式に 対して陰解法による数値計算プログラムが作成できると考えている。
