定常 Stokes 方程式の有限要素解の事後誤差評価と事前誤差評価

(1)

定常 Stokes 方程式の有限要素解の事後誤差評価と事前誤差評価

明治大学大学院

理工学研究科基礎理工学専攻数学系福澤誠人

指導教員桂田祐史助教授

2005 年 2 月 25 日

(2)

を具体的に評価するものだが、不動点定理の条件の成立を数値計算で確認することで、非線形方程式の解の存在を証明することもできる。このように精度保証付き数値計算には二つの意味があるが、今回は前者の意味でのみの精度保証付き数値計算を行う。

従来この精度保証付き数値計算を行うにはたくさんの困難さがあった。当初は区間演算を行うことさえ難しかったが、最近のパーソナルコンピュータやワークステーションでは

IEEE754

規格に準拠した

CPU

が採用されているので、区間演算に必要な丸めのコントロールを比較的簡単にできるようになり、区間演算は容易に実現できる。また基本的アルゴリズムの研究が進み、それに伴って手軽に用いることができるソフトウェアの研究が進んだことも大きい。今回は

MATLAB

で区間演算を実現するツールボックスである

S.M.Rump

による

INT- LAB(http://www.ti3.tu-harburg.de/~rump/intlab/)

を用いたが、これを用いると

MAT- LAB

の通常の計算と同様に簡単に区間演算による計算を行うことができる。

この論文では、中尾充宏、山本野人、渡部善隆による

Stokes

問題に対する有限要素解の精度保証の評価法

[1]

を元に、事後誤差評価式

(3.10)

と事前誤差評価式

(4.5), (4.8)

に現れる右辺を評価する数値実験の追試を行い、さらに区間演算を用いて厳密な意味での精度保証付き数値実験を行うことを目標にして取り組んだ。

この論文の概略は次のようなものである。まず第

1

節で扱う問題と関数空間を設定し、第

2

節で弱解の存在を保証する

inf-sup

条件に関連した定数の数値評価を与える。第

3

節で事後誤差評価式を導き、第

4

節では事前誤差評価式とその構成法を述べる。第

5

節では

Aubin-Nitsche’s

trick

に似た手法を用いて流速に関する

L

² ノルムでの事後誤差評価式と事前誤差評価式を与

える。第

6

節では事前誤差評価をする際に必要となる一般化固有値の絶対値の最大の精度保証付きでの評価法について述べる。第

7

節と第

8

節で事前誤差評価と事後誤差評価それぞれの具体的な計算方法と実験結果について述べる。

1.1 Stokes 方程式

次の凸多角形領域

Ω ⊂ R

²における

Stokes

方程式の

(Dirichlet)

境界値問題を考える。

 



 

−ν4u + ∇p = f in Ω, div u = 0 in Ω,

u = 0 on ∂Ω.

(1.1)

ここで、ν >

0

は粘性係数

(正定数)、u = [u

₁

, u

₂

]

^T は二次元速度ベクトル、f

= [f

₁

, f

₂

]

^T は単位質量に働く外力密度で滑らかな関数、

p

は圧力である。

(5)

1.2 関数空間

H

^k

(Ω)

を

Ω

上での

k

階の

Sobolev

空間、

(·, ·)

を

L

²

(Ω)

での内積、

(∇·, ∇·)

を

H

₀¹

(Ω)

での内積とする。以下のように関数空間を定義する。

H

₀¹

(Ω) := ©

v ∈ H

¹

(Ω) ; v = 0 on ∂ Ω ª , L

²₀

(Ω) := ©

v ∈ L

²

(Ω) ; (v, 1) = 0 ª , S := H

₀¹

(Ω)

²

× L

²₀

(Ω).

L

²

(Ω), H

₀¹

(Ω)

のノルムはそれぞれ

|q|

₀

:= (q, q)

¹²

,

|v|

₁

:= |∇v|

₀ とする。また

H

²

(Ω)

のセミノルムとして

|u|

2

:=

Ã¯ ¯

¯ ¯ ∂

²

u

∂x

²

¯ ¯

2 0

+ 2

¯ ¯

¯ ¯ ∂

²

u

∂x∂y

¯ ¯

2 0

+

¯ ¯

¯ ¯ ∂

²

u

∂y

²

¯ ¯

2 0

!

¹

2

を用いる。そして関数空間

V, V

^⊥を以下のように定義する。

V := ©

v ∈ H

₀¹

(Ω)

²

; div v = 0 ª , V

^⊥

:= ©

v ∈ H

₀¹

(Ω)

²

; (∇v, ∇w) = 0, w ∈ V ª .

最後に、H⁻¹

(Ω)

²を

H

₀¹

(Ω)

²の双対空間とする。

2 inf-sup 条件の定数の評価

この節では

Stokes

方程式の境界値問題

(1.1)

を弱形式で書き直し、弱解の存在を保証する

inf-sup

条件に関連した定数の数値評価を与える。この定数を用いて、

S

のそれぞれの要素に

対してノルム不等式を与える。

2.1 変分法的定式化

Stokes

方程式の境界値問題

(1.1)

に対する弱形式として

Find [u, p] ∈ H

₀¹

(Ω)

²

× L

²

(Ω) s.t.

ν(∇u, ∇v) − (p, div v) = (f, v), ∀v ∈ H

₀¹

(Ω)

²

, (2.1)

−(div u, q) = 0, ∀q ∈ L

²

(Ω)

が得られる。次のように記号の定義を行う。

 

 



 

 

X = H

₀¹

(Ω)

²

, Y = L

²

(Ω), f ∈ L

²

(Ω)

²

,

a(u, v) := ν(∇u, ∇v ) (u, v ∈ V ),

F (v) := (f, v) (v ∈ V ),

Bv := − div v(∈ W ) (v ∈ V ),

g := 0(∈ R(B)) (R(B)

は

B

の値域).

(2.2)

(6)

この定式化から次の定理を用いれば

(2.1)

の解

[u, p]

の存在と一意性を得る

(ただし p

については定数関数分の不定性がある

)

。

¶ ³

Theorem 2.1

X, Y

は実

Hilbert

空間、

a(·, ·) : X×X → R

は有界対称な双

1

次形式で非負、

F (·) : X → R

は有界線形汎関数、B

: X → Y

は有界線形作用素で

R(B)

は

Y

の閉集合、gは

R(B)

の与えられた元とする。さらに

a(·, ·)

は

N (B)

で強圧的

(ここで N(B )

は

B

の零空間)、すなわち正定数

µ

が存在して次式が成立すると仮定する。

a(v, v) ≥ µkvk

²_V

, (∀v ∈ N (B)).

このとき

Find [u, p] ∈ X × Y s.t.

a(u, v) + (Bv, p)

Y

= F (v), ∀v ∈ X, (Bu, q)

_Y

= g, ∀q ∈ Y

の解は存在し、その

1

つを

{u, p} ∈ X × Y

と記すと、uは

X

で一意に定まり、pは

R(B)

内に限れば一意であり、

Y

では

N (B

^∗

)

分だけの付加項の不定性がある、つまり

p

⁰

, p

^∗

∈ Y

を条件を満たす

p ∈ Y

とすると、

p

⁰

− p

^∗

∈ N (B

^∗

)

となる

(

ここで

(·, ·)

_Y は

Y

での内積、また

B

^∗は

B

の共役作用素)。

µ ´

この

Theorem

の証明は菊地

[12]

に載っている。

2

(2.2)

の定式化において

N (B

^∗

)

は定数関数である。なぜなら、q

∈ Y , v ∈ X

に対し

(B

^∗

q, v)

_X

= (q, Bv)

_Y

= −(q, div v)

に注意し、さらに

v ∈ C

₀^∞

(Ω)

²

(

この集合は

X = H

₀¹

(Ω)

² で稠密

)

ととれば、

q ∈ N (B

^∗

)

すなわち

B

^∗

q = 0

となるための

q

に対する条件は、∇q

= 0

となるからである。

Theorem 2.1

において

F

は

X

上の有界線形汎関数であること、

g

は

R(B)

に含まれることのみを要請する。さらに

N (B

^∗

)

は

(2.2)

の定式化において定数関数であるので、

(2.1)

の解の一意存在は次の定理に拡張することができる。

¶ ³

Theorem 2.2

S ∈ H

⁻¹

(Ω)

²

, R ∈ (L

²₀

(Ω))

⁰ とする。このとき、

Find [u, p] ∈ H

₀¹

(Ω)

²

× L

²₀

(Ω) s.t.

ν(∇u, ∇v) − (p, div v ) = Sv, ∀v ∈ H

₀¹

(Ω)

²

, (2.3)

−(div u, q) = Rq, ∀q ∈ L

²₀

(Ω)

の解

[u, p] ∈ H

₀¹

(Ω)

²

× L

²₀

(Ω)

は一意に存在する。

µ ´

次に

S × S

上の双線形形式

L

を次のように定める。

L ([u, p], [v, q]) := ν(∇u, ∇v) − (p, div v ) − (q, div u), [u, p], [v, q] ∈ S . (2.4)

(7)

このとき、(1.1)の変分法的定式化を次で与える。

Find [u, p] ∈ S s.t. L ([u, p], [v, q]) = (f, v), ∀[v, q] ∈ S . (2.5) (2.3)

において

Sv = (f, v) , R = 0

とした問題と

(2.5)

の同値性については明らかであるので、

(2.5)

は

Theorem 2.2

を用いると、S で一意解を持つことがわかる。

また次に述べる定理により、このとき

Ω

のみに依存する正定数

β

_cが存在して

[u, p]

inf

∈S, [u, p]6= 0

sup

[v, q]∈S, [v, q]6= 0

L ([u, p], [v, q])

(|u|

1

+ |p|

0

)(|v |

1

+ |q|

0

) ≥ β

_c

(2.6)

が成り立つ。

¶ ³

Theorem 2.3

X, Y : Hilbert

空間、

a(·, ·) : X × Y → R

を連続双線形形式とする。このとき方程式

Find x ∈ X s.t. a(x, y) = F y, ∀y ∈ Y (2.7)

が任意の

F ∈ Y

⁰に対して一意解を持つためには

a

が次の

2

つの条件を満たすことが必要十分である。

C

₁

:= inf

x∈X

sup

y∈Y

a(x, y)

kxk

_X

kyk

_Y

> 0, (2.8)

sup

x∈X

a(x, y) kxk

X

> 0, ∀y ∈ Y, y 6= 0. (2.9)

さらに

(2.7)

の解

x ∈ X

は、

kxk

_X

≤ kF k

_Y⁰

C

₁ を満たす。

µ ´

この

Theorem

の証明は土屋

[14]

この

Theorem

を適用して

(2.6)

をみたす正定数

β

cの存在を証明するには、

X = Y = S , Y

⁰

= S

⁰とした上で、

∀F ∈ S

⁰

, ∃!S ∈ H

₀⁻¹

(Ω)

²

, ∃!R ∈ (L

²₀

(Ω))

⁰

s.t. F [v, q] = Sv + Rq ([v, q] ∈ S )

となることに注意して

Theorem 2.2

を用いればよい。

(2.6)

から

sup

[v, q]∈S, [v, q]6= 0

L ([u, p], [v, q])

|v |

₁

+ |q|

₀

≥ β

_c

(|u|

₁

+ |p|

₀

), ∀[u, p] ∈ S (2.10)

を得る。[u, p]

∈ S

に対して

δ(u, p) = sup

[v, q]∈S, [v, q]6= 0

L ([u, p], [v, q])

|v|

1

+ |q|

0

(2.11)

(8)

と定義する。この節の残りで

δ(u, p)

を用いて

|u|

₁と

|p|

₀を評価する

(Theorem 2.4)。そのため

に、次の

Lemma 2.1

を用いる。

¶ ³

Lemma 2.1 (Babuˇ ska-Aziz

の不等式

)

Ω

を

R

²の有界

Lipschitz

領域とする。このとき、

∃β > 0, ∀q ∈ L

²₀

(Ω), ∃!v ∈ V

^⊥

s.t.

div v = q, (2.12)

|v|

1

≤ 1

β |q|

0

. (2.13)

ただし、

β > 0

は

Ω

に依存する定数である。

µ ´

この

Lemma

の証明は

Roger Temam [20]

を見よ。

Ω

が滑らかな境界を持つ場合の証明は付録に載せる。

Lemma 2.1

から

inf

p∈L²₀(Ω), p6= 0

sup

u∈H¹₀(Ω)², u6= 0

−(p, div v)

|u|

₁

|p|

₀

≥ β (2.14)

を得る。

(2.14)

は

(2.5)

が

S

で一意解を持つことを保証する

S

の

inf-sup

条件と呼ばれる条件である。もし

Ω

が一般の有界連結領域ならば、βの上界や下界を評価することは難しい。しかしながら、

Ω

が星型の場合、この定数

β

は次の

Horgan

の

Lemma 2.2

により数値的に決定

できる。¶ ³

Lemma 2.2 (Horgan [7])

2

次元の有界

Lipschitz

領域

Ω

は原点に関して星型であり、境界は次式によって極座標で表されるとする。

r = f(θ) on ∂Ω.

また、

F (θ) :=

 



"

1 +

µ f

⁰

(θ) f (θ)

¶

₂

#

¹

2

+ |f

⁰

(θ)|

f(θ)

 



2

とする。このとき

1 β = q

1 + max

θ

F (θ)

とおくと

Lemma 2.1

の不等式

(2.13)

が成り立つ。

µ ´

(9)

Ω = (−1, −1) × (1, 1)

の場合、

f(θ) =

 

 

 

 

 1 cos θ

³

− π

4 ≤ θ ≤ π 4

´ , 1

sin θ µ π

4 ≤ θ ≤ 3π 4

¶ ,

− 1 cos θ

µ 3π

4 ≤ θ ≤ 5π 4

¶ ,

− 1 sin θ

µ 5π

4 ≤ θ < 7π 4

¶

である。このとき、

F (θ) =

 

 



 

 

−1 + 2 1 + sin θ

µ

− π

4 ≤ θ ≤ 0, π ≤ θ ≤ 5π 4

¶ ,

−1 + 2 1 − sin θ

µ

0 ≤ θ ≤ π 4 , 3π

4 ≤ θ ≤ π

¶ ,

−1 + 2 1 − cos θ

µ π

4 ≤ θ ≤ π 2 , 3π

2 ≤ θ < 7π 4

¶ ,

−1 + 2 1 + cos θ

µ π

2 ≤ θ ≤ 3π 4 , 5π

4 ≤ θ ≤ 3π 2

¶

であり、

θ = − π 4 , π

4 , 3π 4 , 5π

4

のとき

F (θ)

は最大値

3 + 2 √

2

をとる。よって、

q 1 + max

θ

F (θ) = q

4 + 2 √ 2

となる。この結果は

Ω

が一般の正方形領域の場合にも成り立つ。

2.2 ノルム不等式

Lemma 2.1

の定数

β

を用いて、次の不等式を得る。

¶ ³

Theorem 2.4

∀[u, p] ∈ S

に対して

δ(u, p)

を

(2.11)

によって定める。このとき、次の評価が成り立つ。

|u|

₁

≤ µ 1

ν

²

+ 1 β

²

¶

¹

2

δ(u, p),

|p|

₀

≤ µ 1

β + ν β

²

¶

δ(u, p).

(2.15)

µ ´

(10)

証明

V

は

H

₀¹

(Ω)

²の閉部分空間であるので、

H

₀¹

(Ω)

²

= V ⊕ V

^⊥ である。よって、

∀u ∈ H

₀¹

(Ω)

²

, ∃!w ∈ V, ∃!u

₀

∈ V

^⊥

, u = w + u

₀

. w 6= 0

とする。

δ(u, p)

の定義式

(2.11)

で

v = w, q = 0

とすると、

δ(u, p) ≥ ν(∇u, ∇w) − (p, div w)

|w|

₁

= ν(∇w, ∇w)

|w|

₁

≥ ν|w|

₁

(2.16)

を得る。この式は

w = 0

でも成り立つ。

一方、u₀

6= 0

であるとき

δ(u, p)

の定義式で

v = 0

とすると、

δ(u, p) ≥ −(q, div u)

|q|

₀

, 0 6= ∀q ∈ L

²₀

(Ω).

よって、q

= − div u

₀とおくと

q ∈ L

²₀

(Ω)

で

β|u

₀

| ≤ |q|

₀ が成り立つ。このとき、

δ(u, p) ≥ |q|

₀

≥ β|u

₀

|

₁

(2.17)

となる。この式は

u

₀

= 0

のときも成り立つ。よって、(2.16), (2.17)を用いると

|u|

²₁

= (∇(w + u

₀

), ∇(w + u

₀

)) = (∇w, ∇w) + (∇u

₀

, ∇u

₀

) = |w|

²₁

+ |u

₀

|

²₁

≤ µ 1

ν

²

+ 1 β

²

¶

δ(u, p)

²

を得る。

(2.15)

の

1

つ目の不等式を示すことができた。

次に

(2.15)

の

2

つ目の不等式を示す。∀p

∈ L

²₀

(Ω)

に対し、Lemma 2.1より

v ∈ V

^⊥ として

div v = −p, |v|

₁

≤ 1

β |p|

₀

(2.18)

を満たすようにとれる。p

6= 0

とする。(∇u,

∇v) = (∇u

₀

, ∇v)

より、

δ(u, p) ≥ ν(∇u

₀

, ∇v) − (q, div u

₀

) + |p|

²₀

|v|

1

+ |q|

0

, ∀q ∈ L

²₀

(Ω).

まず

u

₀

6= 0

のとき、

q ∈ L

²₀

(Ω)

を以下のように定める。

q := K div u

₀

, K := ν(∇u

₀

, ∇v)

| div u

0

|

²₀

.

上式より、

ν(∇u

₀

, ∇v) − (q, div u

₀

) = ν(∇u

₀

, ∇v) − K(div u

₀

, div u

₀

) = 0. (2.19)

さらに、(2.13)より

|u

₀

|

₁

≤ 1

β | div u

₀

|

₀

. (2.20)

(11)

したがって、

|q|

0

= K| div u

0

|

0

= ν(∇u

₀

, ∇v)

| div u

₀

|

₀

≤ ν|u

₀

|

₁

|v|

₁

| div u

₀

|

₀

.

よって、

(2.18), (2.20)

より、

|q|

₀

≤ ν

β

²

|p|

₀

. (2.21)

(2.19), (2.21)

は

q = 0

ととることにより

u

₀

= 0

の場合も成り立つ。ゆえに、

δ(u, p) ≥ ν(∇u

₀

, ∇v) − (q, div u

₀

) + |p|

²₀

|v |

₁

+ |q|

₁

≥ |p|

²₀

1

β |p|

₀

+ ν β

²

|p|

₀

= µ 1

β + ν β

²

¶

₋₁

|p|

₀

.

上式は

p = 0

のときも成り立つ。よって

(2.15)

の

2

つ目の不等式も示すことができた。

2

3 事後誤差評価

この節では流速と圧力を近似する有限要素空間を導入する。そして

Theorem 2.4

を用いて

Stokes

方程式の事後誤差限界を示す。

3.1 有限要素部分空間

T

hを

Ω ⊂ R

²のスケールパラメータ

h > 0

に依存する三角形または四角形からなる分割族とする。また、X_h

⊂ H

₀¹

(Ω) ∩ C(Ω)

を流速

u

を近似する有限要素部分空間、Y_h

⊂ L

²₀

(Ω) ∩ C(Ω) ∩ H

¹

(Ω)

を圧力

p

を近似する有限要素部分空間、

X

_h^∗を

X

_hの基底関数と

∂Ω

上の節点に対応する基底関数で張られる空間、S_h

:= X

_h²とする。

注意中尾、山本、渡辺

[1]

では

Y

_h

⊂ H

¹

(Ω)

と書かれていなかったが、暗に仮定されていて、数値実験で採用された

Y

_hはこの性質を満たしている。

(2.5)

の有限要素解は次で定義される。

Find [u

_h

, p

_h

] ∈ S

_h

× Y

_h

s.t. L ([u

_h

, p

_h

], [v

_h

, q

_h

]) = (f, v

_h

), ∀[v

_h

, q

_h

] ∈ S

_h

× Y

_h

. (3.1)

山本、中尾

[8]

で提案された事後処理の手順を示す。

X

_h

( X

_h^∗

⊂ H

¹

(Ω), X

_h

6= X

_h^∗

であり

X

_h^∗は

H

¹

(Ω)

の要素を近似する能力を持つ。また、

L

²

-

射影

P

0

: L

²

(Ω) → X

h

, L

²

-

射影

P ˆ

₀

: L

²

(Ω) → X

_h^∗

, H

₀¹

-射影 P

₁

: H

₀¹

(Ω) → X

_h を以下のように定める。

(v − P

₀

v, φ) = 0, ∀φ ∈ X

_h

, (v − P ˆ

₀

v, φ) = 0, ∀φ ∈ X

_h^∗

, (∇(v − P

₁

v), ∇φ) = 0, ∀φ ∈ X

_h

.

また、

w

h

∈ X

hに対して、

∇w

h

∈ (X

_h^∗

)

²と

4w

h

∈ L

²

(Ω)

を

∇w

h

:=

µ P ˆ

0

∂w

_h

∂x , P ˆ

0

∂w

_h

∂y

¶

,

4w

h

:= div ∇w

h

(12)

により定める。∀v_h

∈ S

_hに対して、

(−4v

_h

, φ) = (∇v

_h

, ∇φ), ∀φ ∈ H

₀¹

(Ω)

²

, (3.2)

|∇v

_h

− ∇v

_h

|

₀

= inf

wh∈(x^∗_h)²×(x^∗_h)²

|w

_h

− ∇v

_h

|

₀

(3.3)

が成り立つ。以下、上式を示す。まず、(3.2)については、v_h

= [v

_h1

, v

_h2

] ∈ S

_h

, φ = [φ

₁

, φ

₂

] ∈ H

₀¹

(Ω)

² として

(−4v

h

, φ) = (−4v

h1

, φ

1

) + (−4v

h2

, φ

2

)

= (− div ∇v

_h1

, φ

₁

) + (− div ∇v

_h2

, φ

₂

)

= (∇v

_h1

, ∇φ

₁

) + (∇v

_h2

, ∇φ

₂

)

= (∇v

_h

, ∇φ).

よって

(3.2)

は示せた。次に

(3.3)

については、

v

h

= [v

h1

, v

h2

] ∈ S

hとして

|∇v

_h

− ∇v

_h

|

²₀

= |∇v

_h1

− ∇v

_h1

|

²₀

+ |∇v

_h2

− ∇v

_h2

|

²₀

.

ここで、

|∇v

_h1

− ∇v

_h

|

²₀

=

¯ ¯

¯ ¯ P ˆ

₀

∂v

h1

∂x − ∂v

h1

∂x

¯ ¯

2 0

+

¯ ¯

¯ ¯ P ˆ

₀

∂v

h1

∂y − ∂v

h1

∂y

¯ ¯

2

0

= inf

wh1x∈X_h^∗

¯ ¯

¯ ¯ w

h1x

− ∂v

_h1

∂x

¯ ¯

2 0

+ inf

wh1y∈X_h^∗

¯ ¯

¯ ¯ w

h1y

− ∂v

_h1

∂y

¯ ¯

2

0

= inf

wh1∈(X_h^∗)²

|w

_h1

− ∇v

_h1

|

²₀

.

同様に、

|∇v

h2

− ∇v

h2

|

²₀

= inf

wh2∈(X_h^∗)²

|w

h2

− ∇v

h2

|

²₀

.

よって、

|∇v

h

− ∇v

h

|

²₀

= inf

wh∈(X_h^∗)²×(X_h^∗)²

|w

h

− ∇v

h

|

²₀

.

上式で正の平方根をとれば

(3.3)

を得る。

2

ここで、

X

hに次の仮定をおく。

¶ ³

Assumption 3.1

ξ∈X

inf

h

|v − ξ|

1

≤ C

0

h|v|

2

, ∀v ∈ H

₀¹

(Ω) ∩ H

²

(Ω) (3.4)

をみたす

v, h

に依存しない正定数

C

₀が具体的に評価可能である。

µ ´

この仮定は多くの有限要素空間で成り立つ。

X

hが

1

次元の区分

1

次要素の空間では

C

0

= 1

ととれる。また

1

次元の区分

2

次要素のテンソル積として定義される

2

次元矩形要素では、

π C

₀

= 1

2π

とれる

(

中尾、山本、木村

[9])

。射影

P

₁については、

Assumption 3.1

と

Aubin-Nitche’s

trick

から、

(13)

任意の

v ∈ H

₀¹

(Ω)

に対して、

|v − P

₁

v|

₁

≤ |v|

₁

, (3.5)

|v − P

₁

v|

₀

≤ C

₀

h|v|

₁

(3.6)

が成り立つ。

(3.5)

はピタゴラスの定理より明らかに成り立つ。以下

(3.6)

を示す。

v ∈ H

₀¹

(Ω)

とする。ここで

g = v − P

₁

v

とおく。g

∈ L

²

(Ω)

である。次の問題を考える。

Find ψ ∈ H

₀¹

(Ω) s.t. (∇ψ, ∇ϕ) = (g, ϕ), ∀ϕ ∈ H

₀¹

(Ω).

Riez

の定理から、この問題の解

ψ

は存在する。ここで、

|g|²₀

= (g, g) = (g, v

−P₁v) = (∇ψ,∇(v−P₁v)) = (∇(ψ−P₁ψ),∇(v−P₁v))≤ |ψ−P₁ψ|₁|v−P₁v|₁.

Assumption 3.1

と

(3.6)

より

|g|

²₀

≤ C

₀

h|ψ|

₂

|v|

₁ となる。次の

Lemma 3.1

を用いると、

|g|

0

≤ C

0

h|v|

1

を得る。2

¶ ³

Lemma 3.1

Ω

が区分的に滑らかな有界凸領域であり、

(

−4v = g in Ω, v = 0 on ∂Ω

の解

v

が

v ∈ H

²

(Ω)

を満たすとする。このとき、

|v|

₂

≤ |g|

₀ が成り立つ。

µ ´

Ω

が

2

次元領域のときの証明は中尾、山本

[16]

の

p.99

補題

5.2

3.2 事後評価

[u, p], [u

_h

, p

_h

]

をそれぞれ

(2.5), (3.1)

の解とする。また、

e

_h

= u − u

_h

, ε

_h

= p − p

_h とする。まず

δ(e

h

, ε

h

)

の上界を評価する。

任意の

[v, q] ∈ S

に対して、

L ([e

_h

, ε

_h

], [v, q]) = ν(∇e

_h

, ∇v) − (ε

_h

, div v) − (q, div e

_h

) (3.7)

(14)

である。また、(2.5)で

[v, q] = [ξ

_h

, q

_h

]

とした式から

(3.1)

で

[v

_h

, q

_h

] = [ξ

_h

, q

_h

]

とした式を引くと、任意の

[ξ

_h

, q

_h

] ∈ S

_h

× Y

_hに対して、

ν(∇e

_h

, ∇ξ

_h

) − (ε

_h

, div ξ

_h

) − (q

_h

, div e

_h

) = 0 (3.8)

を得る。

(3.8)

で

q

_h

= 0

ととると、任意の

ξ ∈ S

_hに対して、

ν(∇e

h

, ∇ξ

h

) − (ε, div ξ

h

) = 0. (3.9) (3.7)

から

(3.9)

を引くと、

L ([e

_h

, ε

_h

], [v, q]) = ν(∇e

_h

, ∇(v − ξ

_h

)) − (ε

_h

, div(v − ξ

_h

)) − (q, div e

_h

)

= ν(∇(u − u

_h

), ∇(v − ξ

_h

)) − (p − p

_h

, div(v − ξ

_h

)) + (q, div u

_h

)

= ν(∇u

_h

− ∇u

_h

), ∇(v − ξ

_h

)) + ν(∇u − ∇u

_h

, ∇(v − ξ

_h

))

− (p − p

_h

, div(v − ξ

_h

)) + (q, div u

_h

).

ここで、

(2.5)

より

ν(∇u, ∇v) − (p, div v) = (f, v), ∀v ∈ H

₀¹

(Ω)

² となることを用いると、

L ([e

_h

, ε

_h

], [v, q]) = ν(∇u

_h

− ∇u

_h

, ∇(v − ξ

_h

)) + (f + ν4u

_h

− ∇p

_h

, v − ξ

_h

) + (q, div u

_h

)

≤ ν|∇u

h

− ∇u

h

|

0

|v − ξ

h

|

1

+ |ν4u

h

− ∇p

h

+ f|

0

|v − ξ

h

|

0

+ | div u

h

|

0

|q|

0

.

ここで、ξ_h

∈ S

_hを

v = (v

₁

, v

₂

)

の

H

₀¹

-射影であるとする。すなわち、ξ

_h

= (P

₁

v

₁

, P

₁

v

₂

)

^T とする。すると、

(3.5), (3.6)

より、

L ([e

_h

, ε

_h

], [v, q]) ≤ ν|∇u

_h

− ∇u

_h

|

₀

|v |

₁

+ |ν4u

_h

− ∇p

_h

+ f |

₀

C

₀

h|v|

₁

+ | div u

_h

|

₀

|q|

₀

≤ (ν|∇u

h

− ∇u

h

|

0

+ C

0

h|ν4u

h

− ∇p

h

+ f|

0

+ | div u

h

|

0

)(|v|

1

+ |q|

0

).

このようにして次の結果を得る。

¶ ³

Lemma 3.2

[u, p], [u

_h

, p

_h

]

をそれぞれ

(2.5), (3.1)

の解として

e

_h

= u − u

_h

, ε

_h

= p − p

_hとする。さらに

C

0を

Assumption 3.1

の定数とする。このとき任意の

0

でない

[v, q] ∈ S

に対して次式が成り立つ。

L ([e

h

, ε

h

], [v, q])

|v|

₁

+ |q|

₀

≤ ν|∇u

_h

− ∇u

_h

|

₀

+ C

₀

h|ν4u

_h

− ∇p

_h

+ f |

₀

+ | div u

_h

|

₀

.

µ ´

Theorem 2.4

と

Lemma 3.2

から次に述べる

Stokes

方程式の有限要素解の事後誤差限界を得る。

(15)

¶ ³

Theorem 3.1 (

事後誤差限界

)

[u, p], [u

_h

, p

_h

]

をそれぞれ

(2.5), (3.1)

の解とする。また

β

は

(2.13)

を満たす定数とする。

さらに

C

₀を

Assumption 3.1

の定数とする。このとき次の不等式が成り立つ。

|u − u

_h

|

₁

≤ µ 1

ν

²

+ 1 β

²

¶

¹

2

C(u

_h

, p

_h

),

|p − p

_h

|

₀

≤ µ 1

β + ν β

²

¶

C(u

_h

, p

_h

).

(3.10)

ここで、C(u_h

, p

_h

)

は有限要素解を用いて次のように定められる量である。

C(u

_h

, p

_h

) := ν|∇u

_h

− ∇u

_h

|

₀

+ C

₀

h|ν4u

_h

− ∇p

_h

+ f |

₀

+ | div u

_h

|

₀

. (3.11)

µ ´

証明

Theorem 2.4

において

u

を

u − u

h

= e

h

, p

を

p − p

h

= ε

h として用いると、

|u − u

_h

|

₁

≤ µ 1

ν

²

+ 1 β

²

¶

¹

2

δ(e

_h

, ε

_h

),

|p − p

_h

|

₀

≤ µ 1

β + ν β

²

¶

δ(e

_h

, ε

_h

)

を得る。また、

Lemma 3.2

を用いると、

δ(e

_h

, ε

_h

) := sup

[v, q]∈S, [v, q]6= 0

L ([e

_h

, ε

_h

], [v, q])

|v|

1

+ |q|

0

≤ C(u

_h

, p

_h

)

を得る。

2

4 構成的アプリオリ誤差評価

前節では

Stokes

問題の有限要素解の事後誤差限界について述べた。本節では前節と同様の

方法を用いて、

2

種類の構成的アプリオリ誤差の限界の導出方法を考え、アプリオリ定数の評価の計算手順を述べる。

4.1 アプリオリ評価

最初の方法は、前節で述べた事後誤差限界の結果を元にした手法である。任意の

f = [f

₁

, f

₂

] ∈ L

²

(Ω)

²に対して、

P

0

f ∈ S

hを

P

₀

f = [P

₀

f

₁

, P

₀

f

₂

]

^T によって定義する。

L

²

-

射影の性質から、

|f − P

₀

f |

²₀

= |f |

²₀

− |P

₀

f |

²₀

.

(16)

したがって、

0 ≤ θ ≤ π

2

を満たすある

θ

が存在して、

|P

₀

f |

₀

= |f |

₀

sin θ,

|f − P

0

f |

0

= |f |

0

cos θ (4.1)

となる。ここで、

f ∈ L

²

(Ω)

²によらない定数

K

₁

,K

₂

,K

₃が存在して次を満たすと仮定する。

|∇u

_h

− ∇u

_h

|

₀

≤ K

₁

|P

₀

f |

₀

, (4.2)

|ν4u

_h

− ∇p

_h

+ P

₀

f|

₀

≤ K

₂

|P

₀

f |

₀

, (4.3)

| div u

h

|

0

≤ K

3

|P

0

f |

0

. (4.4)

これらの定数は、一般化固有値問題を数値計算を用いて解くことにより決定される。詳しくは次の小節以降で述べる。このとき次の定理を得る。

¶ ³

Theorem 4.1 (

アプリオリ誤差限界

I)

[u, p], [u

_h

, p

_h

]

をそれぞれ

(2.5), (3.1)

β

は

(2.13)

をみたす定数とし、C₀ を

Assumption 3.1

の定数とする。さらに

K

₁

, K

₂

, K

₃を

(4.2), (4.3), (4.4)

をみたす定数とする。このとき任意の

f ∈ L

²

(Ω)

²に対して次式が成り立つ。

|u − u

h

|

1

≤ µ 1

ν

²

+ 1 β

²

¶

¹

2

C

1

(h)|f |

0

,

|p − p

h

|

0

≤ µ 1

β + ν β

²

¶

C

1

(h)|f |

0

.

(4.5)

ここで、

C

1

(h) := p

(νK

1

+ C

0

hK

2

+ K

3

)

²

+ (C

0

h)

²

(4.6)

である。

µ ´

証明

任意の

f ∈ L

²

(Ω)

²に対して、

(3.11)

より

C(u

_h

, p

_h

) = ν|∇u

_h

− ∇u

_h

|

₀

+ C

₀

h|ν4u

_h

− ∇p

_h

+ P

₀

f + f − P

₀

f |

₀

+ | div u

_h

|

₀

.

ここで、

(4.2), (4.3), (4.4)

を用いると、

C(u

h

, p

h

) ≤ νK

1

|P

0

f |

0

+ C

0

h(K

2

|P

0

f |

0

+ |f − P

0

f |

0

) + K

3

|P

0

f |

0

.

さらに

(4.1)

を用いると、

C(u

h

, p

h

) ≤ ((νK

1

+ C

0

hK

2

+ K

3

) sin θ + C

0

h cos θ)|f|

0

≤ p

(νK

₁

+ C

₀

hK

₂

+ K

₃

)

²

+ (C

₀

h)

²

|f |

₀

= C

₁

(h)|f |

₀ を得る。

2

(17)

2

つ目の方法は、∇u_hと

4u

_hを用いることなく直接に導く方法である。(3.7)から

(3.9)

を引くと、任意の

[v, q] ∈ S

と任意の

ξ

_h

∈ S

_hに対して、

L ([e

_h

, ε

_h

], [v, q]) = ν(∇(u − u

_h

), ∇(v − ξ

_h

)) − (p − p

_h

, div(v − ξ

_h

)) + (q, div u

_h

)

を得る。ここで、

ξ

_h

:= v

_h

= [P

₁

v

₁

, P

₁

v

₂

]

^T

, (v = [v

₁

, v

₂

]

^T

)

ととる。すると

H

₀¹

-射影の性質より、

L ([e

h

, ε

h

], [v, q]) = ν(∇u, ∇(v − v

h

)) − (p − p

h

, div(v − v

h

)) + (q, div u

h

).

ここで、(2.5)より

ν(∇u, ∇v) − (p, div v) = (f, v), ∀v ∈ H

₀¹

(Ω)

² となることを用いると、

L ([e

_h

, ε

_h

], [v, q]) = (f − ∇p

_h

, v − v

_h

) + (q, div u

_h

)

≤ |f − ∇p

_h

|

₀

|v − v

_h

|

₀

+ |q|

₀

| div u

_h

|

₀

.

さらに

(3.6)

を用いると、

L ([e

_h

, ε

_h

], [v, q]) ≤ |f − ∇p

_h

|

₀

C

₀

h|v|

₁

+ |q|

₀

| div u

_h

|

₀

≤ (C

₀

h|f − ∇p

_h

|

₀

+ | div u

_h

|

₀

)(|v|

₁

+ |q|

₀

).

ゆえに次の

Lemma

が成り立つ。

¶ ³

Lemma 4.1

[u, p], [u

_h

, p

_h

]

をそれぞれ

(2.5), (3.1)

の解として

e

_h

= u − u

_h

, ε

_h

= p − p

_hとする。さらに

C

₀を

Assumption 3.1

の定数とする。このとき

0

でない任意の

[v, q] ∈ S

に対して次式が成り立つ。

L ([e

h

, ε

h

], [v, q])

|v|

₁

+ |q|

₀

≤ C

0

h|f − ∇p

h

|

0

+ | div u

h

|

0

.

µ ´

ゆえに、Lemma 4.1の結果から、Theorem 3.1の

C(u

_h

, p

_h

)

を以下のようにとることもできる。

C(u

_h

, p

_h

) = C

₀

h|f − ∇p

_h

|

₀

+ | div u

_h

|

₀

.

ここで、

f ∈ L

²

(Ω)

²によらない定数

K

4が存在して次を満たすと仮定する。

| − ∇p

_h

+ P

₀

f |

₀

≤ K

₄

|P

₀

f |

₀

. (4.7)

K

₄を決定するための方法は後で述べる。これからもう

1

つのアプリオリ誤差評価式を得る。

(18)

Theorem 4.2 (

アプリオリ誤差限界

II)

[u, p], [u

_h

, p

_h

]

をそれぞれ

(2.5), (3.1)

β

は

(2.13)

をみたす定数とし、C₀ を

Assumption 3.1

の定数とする。さらに

K

₃

, K

₄を

(4.4), (4.7)

をみたす定数とする。このとき任意の

f ∈ L

²

(Ω)

²に対して次式が成り立つ。

|u − u

h

|

1

≤ µ 1

ν

²

+ 1 β

²

¶

¹

2

C

2

(h)|f |

0

,

|p − p

h

|

0

≤ µ 1

β + ν β

²

¶

C

2

(h)|f |

0

.

(4.8)

ここで、

C

2

(h) := p

(C

0

hK

4

+ K

3

)

²

+ (C

0

h)

²

(4.9)

である。

µ ´

証明

任意の

f ∈ L

²

(Ω)

²に対して、

C(u

h

, p

h

) = C

0

h|f − ∇p

h

|

0

+ | div u

h

|

0

= C

0

h| − ∇p

h

+ P

0

f + f − P

0

f |

0

+ | div u

h

|

0

≤ C

₀

h| − ∇p

_h

+ P

₀

f | + C

₀

h|f − P

₀

f |

₀

+ | div u

_h

|

₀

.

ここで、

(4.4), (4.7)

を用いると、

C(u

_h

, p

_h

) ≤ C

₀

h(K

₄

|P

₀

f |

₀

+ |f − P

₀

f |

₀

) + K

₃

|P

₀

f|.

さらに

(4.1)

を用いると、

C(u

h

, p

h

) ≤ ((C

0

hK

4

+ K

3

) sin θ + C

0

h cos θ)|f |

0

≤ p

(C

₀

hK

₄

+ K

₃

)

²

+ (C

₀

h)

²

|f|

₀

= C(h)|f |

₀ を得る。

2

4.2 定数 C

₁

(h), C

₂

(h) の計算

この節では

(4.6), (4.9)

の定数

C

₁

(h), C

₂

(h)

に現れる

K

₁

, K

₂

, K

₃

, K

₄の具体的な計算方法を述べる。

dim X

_h

= n, dim X

_h^∗

= ˆ n, dim Y

_h

= m

とする。X_h

( X

_h^∗ であるから

n > n ˆ

である。

{φ

_j

}

_1≤j≤nを

X

_hの基底関数、

{ φ ˆ

_j

}

_1≤j≤ˆ_nを

X

_h^∗の基底関数、

{ψ

_j

}

_1≤j≤mを

Y

_hの基底関数とする。

このとき、実係数

{a

⁽¹⁾_j

}

_1≤j≤n

, {a

⁽²⁾_j

}

_1≤j≤n

, {b

_j

}

_1≤j≤m を用いて有限要素解

u

_h

= [u

⁽¹⁾_h

, u

⁽²⁾_h

]

^T

∈ S

_h

, p

_h

∈ Y

_hは次の形に一意に表される。

u

⁽¹⁾_h

= X

n

i=1

a

⁽¹⁾_i

φ

_i

, u

⁽²⁾_h

= X

n

i=1

a

⁽²⁾_i

φ

_i

, p

_h

= X

m

i=1

b

_i

ψ

_i

.

(19)

このとき、X_h

, Y

_hの基底関数を用いて、(3.1)を書き直す。(3.1)より

ν(∇u

⁽¹⁾_h

, ∇v

_h⁽¹⁾

) −

Ã

p

_h

, ∂v

_h⁽¹⁾

∂x

!

= (f

₁

, v

_h⁽¹⁾

), ∀v

_h⁽¹⁾

∈ X

_h

, ν(∇u

⁽²⁾_h

, ∇v

_h⁽²⁾

) −

Ã

p

_h

, ∂v

_h⁽²⁾

∂y

!

= (f

₂

, v

_h⁽²⁾

), ∀v

_h⁽²⁾

∈ X

_h

,

− Ã

q

_h

, ∂u

⁽¹⁾_h

∂x

!

− Ã

q

_h

, ∂u

⁽²⁾_h

∂y

!

= 0, ∀q

_h

∈ Y

_h

となるので、

ν X

n

i=1

a

⁽¹⁾_i

(∇φ

_i

, ∇φ

_j

) − X

m

i=1

b

_i

µ

ψ

_i

, �

定常 Stokes 方程式の有限要素解の 事後誤差評価と事前誤差評価