Explicit a posteriori local error estimation for FEM solutions

TL;DR

This paper introduces an explicit local a posteriori error estimator for FEM solutions of Poisson's equation using the Hypercircle method, effective even without $H^2$ regularity, validated through numerical experiments.

Contribution

It proposes a novel explicit local error estimation technique for FEM solutions of Poisson's equation applicable without $H^2$ regularity assumptions.

Findings

Effective local error estimates demonstrated on 2D and 3D problems.

Method works even when solution lacks $H^2$ regularity.

Numerical experiments confirm the efficiency of the approach.

Abstract

For the finite element solution of Poisson's equation, a local a posteriori error estimation based on the Hypercircle method is proposed. Even for the solution of Poisson's equation without the $H^2$ regularity, this method can provide explicit local error estimation. The efficiency of the proposed method is demonstrated by numerical experiments for the boundary value problem of Poisson's equation defined on the 2D and 3D domains.

Tables

Table 1. Table 1: The relation between mesh size h ℎ h and error estimation

$h$	${\kappa }_h$	$C\left(h\right)$	${\overline{E}}_L$	$Er{r}^{\left(1\right)}$	$Er{r}^{\left(2\right)}$	$Er{r}^{\left(3\right)}$	$E_L$	${\overline{E}}_G$
1/8	0.057	0.070	0.488	0.195	0.337	0.198	0.131	0.546
1/16	0.030	0.036	0.182	0.069	0.123	0.092	0.063	0.264
1/32	0.015	0.018	0.070	0.025	0.044	0.044	0.031	0.129
1/64	0.008	0.009	0.029	0.009	0.016	0.022	0.015	0.064
1/128	0.004	0.005	0.013	0.003	0.006	0.011	0.008	0.032

Table 2. Table 2: The relation between mesh size h ℎ h and error estimation

$h$	${\kappa }_h$	$C\left(h\right)$	${\overline{E}}_L$	$Er{r}^{\left(1\right)}$	$Er{r}^{\left(2\right)}$	$Er{r}^{\left(3\right)}$	${\overline{E}}_G$
1/8	0.073	0.083	1.411	0.494	1.084	0.330	1.800
1/16	0.046	0.050	0.542	0.177	0.425	0.188	0.869
1/32	0.028	0.029	0.208	0.063	0.165	0.088	0.427
1/64	0.017	0.018	0.083	0.022	0.064	0.043	0.211
1/128	0.011	0.011	0.034	0.008	0.025	0.021	0.105

Table 3. Table 3: The condition of local error estimation

Size of $A_M$	Size of $A_C$	${\kappa }_h$	$C_{0}h$	$C\left(h\right)$
1794250	234571	0.044	0.053	0.069

Table 4. Table 4: Local error estimation

${\overline{E}}_L$	$Er{r}^{\left(2\right)}$	$Er{r}^{\left(3\right)}$	${\overline{E}}_G$
0.284	0.080	0.012	0.341

Equations

$$H(\text{mbox}div;\Omega ):=\{\mathbf{q}\in (L^2(\Omega ){)}^2;\text{mbox}div\mathbf{q}\in L^2(\Omega )\}$$

$$-\Delta u=f\text{mbox}in\Omega ,\frac{\partial u}{\partial \mathbf{n}}=g_N\text{mbox}on{\Gamma }_N,u=g_D\text{mbox}on{\Gamma }_D$$

$$\text{mbox}Findu\in V\text{mbox}s.t.(\nabla u,\nabla v)=(f,v)+(g_N,v{)}_{{\Gamma }_N}\forall v\in V_0$$

$$V:=\{v\in H^1(\Omega );v=g_D\text{mbox}on{\Gamma }_D\},V_0:=\{v\in H^1(\Omega );v=0\text{mbox}on{\Gamma }_D\}$$

$$V_0:=\{v\in H^1(\Omega );{\int }_{\Omega }vdx=0\}$$

$$(g_N,v{)}_{{\Gamma }_N}={\int }_{{\Gamma }_N}{g}_{N}vds.$$

$$\Vert \varphi {\Vert }_{\Omega }\le C_p\Vert \nabla \varphi {\Vert }_{\Omega }\forall v\in V_0$$

$$\text{mbox}Find(\lambda ,\varphi )\in \mathbb{R}\times V_0\text{mbox}s.t.(\nabla \varphi ,\nabla v)=\lambda (\varphi ,v)\forall v\in V_0$$

$$\text{mbox}Find({\lambda }_h,{\varphi }_h)\in \mathbb{R}\times V_0^h\text{mbox}s.t.(\nabla {\varphi }_h,\nabla v_h)=\lambda ({\varphi }_h,v_h)\forall v_h\in V_0^h$$

$${\lambda }_1\ge \frac{{\lambda }_{h,1}}{1+(0.1893h{)}^2{\lambda }_{h,1}}$$

$$h:=\underset{K\in {\mathcal{T}}_h}{\max }{h}_K$$

$$V^h$$

$$V_0^h$$

$$\text{mbox}Find{u}_h\in V^h\text{mbox}s.t.(\nabla u_h,\nabla v_h)=(f,v_h)+(g_N,v_h{)}_{{\Gamma }_N}\forall v_h\in V_0^h$$

$$X^h:=\{v\in L^2(\Omega ):v\text{mbox}は各三角形要素上で定数\}$$

$$X^h:=\{v\in L^2(\Omega ):v\text{mbox}は各三角形要素上で定数,{\int }_{\Omega }vdx=0\}$$

$$W^h:=\{{\mathbf{p}}_h\in H(\text{mbox}div;\Omega ):{\mathbf{p}}_h=(a_K+c_{K}x,b_K+c_{K}y)\text{mbox}inK\in {\mathcal{T}}_h\}$$

$$\text{mbox}Find({\mathbf{p}}_h,{\mu }_h)\in W^h\times X^h,{\mathbf{p}}_h\cdot n=g_N\text{mbox}on{\Gamma }_N\text{mbox}s.t.$$

$$({\mathbf{p}}_h,{\mathbf{q}}_h)+(\text{mbox}div{\mathbf{q}}_h,{\mu }_h)+(\text{mbox}div{\mathbf{p}}_h,{\eta }_h)=-(f,{\eta }_h)+(g_D,({\mathbf{q}}_h\cdot n){)}_{{\Gamma }_D}\forall ({\mathbf{q}}_h,{\eta }_h)\in W^h\times X^h$$

$$(f-{\pi }_{h}f,{\eta }_h)=0\forall {\eta }_h\in X^h$$

$$\Vert f-{\pi }_{h}f{\Vert }_{\Omega }\le C_{0}h|f{|}_{H^1(\Omega )}\forall f\in H^1(\Omega )$$

$$-\Delta \varphi =f\text{mbox}in\Omega ,\frac{\partial \varphi }{\partial \mathbf{n}}=0\text{mbox}on{\Gamma }_N,\varphi =0\text{mbox}on{\Gamma }_D$$

$$\text{mbox}Find\varphi \in V_0$$

$$\text{mbox}Find{\varphi }_h\in V_0^h$$

$$(\nabla (v-P_{h}v),\nabla v_h)=0\forall v_h\in V_0^h$$

$$\Vert \nabla P_h\varphi {\Vert }_{\Omega }\le \Vert \nabla \varphi {\Vert }_{\Omega }\le C_p\Vert f{\Vert }_{\Omega }$$

$$\text{mbox}Find({\tilde{\mathbf{p}}}_h,{\sigma }_h)\in W^h\times X^h,{\tilde{\mathbf{p}}}_h\cdot n=0\text{mbox}on{\Gamma }_N\text{mbox}s.t.$$

$$({\tilde{\mathbf{p}}}_h,{\mathbf{q}}_h)+(\text{mbox}div{\mathbf{q}}_h,{\sigma }_h)+(\text{mbox}div{\tilde{\mathbf{p}}}_h,{\tau }_h)=-(f,{\tau }_h)\forall ({\mathbf{q}}_h,{\tau }_h)\in W^h\times X^h$$

$${\kappa }_h:=\underset{f_h\in X^h}{\max }\underset{\begin{array}{c}{v}_h\in V_0^h,{\tilde{\mathbf{p}}}_h\in W^h,\\ {\tilde{\mathbf{p}}}_h\cdot n=0,\text{div }{\tilde{\mathbf{p}}}_h+f_h=0\end{array}}{\min }\frac{\Vert \nabla v_h-{\tilde{\mathbf{p}}}_h\Vert }{\Vert f_h\Vert }$$

$$\text{mbox}div\tilde{\mathbf{p}}+f=0,\tilde{\mathbf{p}}\cdot n=0\text{mbox}on{\Gamma }_N$$

$$\Vert \nabla \varphi -\nabla v{\Vert }_{\Omega }^2+\Vert \nabla \varphi -\tilde{\mathbf{p}}{\Vert }_{\Omega }^2=\Vert \nabla v-\tilde{\mathbf{p}}{\Vert }_{\Omega }^2$$

$$(\nabla \overline{\varphi },\nabla v)$$

Full text

本論文ではHypercircle法を用いてPoisson方程式の有限要素解に対する局所事後誤差評価の手法を提案する． 当該手法の特徴として，解の$H^{2}$正則性によらずに関心のある領域における有限要素解の定量的な局所誤差評価を得ることができる． また、2次元および3次元領域で定義されるPoisson方程式に対して数値実験を行って，提案手法の有効性を示した．

1 はじめに

本研究では，Poisson方程式の境界値問題に対して，Hypercircle法による有限要素解の局所事後誤差評価手法を検討する. 研究の背景として，半導体の抵抗率測定法を支配する偏微分方程式について，近似解の局所誤差評価が要求されることが挙げられる． 具体的に，抵抗率測定に使用される4探針法(Fig.1参照)では，測定試料の表面に4本の探針$A,B,C,D$を接触させ，探針$A,D$間に固定電流を流し，探針$B,C$間の電位差を測定することで抵抗率を測定する[16]． このとき，数値手法によって探針周囲での電位（Poisson方程式の境界値問題の解に対応する）の近似値を求め，近似値の局所誤差を評価することが要求される．

有限要素解の大域的な事後誤差評価については既に多くの研究結果があり，特に近年，Hypercircle法を用いた有限要素解の誤差の具体的な値または誤差の上界評価(定量的な誤差評価)を提供する手法が報告される([8, 5, 7, 11, 2]など)． 特に，劉-大石の論文[8]では菊地の事後誤差評価[5]の拡張として事前誤差評価を提案したことがあり，この評価は本研究の提案手法に重要な役割を果たす． 有限要素解の局所的な誤差評価については，誤差の収束オーダーなどの定性的な誤差評価([1, 14]など)が考えられたが，定量的な有限要素解の局所誤差評価は少なかった．

この論文では，Hypercircle法を利用することで，有限要素解の局所事後誤差評価における世界初の手法を提案する．提案手法の特徴として，L字型領域などの非凸な領域における方程式の厳密解が$H^{2}$正則性を持ってない場合でも，関心のある領域におけるシャープな誤差評価が可能である． 抵抗率測定の場合， 提案手法によって非凸形状の測定試料に対しても，探針周囲での電位の誤差評価が可能となる． 論文では，2次元領域と３次元領域の上で定義されるPoisson方程式の境界値問題に対して数値実験を行い，提案手法の有効性を検討した．

2 準備

2.1 問題設定

この論文を通じて，特別な説明がない限り$\Omega$を$\mathbb{R}^{2}$の有界多角形領域とし，Sobolev空間$H^{m}(\Omega)$ $(m>0)$に関する標準的な記号[10]を用いる． $L^{2}(\Omega)$空間のノルムを$\|\cdot\|_{L^{2}(\Omega)}$または$\|\cdot\|_{\Omega}$とかき， 記号$|\cdot|_{H^{m}(\Omega)},\|\cdot\|_{H^{m}(\Omega)}$は，それぞれSobolev空間$H^{m}(\Omega)$のセミノルム，ノルムを表す． また記号$(\cdot,\cdot)$は$L^{2}(\Omega)$または$(L^{2}(\Omega))^{2}$の内積を表す． Sobolev空間$W^{1,\infty}(\Omega)$は，1階までの弱導関数が$\Omega$上で本質的に有界であるような関数空間である． ベクトル値関数の空間$H(\mbox{div};\Omega)$を以下で定義する．

[TABLE]

次に以下のモデル問題を考える．

[TABLE]

ここで$\Gamma_{N},\Gamma_{D}$は，それぞれ互いに素な$\Omega$の境界$\partial\Omega$の部分集合であり，$\Gamma_{N}\cup\Gamma_{D}=\partial\Omega$を満たす． また$\frac{\partial}{\partial\mathbf{n}}$は$\partial\Omega$上の外向きの単位法線方向微分を表す．上の問題に対する弱形式は以下の式で与えられる．

[TABLE]

このとき，試行関数の関数空間$V$と試験関数の関数空間$V_{0}$は，以下で定義される．

[TABLE]

$\Gamma_{D}$が空集合となる場合，$V_{0}$の定義は以下のように変わる．

[TABLE]

また，

[TABLE]

関数空間$V_{0}$に関して以下のPoincaré不等式が成り立つ．

[TABLE]

Remark 1**.**

式(5)を満たす最適な定数$C_{p}$は固有値問題，

[TABLE]

の最小固有値$\lambda_{1}$に対して$C_{p}=1/\sqrt{\lambda_{1}}$で与えられ，特に$\Omega$が単位正方形$(0,1)^{2}$かつ$\Gamma_{D}=\partial\Omega$の場合は$C_{p}=1/(\sqrt{2}\pi)$となる． 一般的な領域と境界条件の設定の場合， $\lambda_{1}$の厳密な上下界は有限要素法によって得られる([12, 8, 9]など)． 文献[9]によると，Crouzeix-Raviart有限要素空間で定義される離散化固有値問題,

[TABLE]

の最小固有値$\lambda_{h,1}$を用いて，次の$\lambda_{1}$の評価式が成り立つ．

[TABLE]

2.2 有限要素法による近似解の構成

以降では議論を簡明にするためにモデル問題の境界条件に現れる$g_{D},g_{N}$を,それぞれ境界上で区分的な1次関数，定数関数とする． 一般の$g_{D},g_{N}$の場合には，それぞれの関数近似の誤差評価が必要となる．

領域$\Omega$の正則な三角形分割によって得られる三角形要素の族を$\mathcal{T}_{h}$とする． 三角形要素$K\in\mathcal{T}_{h}$に対して$K$の最長辺の長さを$h_{K}$で表し，メッシュサイズ$h$を以下で定義する．

[TABLE]

ここで関数空間$V,V_{0}$に対応する有限要素空間$V^{h},V_{0}^{h}(\subset H^{1}(\Omega))$を以下で定義する．

[TABLE]

このとき(2)の適合有限要素法による定式化は以下で与えられる．

[TABLE]

さらに，以下の有限要素空間を準備する．

• (a)

区分定数関数空間$X^{h}$:

[TABLE]

$\Gamma_{D}$が空集合となる場合，$X^{h}$を以下のように定義する．

[TABLE]

• (b)

[math]次Raviart-Thomas有限要素空間$W^{h}$:

[TABLE]

ここで$a_{K},b_{K},c_{K}$は各三角形要素$K$上の定数である．

上で述べた有限要素空間によって，問題(1)の混合有限要素法による定式化は以下で与えられる．

[TABLE]

ここで，後述の誤差解析に使用される射影作用素$\pi_{h}:L^{2}(\Omega)\to X^{h}$を定義する． 任意の$f\in L^{2}(\Omega)$に対して，$\pi_{h}f\in X^{h}$は

[TABLE]

を満たし，さらに以下の誤差評価が成り立つ．

[TABLE]

ここで三角形分割に依存する定数$C_{0}$は$C_{0}:=\max_{K\in\mathcal{T}_{h}}C_{0}(K)/h$で定義され，明らかな上界を持つ． 先行研究[3, 4, 6]において，最適な$C_{0}(K)$の値は，Bessel関数$J_{1}$の正の最小の根$j_{1,1}\approx 3.83171$によって$C_{0}(K):=h_{K}/j_{1,1}$と与えられることが報告されている．

3 有限要素解の大域的な事前誤差評価

この節では，以下の斉次混合境界値問題に対する有限要素解の大域的な事前誤差評価[8](または[12] の8章，[7]) を紹介する．

[TABLE]

上の問題に対する弱形式と適合有限要素法による定式化は，以下で与えられる．

[TABLE]

ここで，Galerkin射影作用素$P_{h}:V_{0}\to V_{0}^{h}$を紹介する． 任意の$v\in V_{0}$に対して，$P_{h}v$は以下の式を満たす．

[TABLE]

よって，$P_{h}\phi=\phi_{h}$となる．また，安定性に関する以下の評価式が成り立つ．

[TABLE]

ここで，$C_{p}$は式(5)に現れるPoincaré定数である． また，問題(10)の混合有限要素法による定式化は以下で与えられる．

[TABLE]

有限要素解の大域的な事前誤差評価に重要な役割を果たす定数$\kappa_{h}$を以下で定義する．

[TABLE]

これらの定数$C_{0},\kappa_{h}$を利用して，有限要素解の大域的な事前誤差評価を得ることができる．

Theorem 1** (Prager-Syngeの定理[13]).**

関数$\phi$を式(11)の厳密解とし，ベクトル値関数$\tilde{\mathbf{p}}\in H(\mbox{div};\Omega)$は以下を満たす．

[TABLE]

このとき，任意の$v\in V_{0}$に対して，Prager-Syngeの定理[13]に現れる以下のHypercircle式が成り立つ．

[TABLE]

Lemma 2** ([12]の定理8.2, [8]の定理3.2).**

与えられた$f\in L^{2}(\Omega)$に対して，関数$\overline{\phi}\in V_{0},~{}\overline{\phi}_{h}\in V_{0}^{h}$をそれぞれ以下の変分問題の解とする．

[TABLE]

このとき，$\overline{\phi}_{h}$は$\overline{\phi}$の近似解であり，以下の誤差評価が得られる.

[TABLE]

Remark 3**.**

$f$の近似$f_{h}$に対して$\mbox{div }\tilde{\mathbf{p}}_{h}+f_{h}=0$を厳密に満たす$\tilde{\mathbf{p}}_{h}$を構成することが，提案手法の特徴の一つである． 文献[15]は与えられる$f$に対して，要素の境界で$f$に対応する$\tilde{\mathbf{p}}$を構成している．

Theorem 2** (有限要素解の大域的な事前誤差評価[8]).**

与えられた$f\in L^{2}(\Omega)$に対して，$u$と$u_{h}$をそれぞれ変分問題(11)と(12)の解とする． このとき，以下の誤差評価が得られる．

[TABLE]

ここで，$C(h):=\sqrt{\kappa_{h}^{2}+(C_{0}h)^{2}}$である．

Remark 4**.**

方程式の解$u$が$H^{2}(\Omega)$に属するとき、 誤差評価式(20)と(21)の中の$C(h)$の代わりに、 Lagrange補間関数の誤差定数を使って値を与えることができる． 例えば、計算例§5.2では、一様直角三角メッシュに対して、$C(h)$の代わりに$C_{h}=0.493h$を使うことも可能[4, 8]．

Remark 5**.**

定理2において， 以下の事後誤差評価も得られる．

[TABLE]

Remark 6**.**

定数$\kappa_{h}$の計算法:  与えられた$f_{h}\in X^{h}$に対して，変分問題(18)の解$\phi_{h}$を対応させる線形作用素$R_{h}:X^{h}\to V^{h}$と， 変分問題(15)の解$\tilde{\mathbf{p}}_{h}$を対応させる線形作用素$T_{h}:X^{h}\to W^{h}$を定義し，さらに$Q_{h}:=\nabla\circ R_{h}-T_{h}$とする． このとき$\kappa_{h}$は以下のように計算されることが分かる．

[TABLE]

ここで$X^{h},V^{h},W^{h}$は有限次元ベクトル空間なので，上の式は行列の固有値問題に帰着される． より詳細な説明については先行研究[8, 7]を参照する．

Remark 7**.**

式(1)に示した一般の混合境界値問題の場合，与えられるデータ$f,g_{D},g_{N}$から解$u$または有限要素解$u_{h}$への対応はアフィン写像になる．

4 重み付きHypercircle式と有限要素解の局所事後誤差評価

この節では，Hypercircle式を用いた有限要素解の局所事後誤差評価を提案する． 関心のある領域$S$は$\Omega$の部分領域として与えられ，$S$に対応する重み関数(またはCutoff関数)$\alpha$が導入される． 次に，重み関数$\alpha$を用いて重み付きHypercircle式を導入し，前節の技巧を応用して有限要素解の局所事後誤差評価を与える．

4.1 重み関数

関心のある領域$S$に対して，$S$を囲う帯状領域を$B_{S}:=\{x\in\Omega\setminus\overline{S};~{}\mbox{dist}(x,\partial S)<\varepsilon\}$とし(Fig.2-(a), Fig.2-(b)参照)，さらに部分領域$\Omega^{\prime}=S\cup B_{S}$を定義する． 次に，重み関数$\alpha\in W^{1,\infty}(\Omega)$を定義する． 重み関数$\alpha$は$\overline{\Omega^{\prime}}$を台(Support)とする区分的な多項式関数であり，以下の性質を持っている．

[TABLE]

特に，$S$が長方形の場合(Fig.2-(a)参照)，$S$は2つの1次元開区間の直積によって表されるため， 重み関数$\alpha$はそれぞれの開区間に対応する重み関数$\alpha_{1}(x)$と$\alpha_{2}(y)$の積によって構成できる． 開区間$(a,b)$に対応する重み関数$\alpha_{(a,b)}$は以下のように構成できる．また，$\alpha_{(a,b)}$のグラフをFig.2-(c)に示す．

[TABLE]

ただし，$\varepsilon$は帯状領域$B_{S}$の幅である．

次に，重み関数$\alpha$によって，重み付き内積$(\cdot,\cdot)_{\alpha}$と重み付きノルム$\|\cdot\|_{\alpha}$を導入する．

• (a)

重み付き内積$(\cdot,\cdot)_{\alpha}$:  $f,g\in L^{2}(\Omega)$または$f,g\in(L^{2}(\Omega))^{2}$に対して，

[TABLE]

• (b)

重み付きノルム$\|\cdot\|_{\alpha}:$  $f\in L^{2}(\Omega)$に対して，

[TABLE]

このとき，以下の不等式が成り立つことは容易に分かる．

[TABLE]

4.2 重み付きHypercircle式

この節では定理16の拡張として，重み付きノルムに関するHypercircle式を与える．

Theorem 3**.**

関数$u$を式(2)の厳密解とし，ベクトル値関数$\mathbf{p}\in H(\mbox{div};\Omega)$は以下を満たす．

[TABLE]

このとき，任意の$v\in V$に対して，以下の重み付きHypercircle式が成り立つ．

[TABLE]

Proof.

$\|\nabla v-\mathbf{p}\|_{\alpha}^{2}=\|(\nabla v-\nabla u)+(\nabla u-\mathbf{p})\|_{\alpha}^{2}$の展開式を考える．

[TABLE]

上の展開式の交差項$(\nabla v-\nabla u,\nabla u-\mathbf{p})_{\alpha}$について，$v-u$を$w$とおき，交差項を二つの項$(\nabla w,\nabla u)_{\alpha}$, $(\nabla w,\mathbf{p})_{\alpha}$に分けて，それぞれの評価を与える．

まず，$(\nabla w,\nabla u)_{\alpha}$について，$\alpha w$に対して合成関数の微分公式$\alpha\nabla w=-w\nabla\alpha+\nabla(\alpha w)$から以下の式が成り立つことが分かる．

[TABLE]

変分式(2)のテスト関数を$\alpha w$にすると，$(\nabla w,\nabla u)_{\alpha}$について以下の式が得られる．

[TABLE]

次に，$(\nabla w,\mathbf{p})_{\alpha}$について，Greenの公式を用いて以下の式が得られる．

[TABLE]

さらに，定理の仮定$\mbox{div }\mathbf{p}+f=0$を利用することで以下の評価が得られる．

[TABLE]

よって，式$(\ref{eq:alpha-uw}),~{}(\ref{eq:alpha-pw})$の差を取ることで，交差項に関して以下の等式が成り立つ．

[TABLE]

式 (35)の領域$\Omega$を$\Omega^{\prime}$に入れ替えて，さらにHölderの不等式を用いることで，次の評価式が得られる．

[TABLE]

最後に，式(32)に対して上の結果を用いることで結論が得られる． ∎

Remark 8**.**

定理31の式(33), (34)において，関心のある領域$S$と$\Omega$が境界の一部を共有する場合でも，$\Gamma_{N}$上の積分は交差項の評価に影響を与えないことが分かる．

4.3 有限要素解の局所事後誤差評価

有限要素解の局所事後誤差評価を検討するために，与えられた$f\in L^{2}(\Omega)$に対して以下の補助問題を準備する．

[TABLE]

Lemma 9**.**

$|u-\overline{u}|_{H^{1}(\Omega)}$に関して以下の評価が成り立つ．

[TABLE]

Proof.

$\|\nabla(u-\overline{u})\|_{\Omega}$について，$u$と$\overline{u}$の定義から以下の等式が成り立つ．

[TABLE]

射影作用素$\pi_{h}$に関する誤差評価式(9)から以下の評価式を得る．

[TABLE]

試験関数$v$を$u-\overline{u}$とすることで，以下の評価が得られる．

[TABLE]

∎

Lemma 10**.**

与えられた$f\in L^{2}(\Omega)$に対して，変分問題(2)と(6)の解)の解をそれぞれ$u$と$u_{h}$とし， 変分問題(36)と(37)の解をそれぞれ$\overline{u},\overline{u}_{h}$とする. さらに，$\mathbf{p}_{h}\in H(\mbox{div};\Omega)$は以下の条件を満たすベクトル値関数とする．

[TABLE]

このとき$\|\nabla(\overline{u}-u_{h})\|_{\alpha}$に関して，以下の評価式が得られる．

[TABLE]

Proof.

定理31において$f=\pi_{h}f$とすると，以下の重み付きHypercircle式が得られる．

[TABLE]

式(30)の結果により，以下の評価式が得られる．

[TABLE]

∎

Theorem 4** (有限要素解の局所事後誤差評価).**

与えられた$f\in L^{2}(\Omega)$に対して，$u$と$u_{h}$を，それぞれ変分問題(2)と(6)の解とし， $\mathbf{p}_{h}\in H(\mbox{div};\Omega)$は以下の条件を満たすベクトル値関数とする．

[TABLE]

このとき，以下の誤差評価が得られる．

[TABLE]

ただし，

[TABLE]

ここで，$C_{p}$は式(5)に現れるPoincaré定数である．

Proof.

まず，式(30)を用いて，$\|\nabla(u-u_{h})\|_{S}$に関する以下の展開式が得られる．

[TABLE]

式(41), (38)から，以下の評価が得られる．

[TABLE]

定理31において$f=\pi_{h}f,v=u_{h}$とすると，重み付きHypercircle式(31)によって，以下の評価式が得られる．

[TABLE]

式(42),(43)から，$\|\nabla(u-u_{h})\|_{S}$に関して以下の評価式を得る．

[TABLE]

次に，$\|\nabla(u-u_{h})\|_{S}$の評価式(44)にある$\|\overline{u}-u_{h}\|_{\Omega},~{}\|\nabla\overline{u}-\mathbf{p}_{h}\|_{\Omega}$を検討する．

• (a)

$\|\nabla\overline{u}-\mathbf{p}_{h}\|_{\Omega}$の評価:  $\overline{u}\in V$と$\mbox{div }\mathbf{p}_{h}+\pi_{h}f=0,~{}\mathbf{p}_{h}\cdot n=g_{N}\mbox{ on }\Gamma_{N}$を満たすベクトル値関数$\mathbf{p_{h}}\in W^{h}$に対して，以下のHypercircle式が得られる．

[TABLE]

よって，以下の評価式が得られる．

[TABLE]

• (b)

$\|\overline{u}-u_{h}\|_{\Omega}$の評価:  以下の双対問題を定義する．

[TABLE]

双対問題の解$\phi$に対して，$\phi_{h}:=P_{h}\phi$とすると（$P_{h}$は式(13)で定義したGalerkin射影である），$\|\overline{u}-u_{h}\|_{\Omega}$に関して以下の式が得られる． ここで$u_{h}$は$\overline{u}$のGalerkin射影ではないことに注意する．

[TABLE]

さらに式(45)を用いて，以下の展開式が得られる．

[TABLE]

以降では展開式の$(\nabla\phi_{h},\nabla(\overline{u}-u_{h}))$と$\|\nabla(\phi-\phi_{h})\|_{\Omega}$について，それぞれの評価を与える． まず，展開式の$(\nabla\phi_{h},\nabla(\overline{u}-u_{h}))$について，射影作用素$\pi_{h}$に関する誤差評価式(9)から以下の評価式を得る．

[TABLE]

上の不等式の最後に，Galerkin射影作用素に関する安定性の評価式(14)を用いた． 次に，展開式の$\|\nabla(\phi-\phi_{h})\|_{\Omega}$ついて，定理2の事前誤差評価式(20)から以下の評価式を得る．

[TABLE]

よって，$\|\overline{u}-u_{h}\|_{\Omega}$に関して以下の評価式が得られる．

[TABLE]

目標としている式(44)の根号の中の評価は式(43)に対して，評価式(46), (47)を用いることで得られる．

[TABLE]

以上により定理の結論を得る． ∎

Remark 11**.**

定理4の中に現れる$Err^{(1)},Err^{(2)}$および$Err^{(3)}$について，$Err^{(1)},Err^{(2)}$は有限要素解の大域的な誤差に関わる項で，$Err^{(3)}$は有限要素解の局所的な誤差に関わる主要項である． $Err^{(1)},Err^{(2)}$のオーダーは$O(h^{1.5})$，$Err^{(3)}$のオーダーは$O(h^{1})$となることが期待される．特に，メッシュサイズ$h$が十分小さいとき主要項$Err^{(3)}$が支配的になることが期待される．

5 数値実験

5.1 準備

5.2, 5.3節では，2次元領域で定義された，以下の斉次Dirichlet問題に対して数値実験を行う．

[TABLE]

5.4節では，抵抗率測定の問題に応じて，3次元領域で定義された混合境界値問題について検討する．

5.2, 5.3節では，最初に帯状領域$B_{S}$の幅を決定する．$B_{S}$の幅を大きく取ると，誤差評価式(40)に現れる主要項$Err^{(3)}$が大きなる一方で，$B_{S}$の幅を小さく取ると，領域全体の誤差評価に関わる項${Err^{(1)}},{Err^{(2)}}$に含まれる$\|\nabla\alpha\|_{L^{\infty}(\Omega)}$が大きくなる． そこで，$B_{S}$の幅を$0.2\sim 0.375$まで$0.05$ずつ取って，$B_{S}$の幅の影響を調べる．次に，適切な$B_{S}$の幅を取って，メッシュサイズ$h$と誤差評価式(40)の関係を調べる．

5節を通して，誤差評価式(40)に現れる記号${Err^{(1)}},{Err^{(2)}},{Err^{(3)}}$を引き続き使用し，誤差評価式(40)を計算することで得られた値を$\overline{E}_{L}$で表す． すなわち，

[TABLE]

さらに以下の記号を導入する．

[TABLE]

5.2 正方形領域

正方形領域の場合，問題(50)の厳密解$u=\sin(\pi x)\sin(\pi y)$を利用して$\overline{E}_{L}$と$E_{L}$の比較を行うことができる． 数値実験の条件を以下に示す．

1. (a)

領域$\Omega:=(0,1)^{2}$， 部分領域$S:=(0.375,0.625)^{2}$． 2. (b)

一様メッシュのサイズ$\displaystyle h:=\frac{1}{8},~{}\frac{1}{16},~{}\frac{1}{32},~{}\frac{1}{64},~{}\frac{1}{128}$．

帯状領域$B_{S}$の幅の選択 　メッシュサイズ$h=1/64,1/128$の場合に帯状領域$B_{S}$の幅と$\overline{E}_{L}$の関係をFig.3-(a),(b)に示す． この図から最適な$B_{S}$の幅は約$0.3$であり，この付近での幅の$1\%$の変化に対して，$\overline{E}_{L}$の変化が約$1\%$となっていることが分かる． したがって，最適な$B_{S}$の幅の周りでは幅の変化が与える誤差評価への影響は小さいと考えられる．

メッシュサイズ$h$と$\overline{E}_{L}$の関係をTable 2, Fig.4-(a),(b)に示す．ここで，$B_{S}$の幅を$0.3$とした．

Table 2, Fig.4-(a),(b)から，$\overline{E}_{L}$は$E_{L}$の上界を与えていることがわかり， $h\leq 1/32$のときに誤差評価式(40)の主要項$Err^{(3)}$が支配的である事も分かる．したがって，注意11で述べたことが数値実験よって確認できた．

5.3 L字型領域

L字型領域上の問題(50)の解$u$は，L字型領域の非凸な角点での特異性によって$u\notin H^{2}(\Omega)$となることが知られている． 提案する誤差評価式(40)は，解の特異性にも自然に対応できるほか，$S$と$\Omega$が境界の一部を共有する場合でも誤差評価が可能である． この例では，以下の条件で数値実験を行う．

1. (a)

領域$\Omega$: $(0,1)^{2}\setminus[0.5,1]^{2}$，部分領域$S:=(0.375,0.625)^{2}\setminus[0.5,1]^{2}$． 2. (b)

一様メッシュのサイズ$\displaystyle h:=\frac{1}{8},~{}\frac{1}{16},~{}\frac{1}{32},~{}\frac{1}{64},~{}\frac{1}{128}$．

帯状領域$B_{S}$の幅の選択 　メッシュサイズ$h=1/64,1/128$の場合に帯状領域$B_{S}$の幅と$\overline{E}_{L}$の関係をFig.5-(a),(b)に示す． この例の領域の設定によって，$B_{S}$の幅は最大$0.375$まで取れるので，以降の誤差評価の検討では$B_{S}$の幅を$0.375$とする．

メッシュサイズ$h$と$\overline{E}_{L}$の関係をTabel 2,  Fig.6-(a),(b)に示す．

5.4 抵抗率測定法への応用

提案する誤差評価手法の応用例として，前述の四探針法への応用を紹介する．

まず，抵抗率測定法を支配する数理モデルを準備する． 測定試料の占める$\mathbb{R}^{3}$の領域を$\Omega$で表し，抵抗率測定に用いられる4本の探針$A$, $B$, $C$, $D$の位置を，それぞれ$P_{A}$, $P_{B}$, $P_{C}$, $P_{D}$で表す． さらに，探針$A$および$D$の接触部分を，それぞれ$\Gamma_{A}$, $\Gamma_{D}$ ($\subset\partial\Omega$)と表す． このとき，試料内部の電位$u$を支配する数理モデルの一つとして，以下の偏微分方程式の境界値問題が考えられる．

[TABLE]

ここで，$|\Gamma_{A}|,|\Gamma_{D}|$は，それぞれ探針$A,D$の接触面積を表す． このとき、式(4)で定義される関数空間$V_{0}$を使用する．抵抗率測定法の問題(51)で関心が持たれるのは，探針$B$および$C$における電位である． また、有限要素法の計算では、２次元要素ではなく、四面体の分割における３次元有限要素を使用する．

数値実験のために$\Omega,~{}P_{A},~{}P_{B},~{}P_{C},~{}P_{D},~{}\Gamma_{A},~{}\Gamma_{D}$および部分領域$S$について，以下で設定する(Fig.7-(a)参照).

1. (a)

領域$\Omega$:  $(0,8)\times(0,5)\times(0,0.2)$  (直方体領域) 2. (b)

探針位置:

[TABLE] 3. (c)

接触部分$\Gamma_{A},~{}\Gamma_{D}$:  $\Gamma_{A}=B(P_{A};0.02),~{}\Gamma_{D}=B(P_{D};0.02)$ 4. (d)

関心のある部分領域$S$:  $(3.5,4.5)\times(2.25,2.75)\times(0.0,0.2)$

ただし，$B(P;r)$は試料の表面における中心$P$，半径$r$の円である．

次に，問題(51)の局所事後誤差評価について検討する． 問題(51)では，境界条件が区分的な定数関数によって与えられるので，定理4の仮定を満足する． さらに，$f=0$であるため当該誤差評価式(40)は，以下に示す簡便な式となる．

[TABLE]

この例では帯状領域$B_{S}$の幅を0.2とし，有限要素空間として2次の適合有限要素空間および1次Raviart-Thomas有限要素空間が用いられる. さらに$X^{h}$を区分的な1次多項式の関数空間とする． この場合，射影作用素$\pi_{h}$の誤差評価式9については，より良い定数が利用できるが，この例では引き続き$C_{0}$を利用する． 局所事後誤差評価の条件Table.3に，計算に用いたメッシュをFig.7-(b)に，計算結果をTable.4に示す．

ただし$A_{M},~{}A_{C}$は，それぞれ混合有限要素法，適合有限要素法の全体剛性行列である．

6 結論

本論文では有限要素解の定量的な局所事後誤差評価について検討した． 正方形領域とL字型領域の上で定義されたPoisson方程式の境界値問題に対して数値実験を行い，提案手法の有効性について確認できた． さらに，抵抗率測定法に関する数値例を示した．抵抗率測定法に使用する$|u(P_{B})-u(P_{C})|$の誤差評価については，今後の研究で検討する予定である．