数理解析研究所講究録
第 1791 巻 2012 年 206-215

206

A residual bound evaluation of operator equations with
Raviart-Thomas finite element
高安 亮紀 (Akitoshi Takayasu) 1

早稲田大学基幹理工学研究科

Graduate School of Fundamental Science and Engineering, Waseda University
早稲田大学 劉 雪峰 (Xuefeng Liu) 2
Faculty of Science and Engineering, Waseda University,
CREST, JST
早稲田大学 大石 進一 (Shin’ichi Oishi) 3
Department of Applied Mathematics,

Faculty of Science and Engineering, Waseda University,

CREST, JST
Abstract –In this article, a residual evaluation of operator equation is considered in the frameworkof computer-assisted proof. Our computer-assisted approach ensures the existence and local uniqueness of
weak solutions to some nonlinear partial differential equations. Based on Newton-Kantorovich theorem, ournumerical method is a variant of existing methods such as [1, 2, 3, 4]. Residual evaluation for operator
equation plays important role in validating numerical solutions. In order to get accurate residual evaluation,some smoothing techniques have been proposed. Main objective of this article is to obtain a sharp boundevaluation with high order Raviart-Thomas mixed finite element.

Introduction

1

with arbitrary shape. is the set of real numbers. In this article,
be bounded polygonal domain in
we are concerned with Dirichlet boundary value problem of the semi-linear elliptic equation of the form:Let$\mathbb{R}^{2}$$\Omega$$\mathbb{R}$(1)$\{\begin{array}{ll}-\Delta u=f(\nabla u, u, x), in \Omega,u=0, on \partial\Omega\end{array}$is assumed to be Fr\’echet differentiable. For example, $f(\nabla u, u, x)=-b\cdot\nabla u-$satisfies this condition. Verified
andwith $b(x)\in(L^{\infty}(\Omega))^{2},$computation approach will be adopted to explore the existence and local uniqueness of weak solution of(1). Namely, if an approximate solution is given by certain numerical method, we will try to validate theexistence of exact solution in the neighbourhood of the approximation. In the classical analysis of variational
theory, weak solution of Dirichlet boundary problem (1) is defined in variational form:
where

$f$

:

$H_{0}^{1}(\Omega)arrow L^{2}(\Omega)$$cu+c_{2}u^{2}+c_{3}u^{3}+g$

Find$c,$ $c_{2},$ $c_{3}\in L^{\infty}(\Omega)$$u\in H_{0}^{1}(\Omega)$$g\in L^{2}(\Omega)$, satisfying $(\nabla u, \nabla v)=(f(\nabla u, u, x), v)$ ,for all$v\in H_{0}^{1}(\Omega)$

.

(2)

Here,
$( \nabla u, \nabla v):=\int_{\Omega}\nabla u\cdot\nabla vdx$and$(f( \nabla u, u,x), v):=\int_{\Omega}f(\nabla u, u, x)vdx$.

Now we put $V=H_{0}^{1}(\Omega)$ and rewrite $f(\nabla u, u, x)$ as $f(u)$ for simple form. Let us define linear and nonlinear$\mathcal{N}:Varrow V$ ,$:=(\nabla u, \nabla v),$$:=(f(u), v)$ . Furthermore, we define $\mathcal{F}:Varrow V$
Au,
operators. The original problem (1) is equivalent to the following nonlinear operator equation:as$\mathcal{A},$$\mathcal{F}(u)$$($$v)_{V}$$(\mathcal{N}(u), v)_{V}$$:=\mathcal{A}u-\mathcal{N}(u)$Find $u\in V$, satisfying
1 takitoshi Qsuou. waseda. jp

waseda. jp
ishi Owaseda. jp$2_{xf1i}$ uQaoni.
$3_{Q}$$\mathcal{F}(u)=0$

.

(3)

207
is assumed to be Fr\’echet differentiable mapping. Let
be an approximate solution toeq.(3). Fr\’echet derivative of at is denoted by: $Varrow V$ . In order to verify the existence and local
uniqueness of the exact solution in the neighborhood of , we consider to apply the Newton-Kantorovich
theorem [5, 6] to eq.(3).
$\mathcal{F}:Varrow V$$\hat{u}\in V_{h}\subset V$$\mathcal{F}$$\hat{u}$$\mathcal{F}’[\hat{u}]$$\hat{u}$Theorem 1. Assuming Frechet derivative$\mathcal{F}’[\hat{u}]$is nonsingular and satisfies$\Vert \mathcal{F}’[\hat{u}]^{-1}\mathcal{F}(\hat{u})\Vert_{V}\leq\alpha$,

be a closed ball centered at withfor a certain positive . Then, let. Let also
radiusbe an open ball in V. We assume that for a certain positive , it holds:$\overline{B}(\hat{u}, 2\alpha):=\{v\in V:\Vert v-\hat{u}\Vert_{V}\leq 2\alpha\}$$\alpha$$2\alpha$$\omega$$\Vert \mathcal{F}’[\hat{u}]^{-1}(\mathcal{F}’[v|-\mathcal{F}’[w|)\Vert_{V,V}\leq\omega\Vert v-w\Vert_{V},$If$\alpha\omega\leq\frac{1}{2}$$\hat{u}$$D\supset\overline{B}(\hat{u}, 2\alpha)$$\forall v,$$w\in D$.

holds, then there is a solution $u\in V$ of $eq.(3)$ satisfying$\Vert u-\hat{u}\Vert_{V}\leq\rho:=\frac{1-\sqrt{1-2\alpha\omega}}{\omega}$

Furthermore, the solution$u$is unique in$\overline{B}(\hat{u}, \rho)$

.

(4)

.

Remark 1. To apply Newon-Kantorovich theorem, we will calculate the constants below explicitly.
$\Vert \mathcal{F}’[\hat{u}]^{-1}\Vert_{V,V}\leq C_{1}$$\Vert \mathcal{F}(\hat{u})\Vert_{V}\leq C_{2,h}$,$\Vert \mathcal{F}’[v]-\mathcal{F}’[w]\Vert_{V,V}\leq C_{3}\Vert v-w\Vert_{V}$

,

(5)

,

(6)
$\forall v,$.

(7)$w\in D\subset V$

Therefore, if
is confirmed by verified computations, then the existence and local uniquenessof the solution are proved numerically based on Newton-Kantoromch theorem.
$C_{1}^{2}C_{2},{}_{h}C_{3}\leq 1/2$

The main topic of this article is to evaluate the residual bound for$\Vert \mathcal{F}(\hat{u})\Vert_{V}\leq C_{2,h}$$\mathcal{F}(\hat{u}),$

$i.e$

.

.

(8)

In the following, we would like to introduce several ways to evaluate eq.(S). Suppose function
to be
an approximation of exact solution of eq.(3), where. Our aim isis certain finite element subspace
to obtain good estimation of this residual bound. First, we introduce several evaluation methods in Section
2. Second, we show numerical results in Section 3 to demonstrate the efficiency of our proposed method.
For reader convenience, we write down the details for implementation of Raviart-Thomas element method
in appendix.
$\hat{u}\in V_{h}$$V_{h}\subset V$$V_{h}$$s$

Several ways for residual evaluation

2

In this section, we would like to consider the residual evaluation in the form of$\Vert \mathcal{F}(\hat{u})\Vert_{V}=\sup_{0\neq v\in V}\frac{(\mathcal{A}\hat{u}-\mathcal{N}(\hat{u}),v)_{V}}{\Vert v||_{V}}=\sup_{0\neq v\in V}\frac{|(\nabla\hat{u},\nabla v)-(f(\hat{u}),v)|}{||v\Vert_{V}}$in several ways. If an approximate solution satisfies$\hat{u}\in H^{2}(\Omega)\cap V_{h}$, it follows$\Vert \mathcal{F}(\hat{u})\Vert_{V}=\sup_{0\neq v\in V}\frac{|(\nabla\hat{u},\nabla v)-(f(\hat{u}),v)|}{||v\Vert_{V}}=\sup_{0\neq v\in V}\frac{|(-\Delta\hat{u},v)-(f(\hat{u}),v)|}{||v\Vert_{V}}\leq C_{e,2}\Vert\triangle\hat{u}+f(\hat{u})\Vert_{L^{2}}$.

(9)$(2\leq p<\infty)$ for $u\in V$ .
Here,means Sobolev embedding constant, which satisfies
We point out that the evaluation (9) does not work whenfinite element functions, suchis taken as(piecewise linear) oras(piecewise quadratic) elements. This is because
does not belong to$C_{e,p}$$\Vert u\Vert_{L^{p}}\leq C_{e,p}|u|_{H^{1}},$$s$$V_{h}$$P_{1}$$C^{0}$$L^{2}(\Omega)$$\triangle\hat{u}$$P_{2}$

anymore.

To weaken the condition on , we will introduce several methods that do not need the-regularity of
approximate solution. The first method to be introduced is fast but gives little rough bound. The second
one has accurate estimation with smoothing technique. The third one is based on Raviart-Thomas mixed
finite elements [9, 10, 11], which can provide better bound for residue if higher order elements are used.
$\hat{u}$$H^{2}$

208
2.1

Simple bounds

span. Let $u_{h};=\prime P_{h}u\in V_{h}$ be an
be a finite element subspace of , such thatLet$(\nabla(u-u_{h}),\nablav_{h})=0,$
$u\in
V$
In
,
this part, we will show simple
orthogonal projection ofdefined as. From the classicalupper bound of residue. In the following, we denoteby the projection of $i.e$ .
error analysis, such as Aubin-Nitsche’s trick, we have$V$

$V_{h}$$V_{\hslash};=$$\{\phi_{1}, \ldots, \phi_{n}\}$$\forall v_{h}\in V_{h}$$\Vert v-v_{h}\Vert_{L^{2}}\leq C_{M}\Vert v-v_{h}\Vert_{V}$$\Vert v-v_{h}\Vert_{V}\leq\Vert v\Vert_{V}$$P_{h}v$$v,$$v_{\hslash}$and$\Vert v_{h}$

llv

,

(10)

$\leq\Vert v\Vert_{V}$

.

(11)

Here
is a priori error constant for projection . The full discussion of this constant on arbitrary domain, the residual bound of eq.(8) is given using inequalities (10) and (11)is shown in [12]. For$C_{M}$$\mathcal{P}_{h}$$v_{h}\in V_{h}$$\Vert \mathcal{F}(\hat{u})\Vert_{V}$$=$$\sup_{0\neq v\in V}\frac{|(\nabla\hat{u},\nabla v)-(f(\hat{u}),v)|}{\Vert v\Vert_{V}}$$=$$\sup_{0\neq v\in V}\frac{|(\nabla\hat{u},\nabla(v-v_{h}))-(f(\hat{u}),v-v_{h})+(\nabla\hat{u},\nabla v_{h})-(f(\hat{u}),v_{h})|}{\Vert v||_{V}}$$\leq$$\sup_{0\neq v\in V}\frac{|(f(\hat{u}),v-v_{h})|}{||v\Vert_{V}}+\sup_{0\neq v\in V}\frac{|(\nabla\hat{u},\nabla v_{h})-(f(\hat{u}),v_{h})|}{\Vert v\Vert_{V}}$$C_{M}$$\leq$

where the quantity$C_{r}$

llf

is defined by the following procedure$\sup_{0\neq v\in V}\frac{|(\nabla\hat{u},\nabla v_{h})-(f(\hat{u}),v_{h})|}{\Vert v\Vert_{V}}$$=$$\leq$

Let beLet us put
given by$\epsilon_{i}$$\epsilon_{i}$$0^{0\neq v\in V} \sup_{=v_{h}\in V_{h}}\frac{|(\nabla\hat{u},\nabla v_{h})-(f(\hat{u}),v_{\hslash})|}{||v\Vert_{V}}+0^{0\neq v\in V}\sup_{\neq v_{h}\in V_{h}}\frac{|(\nabla\hat{u},\nabla v_{h})-(f(\hat{u}),v_{h})|}{\Vert v_{h}||_{V}}\cdot\frac{\Vert v_{h}\Vert_{V}}{||v\Vert_{V}}$$\sup_{0\neq v_{h}\in V_{h}}\frac{|(\nabla\hat{u},\nabla v_{h})-(f(\hat{u}),v_{h})|}{\Vert v_{h}\Vert_{V}}=:C_{r}$$:=(\nabla\hat{u}, \nabla\phi_{i})-(f(\hat{u}), \phi_{i}),$$c:=(c_{1}, \ldots, c_{n})^{t}$$(\nabla\phi_{i}, \nabla\phi_{j})$(12)$(\hat{u})\Vert_{L^{2}}+C_{r}$and. Since. Let further$(i=1, \ldots, n)$$\epsilon;=(\epsilon_{1}, \ldots, \epsilon_{n})^{t}$. Then, C. follows$v_{h}\in V_{h}$

$D$

be

.

.
, we can expressas$n\cross n$ matrix whose $(i,j)$ -elements are$v_{h}$$v_{h}$$:= \sum_{i=1}^{n}c_{i}\phi_{i}$$C_{r}= \sup_{0\neq v_{h}\in V_{h}}\frac{|(\nabla\hat{u},\nabla v_{h})-(f(\hat{u}),v_{h})|}{||v_{h}\Vert_{V}}=\sup_{c\in R^{n}}\frac{|\sum_{i=1}^{n}q\epsilon_{i}|}{\sqrt{c^{t}Dc}}\leq\sup_{c\in R^{n}}\frac{|c|_{l^{2}}|\epsilon|_{l^{2}}}{\sqrt{c^{t}Dc}}\leq\Vert D^{-1}\Vert_{2}|\epsilon|_{l^{2}}$

.

(13)

From inequalities (12) and (13), we obtain$\Vert \mathcal{F}(\hat{u})\Vert_{V}\leq C_{M}\Vert f(\hat{u})\Vert_{L^{2}}+\Vert D^{-1}\Vert_{2}|\epsilon|_{l^{2}}$

2.2

.

(14)

Accurate bounds with a smoothing technique

The simple bound (14) is a rough bound. overestimation often causes failure in verification. Next, anothermethod for evaluating the residual bound is introduced. This is based on the smoothing technique proposedby smooth function.
by N. Yamamoto et. al. [13]. Here, smoothing means to approximate vector
According to [13], if (piecewise linear) elements are used for approximate solutions, the residual evaluation
becomes almost the same as the rough bound in (14). On the other hand, using higher order element, thissmoothing technique works very well [14]. Let
be a finite element subspace that does not vanish
on boundary of . Let $p_{h}\in(X_{h})^{2}$ be the vector function defined by$\nabla\hat{u}$$P_{1}$$X_{h}\subset H^{1}(\Omega)$$\Omega$$(p_{h}-\nabla\hat{u},v^{*})=0$,$\forall v^{*}\in(X_{h})^{2}$

Namely it is the -projection of
to $p_{h}\in(X_{h})^{2}$ .
$p_{h}[13]$:Further the following Green formula holds for$L^{2}$$\nabla\hat{u}\in(L^{2}(\Omega))^{2}$$p_{h}$.

(15)

makes the quantity$\Vert p_{h}-\nabla\hat{u}\Vert_{L^{2}}$small.

$s$$(p_{h}, \nabla v)+(divp_{h}, v)=0$ ,.

$\forall v\in V$

(16)

209
Therefore, using

and inequalities (10), (11), (13) and eq.(16), we have$p_{h}$$\Vert \mathcal{F}(\hat{u})\Vert_{V}$$=$$\sup_{0\neq v\in V}\frac{|(\nabla\hat{u},\nabla v)-(f(\hat{u}),v)|}{||v\Vert_{V}}$$=$$\sup_{0\neq v\in V}\frac{|(\nabla\hat{u},\nabla(v-v_{h}))-(f(\hat{u}),v-v_{h})+(\nabla\hat{u},\nabla v_{h})-(f(\hat{u}),v_{h})|}{||v||_{V}}$$\leq$$\sup_{0\neq v\in V}\frac{|(\nabla\hat{u},\nabla(v-v_{h}))-(f(\hat{u}),v-v_{h})|}{\Vert v||_{V}}+C_{r}$$\leq$$\sup_{0\neq v\in V}\frac{|(\nabla\hat{u}-p_{h},\nabla(v-v_{h}))+(p_{h},\nabla(v-v_{h}))-(f(\hat{u}),v-v_{h})|}{\Vert v\Vert_{V}}+C_{r}$$\leq$$\sup_{0\neq v\in V}\frac{\Vert\nabla\hat{u}-p_{h}\Vert_{L^{2}}\Vert v-v_{h}\Vert_{V}+\Vert divp_{h}+f(\hat{u})\Vert_{L^{2}}\Vert v-v_{h}\Vert_{L^{2}}}{\Vert v||_{V}}+C_{r}$$\Vert\nabla\hat{u}-p_{h}\Vert_{L^{2}}+C_{M}\Vert divp_{h}+f(\hat{u})\Vert_{L^{2}}+\Vert D^{-1}\Vert_{2}|\epsilon|_{l^{2}}$$\leq$.

(17)

One can use the bound (17) instead of (14). The smoothing element
is obtained by solving an additional
linear equation (15), which takes extra computational costs. Meanwhile, for a certain good approximate
solution, $e.g$ . using(piecewise quadratic) elements, residual bound (17) becomes drastically small [14].
$p_{h}$$P_{2}$

Remark 2. One can consider another evaluation with $H(div, \Omega)$ -smoothing elements $l4/$ . A smoothing
and $divq+f(\hat{u})\approx O$ yields

function $q\in H(div, \Omega)$ satisfying$q\approx\nabla\hat{u}$$\Vert \mathcal{F}(\hat{u})\Vert_{V}\leq\Vert\nabla\hat{u}-q\Vert_{L^{2}}+C_{e,2}\Vert divq+f(\hat{u})\Vert_{L^{2}}$.

One feature of thrs estimation is that it seeks the smoothing function in $q\in H(div, \Omega)\supset(H^{1}(\Omega))^{2}$ , which
can provide better approximation of , compared with the one in $eq.(15)$ .
$\nabla\hat{u}$2.3Raviart-Thomas mixed finite element on triangle element

Inspired by Remark 2, we are concerned with a smoothing technique using mixed finite elements as below.
Here, we would like to introduce Raviart-Thomas mixed finite element [9, 10, 11]. We follow discussions in[10, 11]. Let $H(div, \Omega)$ denote the space of vector functions such that
$H(div, \Omega):=\{\psi\in(L^{2}(\Omega))^{2} :div\psi\in L^{2}(\Omega)\}$Let$K_{h}$

be a triangle element in triangulation of$P_{k}(K_{h})$$R_{k}(\partial K_{h})$Functions of$R_{k}(\partial K_{h})$$RT_{k}(K_{h})$. We define: the space of polynomials of degree less than$:=\{\varphi\in L^{2}(\partial K_{h}) :$\{q\in(L^{2}(K_{h}))^{2}$:$k$, for any edge

\varphi|_{e}:\in P_{k}(e_{i})\}$

are polynomials of degree$;=$$\Omega$.

$\leq k$on each side$e_{i}$

on$e_{i}$

of$K_{h}$,$\partial K_{h}$.

of $K_{h}(i=1,2,3)$ . For $k\geq 0$ , we define$q=(\begin{array}{l}a_{k}b_{k}\end{array})+c_{k}\cdot(\begin{array}{l}xy\end{array}),$$a_{k},$ $b_{k},$$c_{k}\in P_{k}(K_{h})\}$.

The dimension of $RT_{k}(K_{h})$ is $(k+1)(k+3)$ . We now introduce basic result about $RT_{k}(K_{h})$ spaces.
Proposition 1. Letbe subtense of vertex $i(=1,2,3)$ andvector on boundary . For $q\in RT_{k}(K_{h})$ , it follows$e_{i}$$\vec{n}_{|e}:=(n_{1}^{(i)}, n_{2}^{(i)})^{t}$be an outward unit normal$e_{i}$$\{\begin{array}{l}divq\in P_{k}(K_{h}),q\cdot\vec{n}_{|e_{i}}\in R_{k}(\partial K_{h}).\end{array}$Moreover, the divergence operator is surjective from $RT_{k}(K_{h})$ onto $P_{k}(K_{h}),$ $i.e$ . $div(RT_{k}(K_{h}))=P_{k}(K_{h})$ .

210
Proposition 2. For $k\geq 0$ and any $q\in RT_{k}(K_{h})$ , the following relations imply $q=0$ .
$\int_{\partial K_{h}}q\cdot\tilde{n}\varphi_{k}ds=0$$\int_{K}..q\cdot q_{k-1}dx=0$

The Raviart-Thomas finite element space$RT_{k}$$\{p_{h}\in(L^{2}(\Omega))^{2}$

$;=$

:

$RT_{k}$,,$\forall\varphi_{k}\in R_{k}(\partial K_{h})$$\forall q_{k-1}\in(P_{k-1}(K_{h}))^{2}$

A residual bound with

$RT_{k}$$a_{k},$ $b_{k},$$c_{k}\in P_{k}(K_{h})$,is continuous on the inter-element boundaries.

It is a finite dimensional subspace of $H(div, \Omega)$ . Further let us define
It follows $div(RT_{k})=M_{h}$ (cf. Chapter IV.1 of [11]).

2.4.

is given by$p_{h}|\kappa_{h}=(\begin{array}{l}a_{k}b_{k}\end{array})+c_{k}\cdot(\begin{array}{l}xy\end{array}),$$p_{h}\cdot n$,$M_{h}$$:=\{v\in L^{2}(\Omega) :$\}$v|_{K_{h}}\in P_{k}(K_{h})\}$.

element

For the residual bound estimation, the smoothing technique in Subsection 2.2 works well to give accuratebounds. Some general smoothing techniques have been proposed in [2, 4, 13], etc, where smoothing functions
or $H(div, \Omega)$ are often used. One feature of proposal method is that we can use the basic
property of Raviart-Thomas element, $div(RT_{k})=M_{h}$ , for getting effective residual estimation. For given$f_{h}\in M_{h}$ , this property enbables us to define a subspace ofas$p_{h}\in(H^{1}(\Omega))^{2}$$RT_{k}$$W_{fh}=$ $\{ p_{h}\in RT_{k}:divp_{h}+f_{h}=0\}$.

Furthermore, we define $v_{h}\in M_{h}$ by an orthogonal projection of $v\in L^{2}(\Omega)$ such that $(v-v_{h}, w_{h})=0$ ,for $v_{h}\in M_{h}$ is obtained. Also we define. Assuming an error estimate. Finally, inequalities (10) and (11) give the following evaluation of theby the projection of,residual bound using$\forall w_{h}\in$$f_{h}(\hat{u})\in M_{h}$$\Vert v-v_{h}||_{L^{2}}\leq C_{M_{\hslash}}\Vert v\Vert_{V}$$M_{h}$$f(\hat{u})\in L^{2}(\Omega)$$p_{h}\in W_{fh(\hat{u})}$$\Vert \mathcal{F}(\hat{u})\Vert v$$=$$\sup_{0\neq v\in V}\frac{|(\nabla\hat{u},\nabla v)-(f(\hat{u}),v)|}{\Vert v||_{V}}$$=$$\sup_{0\neq v\in V}\frac{|(\nabla\hat{u}-p_{h},\nabla v)+(p_{h},\nabla v)-(f(\hat{u}),v)|}{\Vert v\Vert_{V}}$$\leq$$\sup_{0\neq v\in V}\frac{|(\nabla\hat{u}-p_{h},\nabla v)|}{\Vert v\Vert_{V}}+\sup_{0\neq v\in V}\frac{|(divp_{h}+f(\hat{u}),v)|}{\Vert v\Vert_{V}}$$\leq$$\Vert\nabla\hat{u}-p_{h}\Vert_{L^{2}}+\sup_{0\neq v\in V}\frac{|(divp_{h}+f_{\hslash}(\hat{u})+f(\hat{u})-f_{h}(\hat{u}),v)|}{\Vert v\Vert_{V}}$$=$$\Vert\nabla\hat{u}-p_{h}\Vert_{L^{2}}+\sup_{0\neq v\in V}\frac{|(f(\hat{u})-f_{h}(\hat{u}),v-v_{h})|}{\Vert v\Vert_{V}}$$\leq$$\Vert\nabla\hat{u}-p_{h}\Vert_{L^{2}}+C_{M_{h}}\Vert f(\hat{u})-f_{h}(\hat{u})\Vert_{L^{2}}$

.

(18)

Remark 3. Proposed estimation (18) holds for $k\geq 0$ . If the approximate solution is obtained from ,which has member function to be piecewzse $(k+1)$ -th polynomial. An effective choice of functional spacespanned byelements. The $mte$ of convergence can beandis subspace ofus to choose
and $\Vert f-f_{h}\Vert_{L^{2}}=o(h^{k+1})$ .
expect to be$\hat{u}$$W_{fh}$$W_{fh}$$RT_{k}$$M_{h}$$V_{h}$$P_{k}$$\Vert\nabla\hat{u}-p_{h}\Vert_{L^{2}}=o(h^{k+1})$

3

Computational result

Now we will present numerical results to illustrate our method. All computations are carried out on Mac OS
X 10.6.7,
GHz Quad-Core Intel Xeon (Westmere) with $64GB$ RAM by using MATLAB 20lla with
$2\cross 2.4$

211
a toolbox for verified computations, INTLAB [16]. We use Gmsh [17] http: $//geuz.org/gmsh/)$ to obtain
triangular mesh. Let us treat the following model problem. Here, is assumed to be hexagonal domain,$($$\Omega$$\{\begin{array}{ll}-\triangle u=u^{2}+10, in \Omega,u=0, on \partial\Omega.\end{array}$

There are two approximate solutions
given by finite element method. These are displayed in
Figure 1, 2 with the mesh size
. For the first approximate solution, verification results are shown in
Table 1, 2. Here, we use(piecewise linear) and(piecewise quadratic) elements for getting . We
space for -element andadopt
space for -element.
Comparing two cases in Table 1 and Table 2, we can observe that higher order elements yield improved
result.
$\hat{u}_{1},\hat{u}_{2}\in V_{h}$$2^{-4}$$\hat{u}_{1}$$P_{1}$$RT_{0}$$P_{2}$$RT_{1}$$P_{1}$

Figure 1:

Table 1:$\hat{u}_{1}$$\hat{u}_{1}$:(mesh size

$P_{1},$$\frac{1}{16}$$\hat{u}_{1}$$P_{2}$

Figure 2:)

Table 2:$p_{h}\in RT_{0}$$\hat{u}_{2}$$\hat{u}_{1}$:(mesh size

$P_{2},$$\frac{1}{16}$)$p_{h}\in RT_{1}$Next, we present results with respect tofinite element space. In Table 3, comparison
which is fromof each evaluation (14), (17) and (18) implies our proposed one works well. Numeric values on last column inTable 3 express upper bound of absolute error using (18) residual bounds. Based on Newton-Kantorovich
theorem, we prove that there is a solution in.
$P_{2}$$\hat{u}_{2}$$\rho$$\overline{B}(\hat{u}, \rho)$

Table 3: Residual evaluations for$\hat{u}_{2}$212A

Notes of Raviart-Thomas elements on triangle

In this part, we would like to note representations of the lowest $(RT_{0})$ and lst order $(RT_{1})$ Raviart-Thomas,
are numbered as 1, 2, 3. Their coordinates are. Vertices of
element on a triangle element. Let us denote $a_{i}=x_{j}y_{k}-x_{k}y_{j},$ $b_{i}=y_{j}-y_{k}$ , ci xk–xj where $(i,j, k)$ are eveninpermutation of (1, 2, 3). Here, we put subtense of each vertex as with direction from to . See
Figure 3. Then it follows$(x_{2}, y_{2}),$$(x_{1}, y_{1})$$K_{h}$$K_{h}$$=$$(x_{3}, y_{3})$$j$$e_{i}$

3

Figure 3: Triangle elements$|e_{i}|=(b_{i}^{2}+c_{i}^{2})^{1/2}$

,

$D=$

1
1
$1$

$1$

Furthermore, the unit normal vector$n_{i}$$K_{h}$$x_{1}$$y_{1}$$x_{2}$$y_{2}$$x_{3}$$y_{3}$

and$k$$K_{h}$$(X, Y)$$\tilde{K}_{h}$$=b_{j}c_{k}-b_{k}c_{j}$.

on each side is given by,$n_{i}=(n_{2}^{(i)}n_{1}^{(i)})= \frac{-\sigma}{|e_{i}|}(\begin{array}{l}b_{i}c_{i}\end{array})$

where $\sigma=D/|D|$ is corresponding to the direction of numbering. Namely,$\sigma=\{\begin{array}{ll}1, ( i,j, k :counter clockwise rotation),-1, ( i,j, k :clockwise rotation).\end{array}$

For $q\in RT_{k}(K_{h})$ , degrees of freedom are given by$\int_{\partial K_{h}}q\cdot n\varphi_{k}ds$$\int_{K_{h}}q\cdot q_{k-1}ds$

A.l$RT_{0}$,,$\varphi_{k}\in R_{k}(\partial K_{h})$for $k\geq 0$ ,,$q_{k-1}\in(P_{k-1}(K_{h}))^{2}$,

for $k\geq 1$ .

(19)

(20)

element

For $Ph\in RT_{0}$ , the representation of$RT_{0}$element$p_{h}$

on a triangle$K_{h}$

is given by$p_{h}|_{K_{h}}=(\begin{array}{l}\alpha_{1}\alpha_{2}\end{array})+\alpha_{3}(\begin{array}{l}xy\end{array})$

Let us explain how to determine coefficientsequivalent to (19) in case of $k=0$ .

$\alpha_{i}$. Three freedoms are given by the following form, which is$\gamma_{i}=|e_{i}|p_{h}\cdot n_{i}$213

Figure 4: $RT_{0}(K_{h})$Notice that$[n_{1}^{(3)}n_{1}^{(2)}n^{(1)}1$$p_{h}\cdot n_{i}=p_{h}|_{(x_{j},y_{j})}\cdot n_{i}$$n_{2}^{(3)}n_{2}^{(2)}n_{2}^{(1)}$Figure 5: $RT_{1}(K_{h})$, we have$x_{1}n_{1}^{(3)}+y_{1}n_{2}^{(3)}x_{3}n_{1}^{(2)}+y_{3}n_{2}^{(2)}x_{2}n^{(1)}1+y_{2}n_{2}^{(1)}]\{\begin{array}{l}\alpha_{1}\alpha_{2}\alpha_{3}\end{array}\}=[\gamma_{3/}\gamma_{2}\gamma_{1/}/|\begin{array}{l}e_{1}e_{2}e_{3}\end{array}|]$$\Leftrightarrow\sigma\{\begin{array}{lll}-b_{1} -c_{1} a_{1}-b_{2} -c_{2} a_{2}-b_{3} -c_{3} a_{3}\end{array}\}\{\begin{array}{l}\alpha_{1}\alpha_{2}\alpha_{3}\end{array}\}=\{\begin{array}{l}\gamma_{1}\gamma_{2}\gamma_{3}\end{array}\}$.

Using facts for $i=1,2,3$ ,$\{\begin{array}{ll}\sum a_{i}=D, \sum b_{i}=\sum c_{i}=0,\sum b_{i}x_{i}=D, \sum a_{i}x_{i}=\sum c_{i}x_{i}=0,\sum c_{i}y_{i}=D, \sum a_{i}y_{i}=\sum b_{iy_{i}}=0,\end{array}$and $\sigma D=|D|$ , we have$\alpha_{1}=-\frac{\sum\gamma_{i}x_{i}}{|D|}$

Therefore,$RT_{0}$element on$K_{h}$,$\alpha_{2}=-\frac{\sum\gamma_{i}y_{i}}{|D|}$,can be expressed with freedoms$\alpha_{3}=\frac{\sum\gamma_{i}}{|D|}$.

$\gamma_{i}$,$p_{h}|_{K_{h}}= \sum_{i=1}^{3}\frac{\gamma_{i}}{|D|}(\begin{array}{l}x-x_{i}y-y_{i}\end{array})=\sum_{i=1}^{3}\gamma_{i}\psi_{i}$where$\psi_{i}$are base functions offinite element space.

$RT_{0}$Remark 4. The image of $RT_{0}(K_{h})$ is given in Figure 4. Further for $q\in(L^{2}(\Omega))^{2}$ , let us define linear
functional, $F_{i}(q)=$ I $e_{i}|\{q(x_{j}, y_{j})\cdot n_{i}\}(i=1,2,3)$ . It follows$1\leq i,j\leq 3$$F_{i}(\psi_{j})=\delta_{ij}=\{\begin{array}{l}1, (i=j),0, (i\neq j),\end{array}$A.2$RT_{1}$.

element

Next let us consider lst order Raviart-Thomas finite element. Degrees of freedom are denoted by. For simplicity, we will transform triangle, which has vertices $(0,0),$ $(h, 0),$ $(X, Y)$
toin Figure 3.

$\gamma_{i}\in$$\mathbb{R}(i=1, \ldots, 8)$$\tilde{K}_{h}$$K_{h}$$h=(b_{3}^{2}+c_{3}^{2})^{1/2}$,$(\begin{array}{l}XY\end{array})=\frac{1}{h}(\begin{array}{ll}c_{3} -b_{3}b_{3} c_{3}\end{array}) (\begin{array}{l}-c_{2}b_{2}\end{array})$$n_{1}= \frac{\sigma}{|e_{1}|}(\begin{array}{l}Y-(X-h)\end{array})$In the following, we would like to exPlain,$RT_{1}$$n_{2}= \frac{\sigma}{|e_{2}|}(\begin{array}{l}-YX\end{array})$

element on$\tilde{K}_{h}$.

,$RT_{1}$,$D=hY$,$n_{3}= \frac{\sigma}{|e_{3}|}(\begin{array}{l}0-h\end{array})$

element$p_{h}$$p_{h}|_{K_{h}^{-}}=(\begin{array}{l}\alpha_{1}+\alpha_{2}x+\alpha_{3}y\alpha_{4}+\alpha_{5}x+\alpha_{6}y\end{array})+(\alpha_{7}x+\alpha_{8}y)(\begin{array}{l}xy\end{array})$.

is represented on.

$\tilde{K}_{h}$,214
Coefficientsare obtained by the following method of determination with respect todegrees of freedom are given by (19) and (20),$\alpha_{i}$$\int_{e}.Ph$

where$\phi_{j},$ni$\phi_{k}$$\phi_{j}ds=\gamma_{i}$,$\int_{e}.Ph$

ni $\phi_{k}ds=\gamma_{i+3}$ ,$\int_{K_{h}^{-}}Ph$denote piecewise linear functions on. So that we have$e_{i}$$(\begin{array}{l}l0\end{array})ds=\gamma_{7}$,$\gamma_{i}$. For $i=1,2,3$ ,$(\begin{array}{l}0l\end{array})ds=\gamma s$$\int_{K_{h}^{-}}P\hslash$,, satisfying $\phi_{j}(x_{j}, y_{j})=\phi_{k}(x_{k}, y_{k})=1,$ $\phi_{j}(x_{k}, y_{k})=$$\phi_{k}(x_{j}, y_{j})=0$$\frac{\sigma}{6}|$$-3Y-3Y3_{0}^{0}Y3YO6$$Y(2X+h)Y(X+2h)2(X_{O}^{O}+h)-2_{0}XY-XY$$-2Y^{2}-Y^{2}2Y^{2}2Y\gamma_{0}^{2}00$$-3(X–3(\overline{x_{0}^{3h}-}3h3_{6}X^{-}3Xh)h)$$-(x_{2x_{4^{2X+h)}}^{2}}-h)(X+2h)-t^{x-h)}2(X^{O}+h)-2h^{2}-h^{2}X$$(2X+^{O}h)Y/22h_{0}^{O}Y^{2}hY^{2}Y^{2}0$$-2(x_{0}^{0}-h)Y-(Xh)Y2XYXY2Y0$$h^{2}h^{0}(2^{+x+X^{2}}X+^{0}h)Y/2hY(2x+h)hY(x_{0}+2h)0$

Solving above linear system, we have the value of each coefficients. Then,$p_{h}|_{K_{h}^{-}}= \sum_{i=1}^{8}\gamma_{2}\psi_{i}$where$\psi_{i}$$\ldots,$$8)$$|\begin{array}{l}a_{1}\alpha 2\alpha 3\alpha 4a_{5}\alpha 6a_{7}\alpha 8\end{array}|$$=|\begin{array}{l}\gamma_{l}\gamma_{2}\gamma_{3}\gamma_{4}\gamma_{5}2_{77}/D\gamma_{6}2\gamma_{8}/D\end{array}|$

element is described on.

$\overline{K}_{h}$,,are base functions as following

Remark 5. See Figure 5 for degrees of freedom to, such that
for
$1,$

$RT_{1}$$|$$RT_{1}(\tilde{K}_{h}).$A linear functional is defined by$F_{i}(q),$$(i=$$q\in(L^{2}(\Omega))^{2}$$fi(q)=$$q\cdot n_{l}\phi_{m}ds$,$f\iota_{+3}(q)=\int_{e_{l}}q\cdot n_{l}\phi_{n}\$,$F_{7}(q)= \int_{K_{h}^{-}}q\cdot(\begin{array}{l}l0\end{array})dx$,$F_{8}(q)= \int_{K_{h}^{-}}q\cdot(\begin{array}{l}01\end{array})dx$

51

where $(l, m, n)$ are even permutation of (1, 2, 3). Then, we have$F_{i}(\psi_{j})=\delta_{ij},$$(1\leq i,j\leq 8)$.

215
References[1] M. T. Nakao, A numerical approach to the proof of existence of solutions for elliptic problems, Japan
Journal of Applied Mathematics, 5 (1988), pp.313-332.
[2] M.T. Nakao and Y. Watanabe, Numerical verffication methods for solutions of semilinear elliptic boundary value problems, NOLTA, IEICE, Vol.E94-N, No. 1 pp.2-31, 2011.

[3] M. Plum, Explicit-estimates and pointwise bounds for solutions of second-order elliptic boundary
vaJue problems, Journal of Mathematical Analysis and Applications, 165 (1992), pp.36-61.
$H_{2}$[4] M. Plum, Computer assisted proofs for semilinear elliptic boundary value problems, Japan Journal ofIndustrial and Applied Mathematics, 26 (2009), pp.419-442.
[5] P. Deuflhard and G. Heindl, Affine Invariant Convergence Theorems for Newton Method and Extensions to Related Methods, SIAM Journal on Numerical Analysis, vol.16, no.1, pp. 1-10, February
$s$

1979.

[6] L. V. Kantorovich and G.P.Akilov, Functional Analysis in Normed Spaces, translated from the Russian
by D. E. Brown, Pergamon Press, (1964).
[7] P. Grisvard, Elliptic Problems in Nonsmooth Domain, Pitman, Boston, (1985).
[8] M.T. Nakao, T. Kinoshita, Some remarks on the behaviour of the finite element solution in nonsmoothdomains, Applied Mathematics Letters 21, 12 (2008), pp.1310-1314.

[9] P.A. Raviart and J.M. Thomas, Introduction \’a 1‘Analyse Num\’erique des Equations aux D\’erivees Partielles, Masson, (1983).
[10] D. Braess, Finite elements-Theory, fast solvers, and applications in solid mechanics-, Third Edition,
Cambride University Press, 2007.
[11] F. Brezzi and M. Fortin, Mixed and Hybrid Finite Element Methods, Springer-Verlag, (1991).

[12] X. Liu and S. Oishi, Verified eigenvalue evaluation for elliptic operator on arbitrary polygonal domain,
prepare to publication.

[13] N. Yamamoto and M.T. Nakao, Numerical verifications for solutions to elliptic equations using residual
equations with a higher order finite element, Journal of Computational and Applied Mathematics, 60(1995), pp.271-279.
[14] A. Takayasu and S. Oishi, A method of computer assisted proof for Nonlinear two-point boundary value
problems using higher order finite elements, NOLTA, IEICE, Vol.E94-N, No.1 (2011), pp.74-89.
[15] F. Kikuchi and H. Saito, Remarks on a posteriori error estimation for finite element solutions, Journalof Computational and Applied Mathematics, 199 (2007), pp.329-336.
[16] S.M. Rump, INTLAB-INTerval LABoratory, in Tibor Csendes, editor, Developments in Reliable
Computing, Kluwer Academic Publishers, Dordrecht (1999), pp.77-104.

[17] Gmsh: a three-dimensional finite element mesh generator with built-in pre- and post-processing facilities, ht tp $://www$ . geuz. $org/gmsh/$