θ 法の安定性の解析 - 円盤領域における差分法 56

第 9 章円盤領域における差分法 56

9.9 θ 法の安定性の解析

b[j] -= a[j][i] * b[i];

}

b[n-1] /= a[n-1][n-1];

for (i = n - 2; i >= 0; i--) {

mm = i + m; if (mm > n-1) mm = n-1;

for (j = i + 1; j <= mm; j++) b[i] -= a[i][j] * b[j];

b[i] /= a[i][i];

} }

と変形できることに注意すればよい。

B ≥ 0

とする。任意の

i

について

∑

N j=1

b

_ij

≤ 1

が成り立つことから、

∥ Bv ∥

_∞

≤ ∥ v ∥

_∞

(v ∈ R

)

が得られる。

(

長方形の場合の繰り返しになるが

)

実際

∑

N j=1

b

_ij

v

≤

∑

N j=1

| b

_ij

| | v

| ≤ max

1≤j≤N

| v

|

∑

N j=1

| b

_ij

| = ∥ v ∥

_∞

∑

N j=1

b

_ij

≤ ∥ v ∥

_∞ の

i

についての最大値を取ればよい。

次に

A

⁻¹

v

∞

≤ ∥ v ∥

_∞

(v ∈ R

)

を証明する。θ

= 0

のときは

A = I

であるから明らかなので、0

< θ ≤ 1

とする。w

:= A

⁻¹

v

とおくと

v := Aw

である。番号

i

に対応する格子点の上下左右の格子点がすべて領域内部に属している場合、

v

= [

1 + 2θ (

λ

+ λ

_φ

r

²_i

)]

w

− θλ

_φ

r

_i²

(w

_i+1

+ w

_i₋₁

) − θλ

[(

1 + 1 2i

)

w

_i+p

+ (

1 − 1 2i

) w

_i₋_p

]

が成り立っているので、

| v

| ≥ [

1 + 2θ (

λ

+ λ

_φ

r

²_i

)]

| w

| − θλ

_φ

r

²_i

( | w

_i+1

| + | w

_i₋₁

| ) − θλ

[(

1 + 1 2i

)

| w

_i+p

| + (

1 − 1 2i

)

| w

_i₋_p

| ]

≥ [

1 + 2θ (

λ

+ λ

_φ

r

²_i

)]

| w

| − θλ

_φ

r

²_i

( ∥ w ∥

_∞

+ ∥ w ∥

_∞

) − θλ

[(

1 + 1 2i

)

∥ w ∥

_∞

+ (

1 − 1 2i

)

∥ w ∥

_∞

]

= [

1 + 2θ (

λ

+ λ

_φ

r

_i²

)]

| w

| − 2θ (

λ

+ λ

_φ

r

²_i

)

∥ w ∥

_∞

.

この不等式はすべての

i

について成立する。これからすべての

i

について、

∥ v ∥

_∞

≥ [

1 + 2θ (

λ

+ λ

_φ

r

_i²

)]

| w

| − 2θ (

λ

+ λ

_φ

r

²_i

)

∥ w ∥

_∞

.

ノルム

∥ · ∥

_∞ の定義から、

∥ w ∥

_∞

= | w

_i₀

|

となる

i

₀ が存在する。

i = i

₀ について上の不等式を用いると

∥ v ∥

_∞

≥ [

1 + 2θ (

λ

+ λ

_φ

r

²_i

)]

∥ w ∥

_∞

− 2θ (

λ

+ λ

_φ

r

_i²

)

∥ w ∥

_∞

= ∥ w ∥

_∞

.

すなわち

A

⁻¹

v

∞

= ∥ w ∥

_∞

≤ ∥ v ∥

_∞

.

これから次の定理を得る。

定理

9.9.2 (θ

法の

∥ · ∥

_∞ 安定性解析

) θ = 1

または

(

0 ≤ θ < 1

かつ

τ ≤ 1

2(1 − θ) · h

²_r

h

²_φ

1 + h

²_r

)

を満たすとき、任意の

n

について

U

ⁿ⁺¹

∞

≤ ∥ U

ⁿ

∥

_∞

.

証明

∥ U

ⁿ⁺¹

∥

_∞

= ∥ A

⁻¹

BU

ⁿ

∥

_∞

≤ ∥ BU

ⁿ

∥

_∞

≤ ∥ U

ⁿ

∥

_∞

.

余談

9.9.1 ∥ A

⁻¹

v ∥

_∞ について、結局は長方形領域の場合と同様の評価が得られるわけだが、

証明は少しだけ変更されている。実は、長方形領域

(

というよりも

1

次元開区間

)

の場合の

θ

法の安定性の、筆者が当初思いついた証明は、おそろしく回りくどいものであった

(

印刷されて、かなりの部数配布されてされてしまっていて、一時は穴があったら入りたい気分だった)。

それに比べると、現在の長方形領域の場合の証明は大分すっきりしていると思う。円盤の場合に、その証明はそのままでは通用しないことが分かって、一瞬「困った」と思ったが、ふと考えてみたら、少しだけ先祖返りすると、つまり長方形領域の場合の証明の某バージョンは、そのまま円盤の場合の証明に使えることに気が付いた。世の中捨てたものではない。

さて、上の定理で書いた安定性の十分条件であるが、θ <

1

の場合は、受け入れがたいほど

τ

を小さく取ることになる。残念ながら、数値実験による状況証拠はこの条件が限りなく必要に近いことを示しているようである

(

早い話、

θ = 1

の完全陰解法でないと使い物にならない、

らしい)。本当は必要性の証明も出来れば良いのだが、いまのところは出来ていない。

第 10 ^章 ^{円盤領域における} Laplacian ^の Shortley-Weller ^近似

Shortley-Weller

近似については、山本

[23]

、

[24]

に紹介されている。これは不等間隔の差分法である。

正直、そういうことが出来るのではないかと昔から漠然と考えいたので、「やはりそうか」

という感想を持っている。

[23], [22]

に紹介されていることはなかなか魅力的で、ぜひ原論文にあたって習得したいものである。

10.1 1 次元楕円型境界値問題に対する Shortley-Weller ^近似

Shortley-Weller

近似について、簡単な

1

次元問題で説明しよう

(

これは山本

[23]

による

)

。未知関数

u = u(x)

に関する次の

2

階常微分方程式の境界値問題を考えよう。

− d dx

(

p(x) du dx

)

+ q(x)u = f (x) (x ∈ (a, b)), (10.1)

u(a) = α, u(b) = β.

(10.2)

ただし、

p, q, f

は与えられた関数で、適当な滑らかさを持ち、区間

[a, b]

上で

p > 0

と仮定する。

この区間

(a, b)

を

a = x

₀

< x

₁

< · · · < x

= b

と分割する。刻み幅

h

:= x

− x

_i₋₁

(i = 1, 2, . . . , n)

は定数とは限らない。

h := max

1≤i≤n

h

,

u

:= u(x

) (i = 0, 1, . . . , n),

x

_i+1/2

:= (x

+ x

_i+1

)/2, p

_i+1/2

:= p(x

_i+1/2

) (i = 0, 1, . . . , n − 1)

とおく。

d dx

(

p(x) du dx

)

x=xi

=

p

_i+1/2

u

_i+1

− u

h

_i+1

− p

_i₋_1/2

u

− u

_i₋₁

h

(h

_i+1

+ h

)/2 + O(h

_i+1

− h

) + O

( h

³_i+1

+ h

³_i

h

_i+1

+ h

)

が成り立つ。左辺を右辺第

1

項で近似することを

Shortley-Weller

近似という。

等間隔差分の場合に

d dx

(

p(x) du dx

)

x=xi

≒ p

_i+1/2

u

i+1

− u

h − p

_i₋_1/2

u

− u

i−1

h

という近似をするのは「普通」、「相場」のようだから、不等間隔刻みを採用するのが

Shortley-Weller

近似の肝、ということだろう。

この問題に対してどのような誤差評価ができるかについては、山本

[23]

を参照せよ。

10.1.1 自力で考えたメモ

f (x + h

) = f(x) + f

^′

(x)h

+ 1

2 f

^′′

(x)h

²_r

+ 1

3! f

⁽³⁾

(x)h

³_r

+ 1

4! f

⁽⁴⁾

(x)h

⁴_r

+ · · · , f (x − h

) = f(x) − f

^′

(x)h

+ 1

2 f

^′′

(x)h

²_l

− 1

3! f

⁽³⁾

(x)h

³_l

+ 1

4! f

⁽⁴⁾

(x)h

⁴_l

+ · · ·

であるから、

f(x + h

) − f (x) h

= f

^′

(x) + 1

2 f

^′′

(x)h

+ 1

3! f

⁽³⁾

(x)h

²_r

+ 1

4! f

⁽⁴⁾

(x)h

³_r

+ · · · , f(x − h

) − f (x)

h

= − f

^′

(x) + 1

2 f

^′′

(x)h

− 1

3! f

⁽³⁾

(x)h

²_l

+ 1

4! f

⁽⁴⁾

(x)h

³_l

+ · · · .

これから

f (x

) − f(x)

h

− f(x) − f(x

)

h

= 1

2 f

^′′

(x) (h

+ h

) + 1 3! f

⁽³⁾

(

h

²_r

− h

²_l

) + 1

4! f

⁽⁴⁾

(

h

³_r

+ h

³_l

)

+ · · · .

ゆえに

f

^′′

(x) = 2 h

+ h

( f(x

) − f (x)

h

− f (x) − f (x

) h

)

− h

3 f

⁽³⁾

(x) − 1

12 f

⁽⁴⁾

(x) h

³_r

+ h

³_l

h

+ h

+ · · · .

これから誤差は

O(h

− h

)

であり、特に

h

= h

であれば、

O(h

)

であることが分る。

10.2 2 ^次元 Laplacian ^の Shortley-Weller ^近似

空間

1

次元の場合には、何か特別な理由がないと不等間隔格子点を用いることはないと思うが、2次元領域の場合は、不等間隔格子点の利用は自然で有効だと思われる。

Ω

を

R

² の有界領域とする。刻み幅

h

= h

の等間隔差分格子を作る。境界付近の格子点

P ( ∈ Ω)

に対しては、その上下左右の格子点が

Ω

の外に飛び出してしまうことがある。

この場合、等間隔であることをあきらめて境界上に「格子点」を定める

(カッコつきにしたの

は格子線の交点でないのに格子点と呼ぶのは、少々変な感じがするからである

)

。

P

の上下左右の格子点を

P

, P

_Rとし

(

気持ちは

up, down, left, right

だが、某君のレポートでは代わりに

north, south, west, east

というのを使った)、

h

:= | P

− P | , h

:= | P

− P | , h

:= | P

− P | , h

:= | P

− P |

とおく。h_U

, h

∈ (0, h]

である。

△ u(P ) ≒

u(P

) − u(P )

h

− u(P) − u(P

) h

(h

+ h

)/2 +

u(P

) − u(P )

h

− u(P) − u(P

) h

(h

+ h

)/2

山本

[23]

には次の

3

つの文献が紹介されている。

1. G. H. Shortley and R. Weller, The numerical solution of Laplace’s equation, J. Appl. Physics,

9, pp.334-344 (1938).

2. J. H. Bramble and B. E. Hubbard, On the formulation of finite diﬀerence analogues of the Dirichlet problem for Poisson’s equation, Numer. Math. 4 (1962), pp.313–327.

3. N. Matsunaga and T. Yamamoto, Superconvergence of the Shortley-Weller approxima-tion for Dirichlet problems, J. Comput. Appl. Math. 116, pp.263–273 (2000).

山本先生はその後も以下の論文を出している¹。

1. Tetsuro Yamamoto, Qing Fang and Xiaojun Chen, Superconvergence and nonsupercon-vergence of the Shortley-Weller approximations for Dirichlet problems, Numer. Funct. Anal.

and Optimiz., 22 (3&4), pp.455–470 (2001).

2. Zi-Cai Li, Hsin-Yun Hu, Qing Fang and Tetsuro Yamamoto, Superconvergence of So-lution Derivatives for the Shortley-Weller Deﬀerence Approximation of Poisson’s Equa-tion. II. Singularity Problems, Numer. Funct. Anal. and Optimiz., 24 (3&4), pp.195–221 (2003).

3. Qing Fang, Yukihiro Shogenji and Tetsuro Yamamoto, Error Analysis of Adaptive Finite Diﬀerence Methods Using Stretching Functions for Polor Coordinate Form of Poisson-Type Equation, Numer. Funct. Anal. and Optimiz., 24 ( 1&2), pp.17–44 (2003).

どうやってプログラミングするか²、学生にとって挑戦的な課題であると思う

(

直接参考にするものがないので)。山本先生が紹介している文献は楕円型の問題のみのようだけれど、発展問題の解析はないのかなあ？

ドキュメント内 II ( (ページ 89-94)

θ 法の安定性の解析

第 9 章 円盤領域における差分法 56

9.9 θ 法の安定性の解析

B ≥ 0

i

∑

b

≤ 1

∥ Bv ∥

≤ ∥ v ∥

(v ∈ R

)

(

)

∑

b

v

≤

∑

| b

| | v

| ≤ max

| v

|

∑

| b

| = ∥ v ∥

∑

b

≤ ∥ v ∥

i

A

v

≤ ∥ v ∥

(v ∈ R

)

= 0

A = I

< θ ≤ 1

:= A

v

v := Aw

i

v

= [

1 + 2θ (

λ

+ λ

r

)]

w

− θλ

r

(w

+ w

) − θλ

[(

1 + 1 2i

)

w

+ (

1 − 1 2i

) w

]

| v

| ≥ [

1 + 2θ (

λ

+ λ

r

)]

| w

| − θλ

r

( | w

| + | w

| ) − θλ

[(

1 + 1 2i

)

第 9 章円盤領域における差分法 56