数理解析研究所講究録
1169 巻 2000 年 27-56

27

Numerical Verification Methods for Solutions of Ordinary and Partial
Differential Equations
Mitsuhiro T. Nakao
中尾充宏

Graduate School of Mathematics, Kyushu University 33
Fukuoka 812-8581, Japan-mail: mtnakao@math.kyushu-u.ac.jp$\mathrm{e}$Abstract.
In this article, we describe on a state of the art of validated numerical computations for solutionsof differential equations. A brief overview of the main techniques in self-validating numerics
for initial and boundary value problems in ordinary and partial differential equations including
eigenvalue problems will be presented. A fairly detailed introductions are given for the author’s
own method related to second-order elliptic boundary value problems. Many references which
seem to be useful for readers are supplied at the end of the article.
Key words.
analysis1

Numerical verification, nonlinear differential equation, computer assisted proof in

Introduction
If we denote an equation for unknown

$u$

by

$F(u)=0$(1.1)

then the problem of finding the solution generally implies an dimensional simultaneous system of equations provided is in some finite dimensional ( -dimension) space.
Since very large scale, e.g. several thousand, linear system of equations can be easily
solved on computers nowadays, it is not too surprising that, even if (1.1) is nonlinear,
we can verify the existence and uniqueness of the solution as well as its domain by somenumerical computations on the digital computer. However, in the case where (1.1) isa differential equation, it becomes a simultaneous equation which has infinitely many
unknowns because the dimension of the potential function space containing is infinite.
Therefore, we might naturally feel that it is impossible to study the existence or uniqueness of the solutions by finite procedures based upon the computer arithmetic. Actually,
when the author first ran across an assertion about the possibility of such argumentsby computer, in Kaucher-Miranker [13] about fifteen years ago, he could not believe itand seriously attempted to find their theoretical mistakes. When itbecame clear thatthe grounds of their arguments could not be denied, though the application area seemed
rather narrow, the author was greatly impressed and was also convinced that such a
study must be one of the most ’portant research areas in computational mathematics,
particularly in numerical analysis. Subsequently, the author has learned that there had
been similar researches centered around ordinary differential equations(ODEs). However,with his research background it was natural for this writer to be particularly interested
in numerical verification methods for partial differential equations(PDEs).
$u$

$u$

$n$

$\mathrm{n}$$\mathrm{u}$$\mathrm{i}\mathrm{m}$28
In this exposition, we will survey the state of the art, fairly emphasizing the author’s
own works, about the verification methods for the existence, uniqueness and enclosure ofsolutions for differential equations based on numerical computations.
Many differential equations appearing in the mathematical sciences such as physics ortechnology are numerically (approximately) solved by the use of contemporary supercomputers, and those computed results are supplied for simulations of phenomena, indeuniqueness of solutions. Naturally, there
pendently of the guarantee of existence
are not a few theoretical studies done by mathematicians for such equations. However,
because of the nonlinearity or variety of the problem and so on, applications of unified
theory are quite difficult and it is hard to say that those results sufficiently satisfy demands of the people who are working on the numerical analysis in various areas. Also
even if the existence and uniqueness are already known, in general, we cannot confirm
the order of magnitude of the difference between the computed solutions and the exact
solutions. Actually in such a case, we are obliged to be satisfied with the ambiguous validity which may be interpreted as the comparison with experimental data. On the other
hand, in pure mathematics, the problem proving the existence of solutions for particulardifferential equations can often arise. The principal purpose of numerical verification
methods here is the mathematically and numerically rigorous grasp of solutions of differential equations appearing in various mathematical sciences including pure mathematics.
Therefore, we note that the term ’numerical’ does not mean ’approximate’. Although it
is not too long time since this kind of research was started, a great number of developments will be expected in the 21st century as a new approach in numerical analysis which
exceeds the existing numerical methods in the sense of assurance of numerical qualities
for infinite-dimensional problems.
In the followings, in Section 2, we describe the case studies of numerical verification
methods for ODEs, stressing on Lohner’s method for initial value problems. And, in
Section 3, the methods for partial differential equations for elliptic and evolutional problems, around the author’s method, will be surveyed. Next, in section 4, we will treat the
eigenvalue problems of elliptic operators. Finally, we will give some concluding remarks
in Section 5.
$\mathrm{a}\mathrm{n}\mathrm{d}/\mathrm{o}\mathrm{r}$2Ordinary differential equations(ODEs)

A germ of numerical approaches to the verification of solutions for ODEs already appeared in Cesari [7] in the first half of 1960, and since then not a few of case studies havebeen done. Generally speaking, the functional equation is equivalent to the simultaneos
equations with infinite dimension. Therefore, it should be the essential point to consider
the errors caused by the truncation or discretization of the original problems. The enclosure methods for function space problems, in the author’s opinion, will be classified into
the following three groups according to their verification techniques.
(l)Analytic method:

This is a method such that, by using the estimations of constants, norms of functions
including approximate solutions, and operators appearing in the equation, one prove

29
that they satisfy a certain condition, e.g., the convergence condition in the NewtonKantorovich theorem in the below. Particularly, the norm estimation of an inverse linearlized operator for the original problem is most important and essential task in thismethod. The interval analysis is used only in the simple arithmetic calculations.
(2) Interval method:
An interval can be considered as the set of functions whose ranges are contained thatinterval. For example, $[a, b]$ can be considered as the set:$\{f\in C(\alpha, \beta)|f.(x)\in[a, b], \forall x\in(\alpha, \beta)\}$.

In this method, using this infinite dimensional property of the interval, the truncation
errors are essentially grasped by the interval arithmtic. In order to apply this method,usually, we need a transformation of the original differential equation to an equivalent
integral equation, but there are no norm estimations of the inverse linearlized operator.
(3) Mixed method:
This is an intermediate method between the above two. The interval plays an essential
role in this method, but the analytic arguments in a certain function space as well.
Although some kind of inverse linearlized operator is used, in general, but it is notdirectly evaluated as an infinite dimensional operator.
We now introduce some of the typical studies of these categories.

Analytic methods The following Newton-Kantorovich theorem is used for the exis(e.g., [87] for proof).
tence and local uniqueness of the solution of functional$\mathrm{e}\mathrm{q}\mathrm{u}\mathrm{a}\mathrm{t}\mathrm{i}_{0}\mathrm{n}\mathrm{s}$Theorem 2.1 Let $X,$

be Banach spaces, and let $D\subset X$ be a convex set. For a map$F:D\subset Xarrow Y,$ $F’,$ $F”$ denote the Fr\’echet derivatives. Assume that, for some $u_{0}\in D$and ,and constants $B,$(i) $[F’(u_{0})]^{-1}$ exists, andand $||[F’(u_{0})]^{-1}F(u_{0})||\leq\eta$ .
(ii) $||F’’(u)||\leq\kappa,$ $||u-u0||\leq r$ .
(iii) $h=\eta\cdot B\cdot\kappa<1/2$ .
(iv) $r\geq r_{0}=\eta(1-\sqrt{1-2h})/h$ .
Then, the equation $F(u)=0$ has one and only one solution
in the ball: $||u-u0||\leq r_{0}$ .
$\eta,$$\kappa$$Y$$r$$||[F^{;}(u\mathrm{o})]^{-1}||\leq B$$u^{*}$

applied this theorem to the verification of solutions for nonlinear two pont boundary value problems. On the other hand, in Japan, Urabe [78], [79]developed a useful version of the Kantorovich-like theorem for the verification of solutions for multipoint and periodic boundary value problems. Plum’s method [59] is alsoclassified into this category. In his method, the norm of the inverse linearized operator isbounded by using eigenvalue enclosing technique with homotopy
also Section
3 and 4). Oishi [55] derived a version of the Theorem 2.1 by using some finite dimensional
projection and its a priori error estimates, and proved, as an example, the existence of
periodic solutions for the Duffing equation. Also, in [57], his method was extended forthe finite element projection in one dimension.
$\mathrm{K}\mathrm{e}\mathrm{d}\mathrm{e}\mathrm{m}[15],$$\mathrm{M}\mathrm{C}\mathrm{c}_{\mathrm{a}\mathrm{r}}\mathrm{t}\mathrm{h}\mathrm{y}[18]$$\mathrm{m}\mathrm{e}\mathrm{t}\mathrm{h}_{0}\mathrm{d}(\mathrm{s}\mathrm{e}\mathrm{e}$30Sinai’s work [74] on the verification of the existence of a closed orbit is based upon the
fixed point problem of a Poincar\’e map and the numerical proof of a Kantorovich typeconvergence condition. He verified a closed orbit of the following Lorenz model :

$\{$

$x’$

$=$

$a_{1}x+b1yz+b1^{XZ}$

$y’$

$=$

$a_{2}y-b_{1}yz-b1xz$

$z’$

$=$

$-a_{3}z+(x+y)(b_{2^{X}}+b_{3y})$(2.1)

Also, Nishida [52] solved, applying the theory of pseudo trajectory similar in [74], somebifurcation problems appeared in nonlinear fluid dynamics by numerical enclosing technique of the critical eigenvalues of linearized problem.

Interval methods In Moore’s initial work [24], a basic idea of this method was alreadypresented. Such a primitive technique was, however, not so effective for the practicalproblems. Nowadays, by the view point of the verification principle, the most famous
which we describe an outline in the below. In
method will be Lorner’s
denote the set of real intervals and -dimensional interval
what follows, let $IR$ andvectors, respectively.
Consider the following initial value problem of autonomous system:$\mathrm{t}\mathrm{e}\mathrm{c}\mathrm{h}\mathrm{n}\mathrm{i}\mathrm{q}\mathrm{u}\mathrm{e}[17]$$n$$IR^{n}$(2.2)$u’=f(u)$Here, $u=u(t)\in R^{n},$

$f$

is a

$C^{p}$

map on$R^{n}$

for a fixed positive integer , and the initial$p$value is given as$u(0)\in[u_{0}]=u_{0}+[z_{0}]$

(2.3)

.

point vector. Note that, in case of nondenote the interval vectors,
Here,autonomous, according to add a new dependent variable $u_{n+1}=t$ and an equation$u_{n+1}’=1$ , it can be reduced to the type (2.2). Our purpose is, for fixed $T>0$ , to get aninterval enclosure $[u(t)]$ for a solution of (2.2), (2.3) such that$[u_{0}],$ $[z_{0}]$$u_{0}$$u$$u(t)\in[u(t)]$ ,$t\in[0, T]$.

First, we consider an explicit one step method for (2.2), (2.3) with step size based, where $u(t)$ stands for an exact
andon a function $\Phi=\Phi(u)$ . Settingby $z_{j+1}/h$ ,solution. When we denote the local discretization error between and$h$$u_{j}\equiv u(t_{j})$$t_{j}\equiv jh$$t_{j+1}$$t_{j}$

it holds that$u_{j+1}=u_{i}+h\Phi(uj)+Z_{j}+1$ ,$j\geq 0$

(2.4)

,

can be estimated by and its
where $u_{0}\in[u_{0}]$ . Suppose that is chosen so thatderivative. For example, if the scheme (2.4) is order , then we have, using the Taylor
expansion of $u_{j+1}=u(t_{j}+h)$ at ,$\Phi$$u$$z_{j+1}$$p$$t_{j}$$z_{j+1}= \frac{h^{p}}{p!}u^{(p)}(\hat{t}_{j+1})$,$\hat{t}_{j+1}\in(t_{j}, t_{i+1})$

(2.5)

.

by using (2.5).
is obtained, we can encloseatWhen an assured interval
But such a simple method would yield the monotone and rapid increasing of the width$[u_{j}]$$t_{j}$$u_{j+1}$

31
of the assured interval. The basic idea of Lohner’s method consists in the representation
as a function of the -dimensional vector variable. That is, theofbyisdetermined
all
the
local
up
value
at
discretization
errorsexact
to instead of
only.
and
If is continuously differentiable with respect to $j+2$ variables $z0,$ $z_{1,j+}\ldots,$$Z1$ , then$u_{j+1}=u_{j+1}(z_{0}, Z_{1}, \cdots, z_{j}+1)$ as well. Therefore, by using the mean value theorem arounding representation(so, $s1,$$sj+1$ ), we have the$n$$u_{j+1}$$z_{0},$$Z_{1},$$\cdots,$$z_{j+}1$$j$$t_{j+1}$$z_{j+1}$$u_{j}$$\Phi$

follO.w$\cdots,$$u_{j+1}= \overline{u}_{j+1}+\sum_{k=0}^{j+1}\frac{\partial u_{j+1}(z\wedge)}{\partial z_{k}}(z_{k}-s_{k})$

,

(2.6)

where is usually chosen as the midpoint of . And, is an unknown vector which is deand $s=(s_{0}, s_{1}, \cdots , s_{j+1})$ . Since $\overline{u}_{j+1}=u_{j+1}(s0, s1, \cdots, sj+1)$ ,termined by $z=$
we have$z\wedge$$z_{k}$$s_{k}$$(z_{0}, z_{1}, \cdots , z_{j+1})$$\tilde{u}_{j+1}=\overline{u}_{j}+h\Phi(\overline{u}_{j})+s_{j+1}$,$j\geq 0$

,

(2.7)

where,
and $s_{0}=0$ .
Next, we briefly mention about an interval enclosing algorithm of (2.6). Assume that
we have already got the enclosure until .
$\overline{u}_{0}=u0$$t_{j}$

for $u(t)$ on.
1. First step : Calculation of a rough enclosureThis can be done by computing a constant interval which encloses the solution ofof the form:
the equivalent integral equation on$[t_{j}, t_{j+1}]$$[u_{j+1}^{0}]$$[t_{j}, t_{j+1}]$$u(t)=u_{j}+ \int_{t_{j}}^{t}f(u(s))d_{S}$,$u_{j}\in[u_{j}],$ $t\in[t_{j}, t_{j+1}]$

,

(2.8)

implies the interval enclosure of the solution at $t=t_{j}$ .
where, we defineasFor an interval $X$ such that$[u_{j}]$$Y$$[u_{j}]\subset X$$Y\equiv[u_{j}]+[0, h]\cdot f(X)$

.

(2.9)

If $Y\subset X$ then, by Schauder’s fixed point theorem, it holds that,$u(t)\in Y$

Therefore, we set
small $\in>0$ , setting$[u_{j+1}^{0}]\equiv Y$$\forall t\in[t_{j}, t_{j+1}]$.

. If $Y\subset X$ does not hold, then, for an appropriately(2.10)$X:=(1+\epsilon)Y-\epsilon Y$( -inflation), we retry the computation (2.9) for this $X$ and check the same relation.
If we could not get the desired inclusion within the definite times, then we adopt a
smaller step size.
in the2. Second step: Calculate the local discretization error $[z_{j+1}]$ by the use of
first step.
3. Third step: Compute a new (narrow) enclosure $[u_{j+1}]$ by using $[z_{j+1}]$ .
obtained in the first step by a new
4. Forth step: Replace the rough enclosure$[u_{j+1}]$
repeatstep,thesecond-forth step until we get the
in
and
the
third
intervalinterval as small as possible.
$\epsilon$$[u_{j+1}^{0}]$$[u_{j+1}^{0}]$32
5. Fifth step: Determine $s_{j+1}\in[z_{j+1}]$ and calculate$\overline{u}_{j+1}$.

In [17], this method was applied to enclose a solution of the following Kepler-Ellipse
problem, by using Taylor’s method of order 18 and step size 0.1,,,

$u_{1}’$$=$$u_{3}$$u_{2}’$$=$$u_{4}$$u_{3}’$$=$$-u_{1}(u_{1}^{2}+u_{2})^{-1.5}2$$u_{4}’$

$=$

$\{$

(2.11),$-u_{2}(u_{12}^{2}+u2)^{-1}\cdot 5$ ,.
with initial condition: $u_{1}(0)=1.2,$ $u_{2}(0)=0,$ $u_{3}(0)=0,$
accuracy for each $t\in[0,70]$ . This method can
He got the enclosure of solution within
also be applied to the boundary value problems combining with the shoooting technique.
There are another works on the interval methods for initial value problems, e.g., [2],described a good survey of the techniques including[13], [75] etc. Especially,$u_{4}(0)=\sqrt{\frac{2}{3}}$$10^{-5}$$\mathrm{R}\mathrm{i}\mathrm{h}\mathrm{m}[68]$the fundamental ideas.
, which originally proposed for PDEs,
Mixed methods The author’s method
belongs to this category. The method uses both the constructive error analysis for thefinite element method and the interval coefficient functions of a finite dimensional space,which will be described in detail in the next section. [51] and [73], which uses someproperties of monotone operator, are regarded as another examples for the mixed method.
Oishi [54] proposed a kind of mixed method combining an approximate fundamental
matix with the interval representation of functions. He recently applied the method
to the numerical verification of the existence of solutions for a connecting orbit in thepresented an interesting computer assited. Mrozek etcontinuous
approach on the proof of a chaos property of the Lore.n model besed on the Conleyindex theory.
$(\mathrm{e}.\mathrm{g}.,[33])$$\mathrm{d}\mathrm{y}\mathrm{n}\mathrm{a}\mathrm{m}\mathrm{i}\mathrm{c}\mathrm{S}[56]$$\mathrm{a}1.[25]$$\mathrm{Z}$3Partial differential equations(PDEs)

There has been not so many works on the numerical verification for PDEs. As far
as the author is concerned, it was hard to find any methods except for Plum and theauthor’s own work up to recently. As mentioned before, the former is an analytic method
and the latter a mixed method. There is no interval method for PDEs, for it is difficult totransform a PDE to an equivalent integral equation. In the below, we present the outlineof both methods for second-order elliptic problems, particularly around the author’s
method. Moreover, since, quite recently, some case studies have appeared for computer
assited proof in the actual nonlinear analysis, we will also briefly refer to them.

33
3.1

elliptic problems

We consider the following nonlinear elliptic boundary value problem on a bounded$1\leq n\leq 3$ :inconvex domain$\Omega$$R^{n},$$-\triangle u$$=$$f(x, u, \nabla u)$$x\in$$\Omega$

.

$\{$

$=$

$u$

$0$

,

$x\in\partial\Omega$

(3.1)

,

where the map is assumed to satisfy appropriate conditions on the smoothness. For
an integer $m$ , let $H^{m}(\Omega)\equiv H^{m}$ denote -Sobolev space of order $m$ on . And set{ $\phi\in H^{1}|tr(\emptyset)=0$ on } with the inner product $<\phi,$, wheremeans the inner product on.
In the below, we denoteby. And define$f$$L^{2}$$\Omega$$H_{0}^{1}\equiv$$\partial\Omega$$\psi>\equiv(\nabla\phi, \nabla\psi)$$L^{2}(\Omega)$$(\cdot, \cdot)$$||\cdot|.|_{L^{2}}(\Omega)\equiv||\cdot||_{L^{2}}$$||\cdot||$$| \phi|_{H^{2}}^{2}\equiv\sum_{1i,,j=}^{n}||\frac{\partial^{2}\phi}{\partial x_{i}\partial_{X_{j}}}||^{2}L^{2}$

3.1.1

.

The author’s method

The verification principle of this method is first originated in 1988 by [29] and, in themeantime, several improvements have been done up to now.
In what follows, the map
in (3.1) is assumed to be continuous from the Sobolev
space
into
such that having a bounded image inon a bounded set
in. For example, when $n=2,$ $f(u)\equiv f(x, u, \nabla u):=g_{1}\cdot\nabla u+g_{2}u^{p}$ satisfies above
assumption, whereandare in, and an arbitrary nonnegative
integer. And, for $n=3$ , the same assumption holds for any such that $1\leq p\leq 3$ by the
Sobolev imbedding theorem.
Now letbe a finite dimensional subspace ofdependent on $h(0<h<1)$ . Usually,is taken to be a finite element subspace with mesh size . And let:
denote the-projection defined by$f$$H_{0}^{1}(\Omega)$$L^{2}(\Omega)$$L^{2}(\Omega)$$H_{0}^{1}(\Omega)$$g_{1}=(g_{1}^{1}, g_{1}^{2})$$L^{\infty}(\Omega)$$g_{2}$$p$$p$$(\mathrm{e}.\mathrm{g}., [1])$$S_{h}$$H_{0}^{1}$$S_{h}$$h$$H_{0}^{1}arrow S_{h}$$P_{h}$$H_{0}^{1}$$(\nabla u-\nabla(Phu), \nabla\hat{\emptyset})=0$,$\hat{\phi}\in S_{h}$We now suppose the following approximation property of
For any,$P_{h}$

.

(3.2)

.

$\phi\in H^{2}\cap H_{0}^{1}$$||\phi-P_{h}\emptyset||_{H_{0}^{1}}\leq C_{0}h|\emptyset|_{H^{2}}$

,

(3.3)

where
is a positive constant numerically determined and independent ofassumption holds for many finite elernent subspace ofdefined by piecewise$C_{0}$$h$. This$H_{0}^{1}$$\mathrm{p}\mathrm{o}\mathrm{l}\mathrm{y}\mathrm{l}\mathrm{n}\mathrm{o}\mathrm{m}\mathrm{i}\mathrm{a}\mathrm{l}\mathrm{s}$. For example, it can be taken as
andfor bilinear and biquadratic element, respectively, in case of the rectangular mesh [42].
For the triangular and uniform mesh of the domain in, we can take,$C_{0}=0.81$
for linear element [49], and for the more fine constant, see the arguments in the end of
this sub-subsection.
with quasi-uniform mesh$(\mathrm{e}.\mathrm{g}., [8],[53])$$\frac{1}{2\pi}$$C_{0}=\overline{\pi}$$R^{2}$$\mathrm{e}.\mathrm{g}.,$34
there exists a unique solution

Now, it is well known [11] that for arbitraryof the following Poisson’s equation:$\psi\in L^{2}(\Omega)$$\phi\in H^{2}\cap H_{0}^{1}$$-\triangle\emptyset$$=$$\psi$,$x\in\Omega$

,
(3.4)

$\{$

$\phi$

,

$0$

$=$

When we denote the solution of (3.4) byas well as the following estimate holds:$x\in\partial\Omega$$\phi\equiv K\psi$.

, the map $K$ :$|\phi|_{H^{2}}\leq||\psi||$$L^{2}arrow H_{0}^{1}$is compact

(3.5)

.

Defining the nonlinear map $F(u):=Kf(u),$ $F$ is a compact map on
following fixed point equation of the operator $F$ equivalent to (3.1):$H_{0}^{1}$and we get the(3.6)$u=F(u)$ .

Therefore, if we find a nonempty, bounded, convex and closed subset $U$ in$H_{0}^{1}$satisfying(3.7)$F(U)=\{F(u)|u\in U\}\subset U$ ,then by the Schauder fixed point theorem, there exists an element $u\in F(U)$ such that$u=F(u)$ . Usually, we choose such a set , which is referred a candidate set of solutions,. Here,stands for theandof the form $U=U_{h}\oplus U_{\perp}$ , where. Then, the verification condition can beinorthogonal complement subspace of$U$$U_{\perp}\subset S_{h}^{\perp}$$U_{h}\subset S_{h}$$S_{h}$$S_{h}^{\perp}$$H_{0}^{1}$

written as$P_{h}F(U)$$\subset$$U_{h}$$(I-P_{h})F(U)$$\subset$$U_{\perp}$

(3.8)

$\{$

.

Sometimes we call the quantities $R(F(U)):=PhF(U)$ and $RE(F(U)).–(I-P_{h})F(U)$and the rounding error of $F(U)$ , respectively.
as the rounding into
$S_{h}$

Then (3.8) implies that
(3.9)

$R(F(U))\oplus RE(F(U))\subset U$,

is taken to be a set ofwhich is the basic principle of our verification method. The set
a ball inwith interval coefficients, while
linear combinations of base functions in.
with radius
Namely,$U_{h}$$U_{\perp}$$S_{h}$$S_{h}^{\perp}$$\alpha\geq 0$$U_{h}= \{\phi_{h}\in S_{h}|\phi_{h}=\sum_{=j1}^{M}[\underline{A}j’ j\overline{A}]\phi j\}$(3.10)

and$U_{\perp}=\{\phi\in S_{h}^{\perp}|||\phi||_{H_{0}^{1}}\leq\alpha\}$,(3.11)

35
respectively, where$\{\phi_{j}\}_{j=1}M$

is a basis of$S_{h}$. Here,$\sum_{j=1}^{M}[\underline{A}_{j}, \overline{A}_{j}]\phi_{j}$is interepreted as the setof functions in which each element is a linear combinaion of
whose coefficient of
belongs to the corresponding interval
for each $1\leq j\leq M$ . We denote the set, that is,of such interval functions by$\{\phi_{j}\}_{j=1}^{M}$$[\underline{A}_{j},\overline{A}_{j}]$$\phi_{j}$$S_{h,I}$$S_{h,I}:= \{\hat{\phi}\in S_{h}|\hat{\phi}=\sum_{j=1}^{M}Aj\emptyset j, A_{j}\in IR\}$.

. Note that each element $\phi_{h}\in P_{h}F(U)$ satisfies
Then it can be considered as
the following finite element solution for some $\psi\in U$$S_{h}\subset S_{h,I}$$(\nabla\phi h, \nabla\hat{\emptyset})=(f(\psi),\hat{\emptyset})$,$\forall_{\hat{\phi}}\in S_{h}$.

(3.12)

Therefore, it can be eeasily seen that $P_{h}F(U)$ is directly computed or enclosed fromby solving a linear system of equations with interval right-hand side using
andsome interval arithmetic approaches, e.g., [3], [50]. On the other hand, $(I-P_{h})F(U)$ is
unknown but can be evaluated by the following constructive a priori error estimates forthe finite element solution of Poisson’s equation:$U_{h}$$U_{\perp}$$||(I-P_{h})F(U)||_{H_{0}^{1}} \leq C_{0}h\sup_{u\in U}||f(u)||$.

(3.13)

which is obtained by (3.3) and (3.5).
Thus, the former condition in (3.8) is validated as the inclusion relations of corresponding coefficient intervals, and the latter part can be confirmed by comparing two
nonnegative real numbers which correspond to the radii of balls. In the actual computation, we use an iterative method for both part of $P_{h}F(U)$ and $(I-P_{h})F(U)$ as below.
(1) Verification by the simple iteration method
As stated above, we usually find a candidate set of the form
$U=U_{h^{\oplus}}U\perp$

In the below, we fix an approximate solution

We consider the set of functions$U_{h}\in S_{h,I}$.

$\hat{u}_{h}$(3.14)of (3.1) such that,are intervals, in general, centered at . And the set
satisfying
Then we take an interval vector$0$$A_{i}$.

of the form

$U_{h}= \sum_{1i=}^{n}(u_{i}+A_{i})\phi_{i}$where$\hat{u}_{h}=\sum_{i=1}^{n}u_{i}\phi i\in S_{h}$(3.15)$U_{\perp}$

is same as (3.11).

$(B_{i})$$P_{h}F(U) \subset\sum_{i=1}^{n}B_{i}\phi i$,(3.16)

is usually determined by a solution of the linear system of equations withwhere
matrix $G:=((\nabla\phi_{i}, \nabla\phi_{j}))$ and the $M$intereval right-hand side. Namely, for the$(B_{i})$$M\mathrm{x}M$36
dimensional interval vector $\mathrm{b}:=((f(U), \phi_{i}))$ , the interval vector
as a solution of the following equation$G\cdot(B_{i})=\mathrm{b}$$(B_{i})$can be computed.

(3.17)

Here, $(f(U), \phi_{i})\in IR$ stands for the interval enclosure of the set $\{(f(u), \phi_{i})\in R|u\in U\}$ .

And we set$\beta$$:=C_{0}hu \sup_{\in U}||f(u)||_{L^{2}}$(3.18).

On the actual and detailed computational procedures for the determining the interval
vector and the estimation of the right-hand side in (3.18), refer [29], [30], [34] etc. Now
the computable verification condition is described as$\mathrm{b}$Theorem 3.1 For the sets defined by (3.14), (3.15) and (3.11), if the following conditions
hold$B_{i}\subset u_{i}+A_{i}$$\leq$$\beta$then there exists a solution$u$$\alpha$,$i=1,$$\cdots,$ $n$,(3.19),

of $u=F(u)$ in .
$U$

Based on Theorem 3.1, we obtain the following verification algorithm by using the(cf. [69]). In what follows, we define
simple iteration method withfor a nonnegative real number .
$[\alpha]\equiv\{\phi\in$$\delta- \mathrm{i}\mathrm{n}\mathrm{f}\mathrm{l}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}$$\alpha$$S_{h}^{\perp}|||\phi||_{H_{0}^{1}}\leq\alpha\}$Verification algorithm S-l$:=0$ , initial candidate set is defined by
1. Setting $A_{i}^{(0)}:=[0,0](i=1, \cdots, n)$ and.
and, computeanddetermined by
2. For the candidate set
from (3.17) and (3.18), respectively.
If (3.19) holds, then there exist a desired solution in the set$\alpha^{(0)}$$U^{(0)}:=\hat{u}_{h}$$(A_{i}^{(k)})$$U^{(k)}$$(B_{i}^{(k)})$$\alpha^{(k)}$$U^{(k)}= \hat{u}_{h}+\sum_{i=1}^{M}A\phi_{i}+i(k)[\alpha^{(})k]$$\beta^{(k)}$.

Otherwise, go to the next step.
3. Using some fixed small constant $\delta>0$ , after setting$A_{i}^{(k+1)}$$\alpha^{(k+1)}$

$:=$

$:=$

[-1, 1]$\delta+A_{i}^{(k)}$$(1+\delta)\beta^{(k})$,$i=1,$$\cdots,$ $n$,,

return to the previous step.

Remark 3.1 The above algorithm using Theorem 3.1 could work only for limited case.
However, we can say that not only the implememtation of the procedures is quite simple
and easy but also the essential point of our vemfication principle, $i.e.$ , the direct solving a fnite dimensional problem with addtional error estimates, is clearly shown in thisalgorithm.

37
Numerical example 1 ([30]).
The algorithm S-l has not so wide applications, because it needs that the map $F$
is retractive in some neighborhood of the fixed point to be verified. But, we actually
succeeded the verification for several realistic problems, e.g., as below [30]:$-\triangle u$$=$$f(x)u^{2}+g(x)$$.x\in\Omega.\equiv(0,1)\cross(0,1)$

$u$

$=$

$0$

$x\in\partial\Omega$,(3.20)$\{$,-functions on whose ranges are in [-2, 2] andwhere $f(x)$ and $g(x)$ are arbitrary$[0,7]$ , respectively. Here, as the finite element subspace, we used the bilinear elementson the uniform rectangular mesh with mesh size $h=1/15$ .
$L^{\infty}$$\Omega$$S_{h}$(2) Verification by Newton-like method
In order to apply our verification method for more general problems, we introduce akind of Newton-like method.
First, note that (3.6) can also be rewritten as the following decomposed form inand$S_{h}$$S_{h}^{\perp}:$$P_{h}u$$=$$P_{h}F(u)$$(I-P_{h})u$$=$$(I-P_{h})F(u)$(3.21)$\{$In order to consider the Newton type operator for (3.21), define the nonlinear operator$N$ onby$H_{0}^{1}(\Omega)$$N(u).–u-[P_{h}-P_{h}A’(\hat{u}_{h})]_{h}-1(P_{hh}u-PF(u))$ ,and ’ means the Fr\’echet derivative of at. Here,.
of the restriction operator
denotes the inverse on
The existence of such a finite dimensional inverse operator can be validated by the usual. Also noteinvertibility of a matrix corresponding to the restriction operator
approximate
operator.
by somethat we can replace
We now define
where$A’(\hat{u}_{h})\equiv(-\triangle)^{-1}f’(\hat{u}_{h})$$f$$[P_{h}-P_{h}A’(\hat{u}_{h})]_{h}-1$$\hat{u}_{h}$$(P_{h}-P_{h}A’(\hat{u}_{h}))|s_{h}$$S_{h}$$(\mathrm{e}.\mathrm{g}., [70])$$P_{h}A’(\hat{u}_{h})$$T(u)$$:=$$P_{h}N(u)+(I-P_{h})F(u)$ .

is considered as the Newton-like operator for the former part of (3.21) but
Thenis compact onbythe simple iterative operator for the latter part. Moreover,compactness of $F$ . Furthermore, it can readily be seen that
$T$

$T$

$H_{0}^{1}(\Omega)$

Proposition 3.1 The fixed point equation(3.22)$u=T(u)$is equivalent to (3.6).
Indeed, if (3.21) holds, then we have, by using the former part,$0=[P_{h}-P_{h}A’(\hat{u}_{h})]_{h}-1(P_{h}u-P_{h}F(u))$,(3.23)

38
which yields$P_{h}u=P_{h}N(u)$.

(3.24)

Therefore, $u=Tu$ follows by adding (3.24) to each side of the latter part of (3.21).
Conversely, if $u=Tu$ , then we immediately obtain (3.24), and thus (3.23) follows. The
conclusion would be now straightforward.
$\square$in
If we find a nonempty, bounded, convex and closed subset
satisfying$T(U)=\{T(u)|u\in U\}\subset U$ , then by the Schauder fixed point theorem there exists an,element $u\in T(U)$ such that $u=T(u)$ . If we choose a bounded set such as$U_{h}\in
S_{h,I}$, the verification condition can be written by
and
where
$U$

$H_{0}^{1}(\Omega)$$U$$U=U_{h}\oplus U_{\perp}$$U_{\perp}\subset S_{h}^{\perp}$,$P_{h}N(U)$$\subset$$U_{h}$$(I-P_{h})F(U)$$\subset$$U_{\perp}$$\{$

Notice that if we use the symbol rounding
represented as$R(\cdot)$(3.25).

and rounding error $RE(\cdot)$ , then (3.25) is(3.26)$R(\tau(U))\oplus RE(T(U))\subset U$

which is the corresponding relation to (3.9).
The computational procedure for $P_{h}N(U)(rounding)$ consists of the solving linear system of equations with interval right-hand side which is similar to that in the case of
simple iteration method. But concerned matrix is a Newton type one which is exactly
the same matrix as in the usual simplified Newton method for the discretized problem of(3.1) determined by the following nonlinear system of equations:$(\nabla\hat{u}_{h}, \nabla v_{h})=(f(\hat{u}_{h}), v_{h})$,$v_{h}\in S_{h}$.

(3.27)

We consider the more detailed procedure as below.
Observe that, for arbitrary $u=u_{h}\oplus u_{\perp}\in U=U_{h}\oplus U\perp$ ,$P_{h}N(u)$$=$$P_{h}u-[I-P_{h}A’(\hat{u}_{h})]_{h}-1(P_{h}u-P_{h}F(u))$$=$$[\mathrm{I} -P_{h}A’(\hat{u}_{h})]_{h}-1$$=$$[I-P_{h}A’(\hat{u}_{h})]hP-1Kh(f(u)-f’(\hat{u}_{h})u_{h})$( $P_{h}F(u)-P_{h}A’$ (\^uh)uh),

where we used the fact that $P_{h}u=u_{h}$ , and in t.he last right-hand side, we supposed that.
for simplicity. It is not necessary but, usually, we take$((B_{N})_{i})$Therefore, as in the previous paragraph, we choose the interval vector
satisfying$A’(u_{h})\wedge=f’(u_{h})\wedge$$A’(\hat{u}_{h})\approx f’(\hat{u}_{h})$$P_{h}N(U) \subset\sum_{i=1}^{M}(BN)_{i}\emptyset i$.

(3.28)and the $M$
Actually, if we define the $M\cross M$ matrix$\mathrm{b}_{N}:=((f(U)-f’(u\wedgeh)uh,\phi_{i}))$$((B_{N})_{i})$, thenis determined
dimensional interval vector
by solving the linear equation$G_{N}.--((\nabla\phi_{i}, \nabla\phi_{j})-(f’(\hat{u}_{h})u_{h}, \phi_{j}))$$G_{N}\cdot((B_{N})_{i})=\mathrm{b}_{\mathrm{N}}$.

(3.29)

39
Here,means the interval enclosure of the set {( $f(U)$$\phi_{i})\in R|u=u_{h}\oplus u_{\perp}\in U\}$ as before.
On the other hand, the error bound (rounding error)
is determined exactly same$(f(U)-f’(u_{h})\wedge uh, \phi_{i})\in IR$-$f’(\hat{u}_{h})u_{h},$$\beta_{N}$as (3.18), i.e.,$\beta_{N}=C_{0}h\sup_{\in uU}|1f(u)||_{L^{2}}$.

(3.30)

Then, we get the following computable verification condition of the same type as in
Theorem 3.1.

Theorem 3.2 For the sets (3.14), (3.15) and (3.11), let $((B_{N})_{i})$ be a solution of (3.29)
anda real number defined by (3.30). If the following conditions hold$\beta_{N}$$(B_{N})_{i}\subset u_{i}+A_{i}$$\beta_{N}$then there exists a solution$u$$\leq$$\alpha$,$i=1,$$\cdots,$ $n$,(3.31),

of $u=F(u)$ in the set$V:= \hat{u}_{h}+\sum_{i=1}^{M}(BN)_{i}\emptyset i+[\beta_{N}]$.

By using the above theorem, one can readily obtain a verification procedure based on
the Newton-like iteration, which is similarly described to that of the simple iteration
algorithm S-l.

Numerical example 2 ([80]).
We considered the following two dimensional Allen-Cahn equation which plays an important role in the mathematical biology:$-\triangle u$$=$$\lambda u(u-a)(1-u)$

in

$u$

$=$

$0$

on

$\Omega$,(3.32)$\{$$\partial\Omega$.

Here, $\Omega=(0,1)\cross(0,1)$ and the constant is taken to be $0<a<1/2$ by the reason
in actual problems. And is a positive parameter. For fixed , it is known that theequation (3.32) has two non-trivial solutions for eachwith a certain positive .
But the exact value for the criticalcorresponding to the turning point is unknown by
any theoretical approaches.
We verified several upper and lower bifurcated solutions on the approximate bifurcation
diagram by using the same finite element subspaceas in the previous paragraph. For
example, We enclosed an upper branch solution with the following data:. conditions: $a=0.\mathrm{O}1,$ $\lambda=150$ , mesh size: $h=1/80$ .
. results:, Maximum width of $(B_{N})_{i}\leq 0.025,$
error bound:$a$$\lambda$$a$$\lambda\geq\lambda^{*}$$\lambda^{*}$$\lambda^{*}$$S_{h}$$||\hat{u}_{h}||_{L^{\infty}}\approx 0.96$$.H_{0}^{1}$$\alpha\leq$0.0107.

Accuracy improvement by the residual method In this pargraph, we present
some improvement of accuracy and efficiency of verification by some residual techniquebased on an a posteriori error estimation for the higher order finite element $([84])$ .

40
Let

$u_{h}\wedge$

be an approximate solution of (3.1) satisfying (3.27). We take an element
as the solution to the following Poisson equation:$\overline{u}\in H^{2}(\Omega)\cap H_{0}^{1}(\Omega)$$-\triangle\overline{u}$$=$$f(\hat{u}_{h})$in$\Omega$

$=$

$0$

on

$\partial\Omega$,(3.33)$\{$$\overline{u}$.

We intend to find the exact solution around . Then, notice thatcoincides with . Therefore, while the explicit form ofwhich also implies that
unknown, the norm can be estimated as follows, by using (3.3) and (3.5):$u$$\hat{u}_{h}$$P_{h}\overline{u}$,is$v_{0}\equiv\overline{u}-\hat{u}_{h}\in S_{h}^{\perp}$$\overline{u}$$||v_{0}||_{H_{0}}1\leq C_{0}h||f(\hat{u}h)||_{L^{2}}$.

$v_{0}$(3.34)

Thus, in this case, we consider the candidate set $U$ of the form
$U= \hat{u}_{h}+v_{0}+\sum^{M}Wi\phi_{i}+[\alpha i=1]$,

where $W_{i}\in IR$ .
Setting $u=\overline{u}+w,$ $(3.1)$ is rewritten as the following residual form finding$-\triangle w$$=$$f(\hat{u}_{h}+v_{0}+w)-f(\hat{u}_{h})$

in

$w$

$=$

$0$

on

$\Omega$$w$,(3.35)$\{$$\partial\Omega$.

.

be a finite element subspace whose basis consists of the basis ofNow, let. We define theand the base functions having nonzero values on the boundary-projection ofby thetwo dimensional vector valued.
and setinto
Then we can easily obtain the llowing identity:$S_{h}^{*}\subset H^{1}$$\partial\Omega$$S_{h}$$L^{2}$$\mathrm{f}\mathrm{u}\mathrm{n}\mathrm{c}\mathrm{t}\mathrm{i}_{0}\mathrm{n}\overline{\nabla}\hat{u}_{h}\in S_{h}^{*}\cross S_{h}^{*}$$L^{2}(\Omega)\cross L^{2}(\Omega)$$\nabla\hat{u}_{h}\in$$\overline{\triangle}\hat{u}_{h}\equiv\nabla\cdot\overline{\nabla}\hat{u}_{h}$$S_{h^{\mathrm{X}S}h}^{*}*$$\mathrm{f}o$$(\overline{\nabla}\hat{u}_{h}, \nabla\phi)+(\overline{\triangle}\hat{u}_{h}, \phi)=0$,$\forall_{\phi}\in H_{0}^{1}(\Omega)$(3.36).

By using the above equality, we have$(\nabla v_{0}, \nabla\emptyset)=(\overline{\nabla}\hat{u}_{h}-\nabla\hat{u}h, \nabla\phi)+(\overline{\triangle}\hat{u}_{h}+f(\hat{u}_{h}), \emptyset)$

Choosing $\phi=v_{0}$ and using the Aubin-Nitsche inequlity
following estimate:$||v_{0}||_{H_{0}^{1}}$$\leq$,$\forall\phi\in H_{0}1(\Omega)$$||v_{0}||_{L^{2}}\leq C_{0}h||v0||_{H_{0}^{1}}$$||\overline{\nabla}\hat{u}_{h}-\nabla\hat{u}_{h}||_{L^{2}}+C_{0}h||\overline{\triangle}\hat{u}h+f(\hat{u}_{h})||_{L^{2}}$.

.
, we get the(3.37)

Notice that, if we use the higher order element, the right-hand side of (3.37) will converge
to with higher than $O(h)$ . Thus we can expect that the a posteriori error estimates
for the right-hand side can be actually smaller than the a priori estimates (3.34), i.e.,$O(h)$ . Therefore, it will be possible that, if we use a higher order finite element to obtain, then it will enable us to verify with high accuracy even for thethe approximation
relatively rough mesh. Moreover, this technique can also be extended for the nonconvex
and nonsmooth domain such as -shape domain in which the solution has low regularity
and we could no longer apply the a priori estimates ([83]).
$0$$\hat{u}_{h}$

$L$

41
Remark 3.2 We consider the meaning of the initial error
by $H^{-1}$ and $lei\ll.,$ $\cdot\gg be$ the duality pairing onof, observe that
$H_{0}^{1}$$v_{0}$. Let denote the dual space. Taking account of$H^{-1}\mathrm{x}H_{0}^{1}$$\triangle\hat{u}_{h}\in H^{-1}$$||\triangle\hat{u}_{h}+f(uh)\wedge||_{H}-1$$=$$\sup_{0\neq\phi\in H_{0}^{1}}\frac{\ll\triangle\hat{u}_{h}+f(\hat{u}_{h}),\phi\gg}{||\phi||_{H_{0}^{1}}}$$=$$0 \neq\phi\in H\mathrm{s}\mathrm{u}\mathrm{p}\frac{(\nabla(-\hat{u}_{h}+\overline{u}),\nabla\phi)}{||\phi||_{H_{0}^{1}}}0^{1}$$=$$||-\hat{u}_{h}+\overline{u}||H^{1}0$$=$$||v_{0}||_{H}0^{1}$.

Therefore, we can say thatstands for an element in
Riesz representation theorem from the residual functional$H_{0}^{1}(\Omega)$$v_{0}$

which is determined by the.

$\triangle\hat{u}_{h}+f(\hat{u}_{h})\in H^{-1}$

Numerical example 3 ([84]).
We verified a solution of the following Emden’s equation on the unit square in$-\triangle u$$=$$u^{2}$

in

$u$

$=$

$0$

on

$\Omega$,.

(3.38)$\{$$\partial\Omega$$R^{2}$.

on the uniform rectangular mesh withrnesh size . Therefore, we can use the constant:. We show the verification
results in Table 1. for several mesh sizees $h=1/12,$ . , 1/20. Note that, in case of thelinear element with some other residual technique, we needed at least $h=1/80$ for theverification due to the large righthandi.e.,, which confirms us
the effectiveness of the present a posteriori technique. On furher comparison with linear
a posteriori method, see [84].
We used the biquadratic finite element subspace$h$$S_{h}$$C_{0}= \frac{1}{2\pi}$$\cdot’$$\mathrm{s}\mathrm{i}\mathrm{d}\mathrm{e}\rangle$$||\hat{u}_{h}^{2}||_{L\infty}\approx 900([80])$

Table 1. Verifiation results for Emden’s equation

Some other extentions:of the present method.
(i)

We briefly mention about the improvements or extentions

error bounds This can be done by using the following constructive a priori-projection of the Poisson equation (3.4) :error estimates for the$L^{\infty}$$L^{\infty}$$H_{0}^{1}$$||\phi-P_{h}\emptyset||_{L^{\infty}(\Omega)}\leq c_{0^{\infty}}^{()_{h|}}|\psi||$.

(3.39)For example, it can be taken as $C_{0}^{(\infty)}=1.054([35])$ and 0.831([41]) for the uniform
bilinear and biquadratic element with rectangular mesh, respectively. For the triangular case, we can take, e.g., $C_{0}^{(\infty)}=1.818$ for the linear element with uniform mesh

42
error estimates was also presented([35]). The residual method using a posterioriin [41]. As a numerical example, in [41], the solution for Emden’s equation (3.38)
was verified with the following error bound in the neighborhood of the approximate
solution :$L^{\infty}$$u$$\hat{u}_{h}$$||u-\hat{u}_{h}||_{L^{\infty}}\leq 0.814$$( \frac{||u-\hat{u}_{h}||_{L^{\infty}}}{||\hat{u}_{h}||_{L^{\infty}}}\leq 0.0275)$,

where the same finite element subspace as in Numerical example 3. was adopted
with $h=1/20$ .
(ii) Verification with local uniqueness Based on the verification method describedabove, one can formulate a method to prove the existence and local uniqueness byusing the Banach fixed point theorem. We now, by using the Newton-like operator, write the residual equation (3.35) as:$T$$w=T(w)$ .

Then we consider the following set of functions$T(0)+T’(W)W$$:=$$\{v\in H_{0}^{1}(I)|v=T(0)+T’(\tilde{w})w, \tilde{w}, w\in W\}$,

where $T’$ stands for the Fr\’echet derivative of . By appropriate bounding of thisset, one can prove that is a contraction map on the candidate set $W=W_{h}\oplus W_{\perp}$
satisfying $T(W)\subset W$ . Therefore, by Banach’s fixed point theorem, we obtain the
inclusion result that there exists a unique solution $w\in W$ with $w=T(w)$ . For
details and examples, see [85].
(iii) Parameter dependent equations For the parameter dependent elliptic problems
such as (3.32), it is important to enclose the solution curve itself or to verify theexistence of some singular points, e.g., turning point or bifurcation point. It is also
possible to apply our verification method to these problems by using some additional
techniques such as the solution curve enclosing for regular solutions [43] and the bordering technique for turning p\‘Oints [76], [21]. The verification for simple bifurcation
points would also be possible by applying the similar method to that in $[77],[66]$ .
(iv) Navier-Stokes equations In [45], an a posteriori and a constructive a priori error
estimates for finite element solutions of the Stokes equations were presented. By
using these results a prototype verification has been derived in [82] for the solutionsof Navier-Stokes problems with small Reynolds number and small solutions. Recently, the present method was applied to verify the bifurcated periodic solutions of
the Rayleigh-B\’enard problem for the heat convection in a two dimensional domain[47]. By these validated computational results, several properties which were not yet
proved by the theoretical approaches were numerically verified.
(v) Variational inequalities Our method can also be applied to the enclosing solutionsof the variational inequalities with nonlinear right-hand side of the form [71]:
$T$

$T$

Find $u\in K$ such that
(3.40)$\{$$(\nabla u, \nabla(v-u))\geq(f(u), v-u),$ $\forall v\in K$,

43which appears in arguments on the obstacle problems $([10])$ . Here, $K=\{v\in$. Furthermore, in [46] and [72], somewhat different kind of variational inequalities are treated.
(vi) Estimation of the optimal constant When we apply our verification method tothe boundary value problems in non-rectangular polygonal domains, we need thein (3.3) as small as possible. Therefore, it is important to estimate theconstant. The problem of finding such a constant is reduced
optial value of the constantto the following an eigenvalue-like problem of the Laplacian:$H_{0}^{1}(\Omega)|v$$\geq 0\}$

$C_{0}$$C_{0}$$-\triangle u$$\frac{\partial u}{\partial n}$$\{$$\int_{\Gamma_{3}}ud_{S}$$=$$\lambda u+\lambda\psi+\triangle\psi$

$=$

$0$

$=$

$0$

$on\partial\Omega$in$\Omega$,,(3.41),

whose vertices are points $(0,0),$ $(0,1)$ andis the standard triangle inwhere$(1, 0)$ , andon the -axis. When we denoteis the edge from (0.0) to $(0,1)$ in
the minimal value in the set of the all solutions of (3.41) by , the optimal constant$R^{2}$$\Omega$$\partial\Omega$$\Gamma_{3}$$\mathrm{y}$$\lambda$$C_{0}$$\lambda_{0}$looking for satisfies$C_{0} \leq\frac{1}{\sqrt{\lambda_{0}}}$. By applying our numerical verification method forthe solution of (3.41), we obtained the estimate $C_{0}\leq 0.494[48]$ which is considered
as a really significant improvement of the existing best constant $C_{0}\leq 0.81$ . On the, we can say this would be ‘nearly’
other hand, since it is known that 0.467
optimal estimate of .
$\leq C_{0}$$C_{0}$

3.1.2

Plum’s method

Plum presented a verification principle, which belongs to the category of analytic
method, for the solutions of nonlinear elliptic problems in [61], [62] incorporating withhis results on the enclosing eigenvalues for elliptic operators.
Let assume that the nonlinear function in (3.1) is sufficiently smooth in with respect
to each variable.
smoothness is obAlso assume that an approximate solution of (3.1) with-element when the finite element approximation
tained. This means that we use theis adopted for (3.1). And the residual is bounded by a sufficiently small as follows:$f$$H^{2}(\Omega)$$\omega$$C^{1}$$\delta$.

(3.42)$||-\triangle\omega+f(\omega)||_{L^{2}}\leq\delta$

Also, suppose that there exists a constant $K$ satisfying$||u||_{L}\infty\leq K||Lu||_{L^{2}}$where$Lu \equiv-\triangle u+\frac{\partial f}{\partial u}(\omega)u$,$\forall u\in H^{2}$,(3.43). The above constant $K$ is determined by the arguments

described in later.
Next, assume that the following non-decreasing function$g$: $[0, \infty)arrow[0, \infty)$ can bechosen as$|f(\omega(X)+y)-f(\omega(X))-J(x)y|\leq g(|y|)(_{X\in\Omega,y\in}R)$ ,(3.44)

44
where $g(t)=o(t)$ as$tarrow \mathrm{O}$and$J(x) \equiv\frac{\partial f}{\partial u}(\omega(x))$.

Moreover, let us suppose that there exist positive constants
for any $u\in H^{2}$ ,$||u||_{L^{2}}\leq K_{0}||Lu||L^{2}$,$||\nabla u||_{L^{2}}\leq K_{1}||Lu||_{L^{2}}$,$K_{i},$$i=0,1,2$ such that,$|u|_{H^{2}}\leq K_{2}||Lu||_{L^{2}}$.

(3.45)

Here,can be determined by the later consideration. Then the following theorem
provides the verification condition of the exact solution to (3.1) in a neighborhood of$K_{i}$$u$$\omega$.

Theorem 3.3 For$\alpha\geq 0$, if the inequality$\delta\leq\frac{\alpha}{K}-\sqrt{|\Omega|}\cdot g(\alpha)$

holds, then there exists a solution
the measure of .

$u$,

(3.46)of (3.1) satisfying$||u-\omega||_{L^{\infty}}\leq\alpha$. Here,$|\Omega|$means$\Omega$defined by the residual simpliThe condition (3.46) implies that the map on
is retractive on the set $D\equiv\{u\in$fied Newton operator for the equation (3.1) at. Therefore, the verification is based on the Schauder fixed point. In order to verify the solution by this method, we need a fine approxtheorem inas well as the estimates of the constants appeared in the above.
imationParticularly, the constant $K$ in (3.43), which provides the inverse norm for the linearlized
operator of the original equation, plays the most essential role. This constant is estimatedas follows (see [62] for details):
First, we need an estimate of the following positive constant a.
$L^{\infty}(\Omega)$$\omega$$L^{\infty}(\Omega)|||u|\}_{L^{\infty}}\leq\alpha\}$$L^{\infty}(\Omega)$$\omega\in H^{2}\cap H_{0}^{1}$$\sigma\leq\min${$|\lambda||\lambda$

is the eigenvalue of operator

$L$

on

$H^{2}\cap H_{0}^{1}$}.

In order to get a lower bound of the eigenvalues for , he uses the homotopy method
starting with some simple operator, e.g., $L=\triangle$ , whose lower bound of eigenvalue isalready known ([58]).
in (3.45) are determined as follows.
Next, using this , the constants.
We only describe here for the case that is independent ofsuch thatFirst, it can be taken as $K_{0}=1/\sigma$ . And defining the$L$$K_{i}$$\sigma$$\nabla u$$f$$\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{t}_{\mathrm{S}}\underline{c}\mathrm{a}\mathrm{n}\mathrm{d}_{\overline{C}}$$\underline{c}\leq J(x)\leq\overline{c}$$K_{1}$,$\forall_{X\in\overline{\Omega}}$,(3.47)

is given by$[K_{0}(1-\underline{C}K0)]^{1/2}$$K_{1}=\{$$\frac{1}{2\sqrt{\underline{c}}}$Also,$K_{2}$, if$\underline{c}K_{0}\leq 1/2$

(3.48)

otherwise.

is decided as:$K_{2}=1+K0^{\cdot} \max\{-\underline{C}, \frac{1}{2}(\overline{c}-\underline{C})\}$

,

.

45
Thus, if the constants$C_{j},$$j=0,1,2$ in the Sobolev inequality:(3.49)$||u||L\infty\leq C_{0}||u||_{L^{2}}+C_{1}||\nabla u||_{L^{2}}+C_{2}|u|_{H^{2}}$can be numerically estimated, then the desired constant $K$ in (3.43) is obtained by$K\equiv C_{0}K_{0}+C_{1}K_{1}+C_{2}K_{2}$ . On the other hand,in (3.49) is computed as follows [62].
$C_{j}$,$j=0,1,2$ ,$C_{j}= \frac{\gamma_{j}}{|\Omega|}\cdot[\mathrm{m}\mathrm{a}_{\frac{\mathrm{x}}{\Omega}}x\mathrm{o}\in\int_{\Omega}|x-x0|2\mathcal{U}dX]^{1}/2$(3.50), 1.1548, 0.22361 for $n=2$ , andwhere,(1.0708, 1.6549, 0.41413) for $n=3$ .
As the numerical examples, in [61], he verified the solutions, for several parameter , on the unit square inwith homogeneous Dirichlet condition. The
approximate solution is constituted as the piecewise biquintic-spline on the 8
finite element mesh.
In the meantime, he extended the method to the equation having turning points andbifurcation points in [65] and [66], respectively, as well as the problem on the nonconvex
nonsmooth domain in [64].
$(\gamma 0, \gamma 1, \gamma 2)=(1$$(\gamma 0, \gamma 1, \gamma 2)=$$)$$\mathrm{o}\mathrm{f}-\triangle u=\lambda e^{u}$$R^{2}$$\lambda$$C^{1}$$\omega$$\cross 8$

Heywood’s method

3.1.3

In [12], Heywood et al. presented a numerical verification method for the spatiallyperiodic solutions of the steady Navier-Stokes equations by usihg the spectral techniques.
They considered the following problem on the fundamental domain $\Omega=(0,1)\mathrm{x}(0,1)$ ofperiodicity :,$x\in\Omega$$-I\text{ノ}\triangle u+(u\cdot\nabla)u+\nabla p-f$

$=$

$0$

$\nabla\cdot u$

$=$

$0$

$u(x+K)$$=$$u(x)$$\int_{\Omega}u(x)dx$

$=$

$0$

$\{$

,

$x\in\Omega$,,$\forall x\in R^{2},$$K\in Z^{2}$,(3.51),

where is a prescribed force satisfying the same periodic condition as , and l ノ is the
kinematic viscosity.
Under some appropriate setting of the spatially periodic function space, i.e.,
and divergence free}, the original equation (3.51) can be written as :{such that
find
$f$

$u$

$V\approx$$H^{2}(\Omega)\cap$$\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{i}\mathrm{o}\mathrm{d}\mathrm{i}\mathrm{c}$$\exists u\in V$$F(u)\equiv-I\text{ノ}\triangle u+(u\cdot\nabla)u-f=0$.

(3.52)

Now let $M$ be a constant satisfying$||(u\cdot\nabla)v||_{L^{2}}\leq M||u||_{V}||v||_{V}$,$\forall u,$.

$v\in V$(3.53)For some approximate solution \^u of (3.52), denoting the Fr\’echet derivative of $F$ at \^u by$L\equiv F’$(\^u) : $Varrow L^{2}$ , the following theorem is obtained.

46
Theorem 3.4 If$L$is invertible with a bound $||L^{-1}||$ and$4M||L^{-}1||^{2}||F(\hat{u})||_{L}2\leq 1$

holds, then (3.52) has a locally unique solution$u$satisfying$||u-\hat{u}||_{V}\leq 2||L^{-1}|||F(\hat{u})||_{L^{2}}$.

This theorem can be proved by the fact that the simplified Newton operator at \^u definesa contraction mapping on a neighborhood of it under the above condition.
The main task of the verification procedure by the application of Theorem 3.4 is concerned with estimating the norm $||L^{-1}||$ . This is done by a reduction in two steps. The
first step is to approximate with a simpler, but still infinite dimensional, linear operahaving the same finite part as , where $N$ is a positive interger, and converging
tor
is introduced as the spectral
to as $Narrow\infty$ . Next, a finite dimensional operator
to some finite dimensional
Galerkin approximation of , which is the restriction ofis invertible andis bounded through a reduction
space, and it is shown that
. Thus, the general perturbation theorem yields theto the property with respect to$||L^{-1}||$by using the error estimate $||L-L_{N}||$ .
inveritibility of and the bound
They applied the method to the case that $l\text{ノ}=0.001$ and $f(x)\equiv f(x_{1}, x_{2})$ is a simple
and got some verified results by constructing
numerically whichsine function inapproximationtoaof restictedfinite dimensional space withis the spectral Galerkindimension 2,820. Then, the quantities in Theorem 4 were estimated as follows:
$L$

$L$

$L_{N}$$L$$\overline{L}_{N}$$L$$L_{N}$$||L_{N}^{-1}||$$L_{N}$$\overline{L}_{N}$$L$$\overline{L}_{120}$$x_{2}$$L$$M=15.42493$ ,$||L^{-1}||<390.16$ ,$||F(\hat{u})||_{L^{2}}=3.421\cross 10^{-8}$.

and thus they obtained$M||L^{-1}||^{2}||F(\hat{u})||_{L^{2}}\leq 0.32132$,

which implies by the theorem that there exists a locally unique solution
the error bound$||u-\hat{u}||_{V}\leq 2.7\cross 10^{-5}$$u$of (3.52) with.

In these numerical computations, they neglected the round-off error of floating point
arithmetic with double precision as well as used commercially available software concernedwith linear algebra without verification.

3.2

Evolution equations

The study for the numerical verification method for evolution problems has been still
made less progress than for the elliptic case.
We consider the following parabolic problems$=$$f(x, t, u)$ ,$u(x, t)$

$=$

$0$

$u(x, 0)$

$=$

$0$

$\frac{\partial u}{\partial t}-\triangle u$$(x, t)\in\Omega\cross J$

$\{$

,
,

,

$(x, t)\in\partial\Omega\cross J$$x\in\Omega$,,(3.54)

47
is a convex domain inor
and $J=(0, T)(T>0)$ . If the timeisfixed, then the equation (3.54) can be treated as a kind of stationary problem. Thus,
for example, the the author’s method in the previous subsection can also be applied, inprinciple, to the verification of solutions of this problem. In such applications, the simple
linear problem which corresponds to the Poisson equation in the elliptic case is as follows:where$R^{1}$$\Omega$$R^{2}$$\frac{\partial\phi}{\partial t}-\triangle\emptyset$

$=$

$g$

$\phi(x, t)$

$=$

$0$

$\phi(x, 0)$

$=$

$0$

$T$

$(x, t)\in\Omega\cross J$,(3.55)

$\{$

,
,

,

$(x, t)\in\partial\Omega\cross J$$x\in\Omega$,where is the prescribed function. Therefore, if one obtain a fixed point formulationin the appropriate function space and the constructive a priori error estimates for thefinite element solution of (3.55), then the arguments in the elliptic case can be applied
to the verification for the problem (3.54). In [32] and [38], the verification examples were
presented based on the simple iteration method for one and two space dimensional cases,
respectively.
On the other hand, noting that it is read\’ily possible to estimate the norm of the inverseof a linearized operator for (3.54), Plum’s method can also be applied to the same problem
without any complicated work concerning the eigenvalue enclosing for the operator. In
line with this direction, Minamoto [19], [22] formulated and presented some verification
examples based on the residual Newton type operator. Also for the case of the following
hyperbolic equations, in [20], [23], by using the similar arguments to parabolic case,several verification results were obtained.
$g$

$=$

$f(x, t, u)$ ,$u(x, t)$

$=$

$0$

$u(x, 0)$

$=$

$0$

$=$

$0$

$\frac{\partial^{2}u}{\partial t^{2}}-\triangle u$

$\{$

,
,

$(x, t)\in\Omega\cross J$,,$(x, t)\in\partial\Omega\cross J$$x\in\Omega$,$x\in\Omega$.

(3.56)$\partial u$$\overline{\partial t}(x, 0)$,

Since, in these verification procedures by analytic methods, the residual Newton method
are utilized, the accuracy of the appriximate solution essentially affects the success ofverification. Usually, due to the difficulty to improve the accuracy for time direction, it is
not so easy to compute an approximation with sufficiently small residue for the practicalproblems.
Recently, in [14], Kawanago proposed a method to verify a time periodic solution forsome bifurcation problem on the following semilinear dissipative wave equation:$u_{tt}-u_{xx}+u_{t}+u^{3}$$=$$\lambda\sin x\cos t$$u(\mathrm{O}, t)=u(\pi, t)$

$=$

$0$

,

$(x, t)\in(0, \pi)\cross(0, \infty)$,(3.57)$\{$,$t\in(0, \infty)$,

where is a positive parameter. He numerically proved that the above equation has aaround an approximate solution
symmetry-breakung pitchfork bifurcation pointof
Newtontype by using spectral basis
analytic
method
.
used
a
kind
of
He
pair
for both of space and time.
$\lambda$$(\lambda_{0}, u_{0})$$(\tilde{\lambda}_{0},\tilde{u}_{0})$48
Eigenvalue problems

4

Consider the following self-adjoint elliptic eigenvalue problem:,$Au\equiv-\Delta u+qu$$=$$\lambda u$

$u$

$=$

$0$

$R^{2}$

and$q\in L^{\infty}(\Omega)$$x\in\Omega$

,
(4.1)

$\{$

is a bounded convex domain in
where(4.1) as :$\Omega$find$\exists(u, \lambda)\in H_{0}^{1}(\Omega)\cross R$,$x\in\partial\Omega$,. We normalize the problem

satisfying$-\triangle u+(q-\lambda)u$

$=$

$0$

,

$x\in\Omega$

,
(4.2)

$\{$

$=$

$\int_{\Omega}u^{2}$

1.

As in the usual elliptic problems, it is readily seen that the problem (4.2) can be rewritten, and that, in order to enclose
as a fixed point equation of a compact map on, we can also apply
around an approximate pairthe eigenpair
the author’s verification method described in the previous section (see [44], [26], [27]for details). Namely, the exact eigenvalue and the corresponding eigenfunction areof the form $U=$
enclosed in an interval A $\in IR$ and in a candidate set$\hat{u}_{h}+v_{0}+U_{h}+[\alpha]$ , respectively.
Here, as in the previous section, $U_{h}\in S_{h,I},$. Furthermore, in the present case, we can assert the uniquenessandof the eigenvalue in A by the almost same verification condition as in Theorem 3.1 or
3.2([27]). Note that this method gives the enclosure of eigenvalues as well as verifies thecorresponding eigenfunctions in contrast to other methods, described later, which present
only eigenvalue enclosing.
We now briefly remark on the method to enclose the eigenvalues in order of magnitude.
In such a case, an eigenvalue excluding procedure plays an essentially important role. Thiscan be done as follows.
,
We consider an sufficiently narrow interval , and set, for a$H_{0}^{1}(\Omega)\cross R$$(\hat{u}_{h}, \lambda_{h})\in S_{h}\cross R$$(u, \lambda)$$\lambda$$u$$U\subset H_{0}^{1}(\Omega)$$v_{0}\in$$S_{h}^{\perp}$$[\alpha]\subseteq S_{h}^{\perp}$$\lambda\in\Lambda$$\Lambda$$L(\lambda)\equiv-\triangle u+(q-\lambda)u$

Then, since$u=0$ :

$L(\lambda)$.

is a linear elliptic operator, the following equation has a trivial solution$L(\lambda)u$

$=$

$0$

,

$u$

$=$

$0$

,

$x\in\Omega$

,
(4.3)

$\{$

$x\in\partial\Omega$.

Therefore, if we validate the uniqueness of the solution in (4.3), then it implies that is
not an eigenvalue of (4.1). The uniqueness property can be proved, taking into accountis linear, by the method analogous to that in the previous section.
that the operator
Thus, by using some interval computing techniques, it is easily to verify that there is no
eigenvalue of (4.3) in A. Thus the eigenvalue excluding process advances from the one to
the next, backward or forward to the adjacent intervals.
We now note that, by some eigenvalue shift, we can easily present the lower bound of
the spectrum of (4.1). Therefore, by appropriately combining this excluding procedure$\lambda$$L(\lambda)$

49
with the enclosing technique described before, we make the eigenvalue ordering as far as
each eigenvalue is geometrically simple $([28])$ .

A numerical example:
Particularly, the eigenvalue with smallest absolute value(ESAV), which is important toverify the solution of the corresponding nonlinear elliptic equations, was verified for thefollowing problem $([44])$ :$-\triangle u+l\text{ノ}(3v_{h}^{2}-2\wedge(a+1)v_{h}\wedge+a)u$$=$$\lambda u$

$u$

$=$

$0$

$x\in\Omega$,$x\in\partial\Omega$

.

,

(4.4)

$\{$

,

This is the eigenvalue problem for the linearlized operator of the Allen-Cahn equation(3.32) at . Hereis an approximate lower branch solution of the original problem inthe finite element subspaceof biquadratic polynomials as in the numerical example 2
in the previous section with $l\text{ }=150,$
and mesh size $h=1/20$ . The ESAV of(4.4) were enclosed using the present method in the interval:$v_{h}\wedge$$v_{h}\wedge$$S_{h}$

ノ$a=0.\mathrm{O}1$$\Lambda=[-16.67017, - 16.55062]$

and the corresponding eigenfunction was enclosed in the set : $U=\hat{u}_{h}+v_{0}+U_{h}+[\alpha]$ ,$||v_{0}||_{L^{2}}\leq 0.0001372$ , [maximum width of the coefficient intervals
where$\alpha\leq0.00105$ .
$\leq
0.00321$in ]and$||\hat{u}_{h}||_{L^{\infty}}\approx 2.308,$$U_{h}$

Remark 4.1 It is seen that the above method can also be applied to the non-selfadjoint
eigenvalue problems of the form:$-\triangle u+p\cdot\nabla u+qu$$=$$\lambda u$

$u$

$=$

$0$

,

$x\in\Omega$

,
(4.5)

$\{$

,

$x\in\partial\Omega$.

Now, we will mention about other methods. The following proposition is well known
as the Weinstein bounds for the eigenvalues of the form $Au=\lambda u$ :Proposition 4.1 Let$D(A)$ ,$(\tilde{u},\tilde{\lambda})$be an approximate eigenpair for $Au=\lambda u$ such that

where $D(A)$ is the domain of$A$, and let

least one exact eigenvalue in the interval$\epsilon:=\frac{||A\tilde{u}-\tilde{\lambda}\tilde{u}||}{||\tilde{u}||}$$[\tilde{\lambda}-\epsilon,\tilde{\lambda}+\epsilon]$$\tilde{u}\in$. Then there exists at.

This proposition was extended to more practical version known as Kato’s bound $([67])$ .
On the other hand, Lehman-Goerisch method [9], [4] is well-known, and is based on
the Rayleigh-Ritz method which gives only upper bounds to the eigenvalues. We now
describe the outline of this method. Let denote the N-th eigenvalue byordered frombe the largest eigenvalue of the discretized problem for (4.1)
the smallest one, and let
as a set of trial functions. We assume
by using the approximate eigenfunctions
that some $\rho\in R$ is known satisfying$\lambda_{N}$$\Lambda_{N}$$\{\tilde{u}_{i}\}_{i=1}^{N}$$\Lambda_{N}<\rho<\lambda_{N+1}$

.

(4.6)

50
Moreover, define the following matrices:$(A_{1})_{ij}:=(A\overline{u}i, A\tilde{u}j)$,$(A_{2})_{ij}:=((A-\rho)\tilde{u}_{i},\tilde{u}j)$,$(A_{3})_{ij}:=((A-\rho)\overline{u}i, (A-\rho)\tilde{u}_{j})$

(4.7)
,

inner product on . Then, by the asstands for theandwhere $j=1,$
sumption (4.6), the generalized eigenvalue problem: $A_{2}x=\mu A3X$ has negative eigenvalues. Then we have$\cdots,\dot{N}$$i,$$L^{2}$$(\cdot, \cdot)$$\Omega$$\{\mu_{i}\}_{i=1}^{N}$

Theorem 4.1 Assuming the above conditions, it holds thatfor$\lambda_{N+1-i}\geq\rho+\frac{1}{\mu_{i}}$$i=1,$$\cdots,$$N$.

By using this theorem one can determine the lower bounds of eigenvalues provided thatthe spectral parameter is a priori known by some other consideration.

An example:
In [6], Behnke applied some extended method of the above to the following eigenvalue
problem of forth order related to the vibrations of a clamped plate:$x\in\Omega$$\frac{\partial^{4}}{\partial x^{4}}u+P\frac{\partial^{4}}{\partial x^{2}\partial y^{2}}u+Q\frac{\partial^{4}}{\partial y^{4}}u=\lambda u$

$u=0$

and

$x\in\partial\Omega$$\frac{\partial u}{\partial n}=0$,,

are positive constants and $\Omega=\cross(-\frac{b}{2}, \frac{b}{2})$ . He considered the eigenwhere $P,$
values as functions of $s=a/b$ and numerically proved a curve veering phenomena $(Cf$.
[5] .
Now, the Lehman-Goerisch method a priori needs a spectral parameter in (4.6). It
is, in general, not necessarily easy to decide such a parameter. Plum [58], [60], [67]introduced a homotopy method which overcomes this difficlty by connecting the given
problem (4.1) with a simple problem whose spectrum is already known.
For example, the problem (4.1) is connected to the simple eigenvalue problem for the
Laplacian by the following homotopy, for $t\in[0,1]$ ,
$Q$

$)$

$\rho$$A_{t}u\equiv-\triangle u+tqu$$=$$\lambda u$,$x\in\Omega$

,
(4.8)

$\{$

$u$

$=$

$0$

,

$x\in\partial\Omega$.

the homotopy is set such that
When we denote the n-th eigenvalue of (4.8) byis monotonically non-decreasing with respect to . Therefore, basically both spectra
have invariant structure. Starting at some $t=t_{1}\in(0,1]$ , by repeated applications of
the Lehmann-Goerisch method or other bounding methods with additional procedures,are enclosed, then increase , e.g., $t=t_{2}\in(t_{1},1]$ , little by little
ordered eigenvalues
up to the final stage: $t=1$ .
$\lambda_{n}^{t}$$n$$\lambda_{n^{1}}^{t}$$t$$\lambda_{n}^{t}$

51
Examples:In [60], he enclosed several eigenvalues of the following problems on the rectangulardomain $\Omega=(0, \pi)\cross(0, \pi)$ with Dirichlet boundary condition:

1. $- \triangle u+10\cos^{2}[\frac{1}{6}(2x+y)]u=\lambda u$ .
$2.- \triangle u=\frac{\lambda u}{1+\frac{1}{\pi^{2}}(_{X^{2}}1^{+2_{X_{2}}})2}$

5.

Conclusions

We have surveyed numerical verification methods for differential equations, especiallyaround PDEs and the author’s works. But the period of this research is shorter than thehistory of the numerical methods for differential equations by computer and we can say it
is still in the stage of case studies. Indeed, recently, this kind of studies have been referredlittle by little for practical applications in PDEs but there are many open problems to beresolve. Therefore, we can make no prediction that these approaches will grow into really
useful methods for various kinds of equations in mathematical analysis. Also, since theprogram description of the verification algorithm is very complicated in general, there
is another problem like software tecnology associated with assurance for the correctnessof th everification prograqm itself. Actually, some of the mathematician would not give
credit the computer assisted proof in analysis as correct as they believe the theoretical
proof, which might cause a kind of seriously emotional problemm in the methodology ofmathematical sciences. And there is another difficulty from the huge scale of numericalcomputations which often exceed the capacity of the concurrent computing facilities.
However, in the 21st century, the computing environment would make more and more
rapid progress, which should be beyond conception in the present state. The author
believes that the above difficulties should be overcome by this evolution and that thecomputer would greatly contribute to the theory of analysis in mathematics and create
a new research area which should be called computer aided analysis. On the other hand,
numerical methods with guaranteed accuracy for differential equations would highly improve the reliability in the numerical simulation of the complicated phenomena in sienceand technology.

$\#\vee \mathit{4}’\#\mathrm{X}\mathrm{f}\mathrm{f}\mathrm{i}^{\backslash }$[1] Adams, R.A., Sobolev spaces, Academic Press, New York, 1975.
[2] Adams, E. &Kulisch, U., Scientific computing with automatic result verification,
Academic Press, San Diego, 1993.
[3] Alefeld, G. &Herzberger, J., Introduction to interval computations, Academic Press,
New York, 1983.

[4] Behnke, H., The determination of guaranteed bounds to eigenvalues with the useof variational methods , Computer arithmetic and self-validating numerical meth(Ullrich, C. ed.), Academic Press, San Diego, 1990, 155-170.
$1\mathrm{I}$$\mathrm{o}\mathrm{d}\mathrm{s}$52[5] Behnke, H., A numerically rigorous proof of curve veering in an eigenvalue problem
for differential equations, Zeitschrift f\"ur Analysis und ihre Anwendungen 15 (1996),

181-200.

[6] Behnke, H., On the eigenvalue problem of the rectangular clamped plate, private

communication.

[7] Cesari, L., Functional analysis and periodic solutions for one dimensional problems,
Contrib. Differential Equations 1 (1963), 149-187.

[8] Ciarlet, P.G., The finite element method for elliptic problems, North-Holland, Amsterdam, 1978.
[9] Goerisch, F. &He, Z., The determination of guaranteed bounds to eigenvalues withthe use of variational methods I, Computer arithmetic and self-validating numericalmethods(Ullrich, C. ed.), Academic Press, San Diego, 1990, 137-153.
[10] Glowinski, R., Numerical Methods for Nonlinear Variational Problems, Springer,
New York, 1984.
[11] Grisvard, P., Elliptic problems in nonsmooth domain, Pitman, Boston, 1985.
[12] Heywood, J.G., Nagata, W. &Xie, W., A numerically based existence theorem for
the Navier-Stokes equations, Journal of Mathematical Fluid Mechanics 1 (1999), 5-23.

. &Miranker, W.L., Self-validating numerics for function space prob[13] Kaucher,lems, Academic Press, New York, 1984.
$\mathrm{W}.\mathrm{E}$[14] Kawanago, T., Generalization of bifurcation theorems aiming for the study of nonlinear forced vibrations described by PDEs, preprint.
[15] Kedem, G., A posteriori error bounds for two-point boundary value problems, SIAM
J. Numer. Anal., 18 (1981), 431-448.
[16] Lehmann, R., Computable error bounds in the finite-element method, IMA J. Numerical Analysis 6 (1986), 265-271.

[17] Lohner, R., Enclosing the solutions of ordinary initial and boundary value problems,. Teubner Stuttgart, (1987), 255-286.
Computerarithmetic ( . Kaucher, E. et al.),$\mathrm{e}\mathrm{d}\mathrm{s}$$\mathrm{B}.\mathrm{G}$-error bounds for the ap. &Tapia, R.A., Computable a priori[18]
J. Numer. Anal. 12
SIAM
problems,
valueboundary
two-point
proximate solution of(1975), 919-937.
$\mathrm{M}\mathrm{c}\mathrm{C}\mathrm{a}\mathrm{r}\mathrm{t}\mathrm{h}\mathrm{y},$$\mathrm{M}.\mathrm{A}$$L^{\infty}$[19] Minamoto, T. &Nakao, M.T., Numerical verifications of solutions for nonlinear
parabolic equations in one-space dimensional case, Reliable Computing 3, (1997), 137-

147.

[20] Minamoto, T., Numerical verifications of solutions for nonlinear hyperbolic equations
in one-space dimensional case, Applied Mathematics Letters 10 (1997), 91-96.
[21] Minamoto, T., Yamamoto, N. &Nakao, M.T., Numerical verification method forsolutions of the perturbed Gelfand equation, to appear in Methods and Applicationsof AnalySiS.
[22] Minamoto, T., Numerical existence and uniqueness proof for solutions of semilinear
parabolic equations, to appear in Applied Mathematics Letters.

[23] Minamoto, T., Numerical existence and uniqueness proof for solutions of nonlinear
hyperbolic equations, to appear in Journal of Computational and Applied Mathematics.
, 1966.
[24] Moore, R., Interval analysis, Prentice-Hall, Englewood Cliffs,$\mathrm{N}\mathrm{J}$53[25] Mishaikow, K. &Mrozek, M., Chaos in the Lorenz equations: A computer assisted
proof. Part II: Details, Mathematics of Computation 67 (1998), 1023-1046.

[26] Nagatou, K., Yamamoto, N. &Nakao, M.T., An approach to the numerical verification of solutions for nonlinear elliptic problems with local uniqueness, Numerical
Functional Analysis and Optimization 20 (1999), 543-565.

[27] Nagatou, K., A numerical method to verify the elliptic eigenvalue problems including
a uniqueness property, Computing 63 (1999), 109-130.
[28] Nagatou, K. &Nakao, M.T., An enclosure method of eigenvalues for the elliptic
operator linearlized at an exact solution of nonlinear problems, to appear in Linear
Algebra and its Applications.
[29] Nakao, M.T., A numerical approach to the proof of existence of solutions for elliptic
problems, Japan Journal of Applied Mathematics 5, (1988) 313- 332.
[30] –, A computational verification method of existence of solutions for nonlinear elliptic equations, Lecture Notes in Num. Appl. Anal., 10, (1989) 101-120. In proc. Recent
Topics in Nonlinear PDE 4, Kyoto, 1988, $North-Ho\iota\iota and/Kinokuniya$ , 1989.
[31] –, A numerical approach to the proof of existence of solutions for elliptic problems
I, Japan Journal of Applied Mathematics 7 (1990), 477-488.

[32] Nakao, M.T., Solving nonlinear parabolic problems with result verification: Part I,
in Journal of Computational and Applied Mathematics 38 (1991) 323-334.
[33] –, A numerical verification method for the existence of weak solutions for nonlinear
boundary value problems, Journal of Mathematical Analysis and Applications 164(1992), 489-507.
[34] –, Solving noniinear elliptic problems with result verification using an $H^{-1}$ residual
iteration, Computing Suppl. 9 (1993), 161-173.

[35] –, Computableerror estimates in the finite element method with application
to nonlinear elliptic problems, Contributions in Numerical Mathematics.
ed.), World Scientific, Singapore, 1993, 309-319.
$L^{\infty}$$(\mathrm{A}\mathrm{g}\mathrm{a}\mathrm{r}\mathrm{w}\mathrm{a}\mathrm{l},$$\mathrm{R}.\mathrm{P}$. &Yamamoto, N., Numerical verifications of solutions for elliptic equa[36] Nakao,
tions with strong nonlinearity, Numerical Functional Analysis and Optimization 12(1991), 535-543.
$\mathrm{M}.\mathrm{T}$[37] –, Numerical verifications of solutions for nonlinear hyperbolic equations, Interval
Computations 4 (1994), 64-77.

. &Watanabe, Y., On computational proofs of the existence of solutions[38] Nakao,
to nonlinear parabolic problems, Journal of Computational and Applied Mathematics
50 (1994), 401-410.
$\mathrm{M}.\mathrm{T}$. &Yamamoto, N., A simplified method of numerical verification for non[39] Nakao,
linear elliptic equations, in the Proceedings of International Symposium on Nonlinear
Theory and its Applications (NOLTA’95), Las Vegas, USA, (1995), 263-266.
$\mathrm{M}.\mathrm{T}$[40] Nakao, M.T., Yamamoto, N. &Watanabe, Y., Guaranteed error bounds for finite element solutions of the Stokes problem, in Scientific Computing and Validated Numerics), Akademie Verlag, Berlin (1996), 258-264.
(G. Alefeld et al.
$\mathrm{e}\mathrm{d}\mathrm{s}.$. &Yamamoto, N., Numerical verification of solutions for nonlinear[41] Nakao,
residual method, Journal of Mathematical Analysis andelliptic problems using
Applicatons 217, (1998), 246-262.
$\mathrm{M}.\mathrm{T}$$L^{\infty}$[42] Nakao, M.T., Yamamoto, N. &Kimura, S., On best constant in the optimal error-projection into piecewise polynomial spaces, Journal of Approxiestimates for themation Theory 93, (1998), 491-500.
$H_{0}^{1}$

54[43] Nakao, M.T., Yamamoto, N. &Nishimura, Y., Numerical verification of the solution
curve for some parametrized nonlinear elliptic problem, in Proc. Third China-Japan
Seminar on Numerical Matehmatics, Aug. 26-30, Dalian, China, 1996 ( . Shi, Z.-C.
&Mori, M.), Science Press, Beijing , (1998), 238-245.
$\mathrm{e}\mathrm{d}\mathrm{s}$[44] Nakao, M.T., Yamamoto, N. &Nagatou, K., Numerical verifications of eigenvalues ofsecond-order elliptic operators, Japan Journal of Industrial and Applied Mathematics
16 (1999), 307-320.
[45] Nakao, M.T., Yamamoto, N. &Watanabe, Y., A posteriori and constructive a priori
error bounds for finite element solutions of Stokes equations, Journal of Computationaland Applied Mathematics 91 (1998), 137-158.
. &Ryoo, C.S., Numerical verification of solutions for elasto[46] Nakao, M.T., Lee,plastic torsion problems, Computers&Mathematics with Applications 39 (2000), 195$\mathrm{S}.\mathrm{H}$204.

[47] Nakao, M.T., Nobito Yamamoto, Yoshitaka Watanabe and and Takaaki Nishida,
A numerical verification of bifurcated solutions for the heat convection problems, in
preparation.

. &Yamamoto, N., Computer assisted enclosure of the optimal constant[48] Nakao,in the error estimation of a finite element projection, in preparation.
[49] Natterer, F., Berechenbare Fehlerschranken f\"ur die Finite Elemente Methode, ISNM
vol. 28, Birkh\"auser Verlag, Basel (1975), 109-121.
$\mathrm{M}.\mathrm{T}$[50] Neumaier, A., Interval methods for systems of equations, Cambridge Univ. Press,
Cambridge, 1990.
[51] Nickel, K., The construction of a priori bounds for the solution of a two point
boundary value problem with finite elements I, Computing 23 (1979), 247-265.
[52] Nishida, T., Teramoto, Y. &Yoshihara, H., Bifurcation problems for equations of
fluid dynamics and computer assisted proof, Lecture Notes in Num. Appl. Anal., 14(1995), 145-157.

. &Reddy J.N., An introduction to the mathematical theory of finite[53] Oden,elements, John Wiley&Sons, New York, 1976.
$\mathrm{J}.\mathrm{T}$[54] Oishi, S., Two topics in nonlinear system analysis through fixed point theorems,
IEICE Trans. Fund. E77-A (1994), 1144-1153.

[55] –, Numerical verification of existence and inclusion of solutions for nonlinear operator equations, Journal of Computational and Applied Mathematics 60 (1995), 171-185.
[56] –, Numerical analysis with guranteed error bounds of connecting orbits in continuous dynamics, Kokyuroku RIMS Kyoto University 928 (1995), $14- 19$ ( Japanese).
$\mathrm{i}\mathrm{n}$[57] –, Method of numerical inclusion for solutions of fixed point type equations,preprint.
[58] Plum, M., Eigenvalue inclusions for second-order ordinary differential operators by
a numerical homotopy method, ZAMP 41 (1990), 205-226.
[59] –, Computer-assisted existence proofs for two-point boundary value problems, Computing 46 (1991), 19-34.
[60] –, Bounds for eigenvalues of second-order elliptic differential operators, ZAMP 42(1991), 848-863.

[61] –, Existence proofs in combination with error bounds for approximate solutions of
weakly nonlinear second-order elliptic boundary value problems, ZAMM 71 (1991) 6,
T660- T662.

55[62] –, Explicit-estimates and pointwise bounds for solutions of second-order elliptic
boundary value problems, Journal of Mathematical Analysis and Applications 165(1992), 36-61.
$H_{2}$[63] –, Numerical existence proofs and explicit bounds for solutions of nonlinear elliptic
boundary value problems, Computing 49 (1992), 25-44.
[64] –, Enclosures for weak solutions of nonlinear elliptic boundary value problems,
WSSIAA, World Scientific Publishing Company, 3 (1994), 505-521.
[65] –, Existence and enclosure results for continua of solutions of parameter-dependent
nonlinear boundary value problems, Journal of Computational and Applied Mathematics 60 (1995),187-200.

[66] –, Enclosures for two-point boundary value problems near bifurcation points, in), Akademie Verlag,
Scientific Computing and Validated Numerics (G. Alefeld et al.
Berlin (1996), 265-279.
$\mathrm{e}\mathrm{d}\mathrm{s}.$[67] –, Guaranteed Numerical Bounds for Eigenvalues, In: Spectral Theory and Computational Methods of Sturm-Liouville Problems (Eds. D. Hinton and P. W. Schaefer),
Dekker, New York (1997), 313-332.
[68] Rihm, R., Interval methods for initial value problems in ODE’s, in Topics in Validated Computations, Herzberger, J.(ed.), North-Holland, Amsterdam (1994), 173-207.
[69] –, Solving algebraic problems with high accuracy, A new approach to scientificcomputation ( . Kulisch, U. &Miranker, W.L.), Academic Press, New York, 1983.
$\mathrm{e}\mathrm{d}\mathrm{s}$[70] –, Verification methods for dense and sparse systems of equations, in Topics in
Validated Computations, Herzberger, J.(ed.), North-Holland, Amsterdam (1994), 63-

136.

[71] Ryoo, C.-S. &Nakao, M.T., Numerical verification of solutions for variational inequalities, Numerische Mathematik 81 (1998), 305-320.

[72] Ryoo, C.-S., Numerical verifications of solutions for simplified Signorini problems,
to appear in Computers&Mathematics with Applications.
[73] Schr\"oder, J., A method for producing verified results for two-point boundary value
problems, Computing Suppl., 6 (1988), 9-22.
[74] Sinai Ja. G. &Vul, E.B., Discovery of closed orbits of dynamical systems with theuse of computers, J. Statistical Physics, 23 (1980), 27-47.
[75] Stetter, H.J., Validated solution of initial value problems for ODE, Computer arithmetic and self-validating numerical methods(Ullrich, C. ed.), Academic Press, San
Diego, 1990, 171-187.

[76] Tsuchiya, T. &Nakao, M.T., Numerical verification of solutions of parametrized
nonlinear boundary value problems with turning points, Japan Journal of Industrial
and Applied Mathematics 14 (1997), 357-372.
[77] Tsuchiya, T., Numerical verification of simple bifurcation points, Kokyuroku RIMS
Kyoto University 831 (1993), 129-140.
[78] Urabe, M., Galerkin’s procedure for nonlinear periodic systems, Arch. Rat. Mech.
Anal., 20 (1965), 120-152.
[79] –, The Newton method and its application to boundary value problems with nonlinear boundary conditions, in Proc. US-Japan Seminar on Differential &Functional
Equations, 1967, 383-410.

56[80] Watanabe, Y. &Nakao, M.T., Numerical verifications of solutions for nonlinear
elliptic equations, Japan Journal of Industrial and Applied Mathematics 10 (1993),

165-178.

. &Yamamoto, N., Verified computation of solutions[81] Watanabe, Y., Nakao,
equationsrelated to MHD equilibria, Nonlinear Analysis,
elliptic
for nondifferentiable(1997),
28,577-587.
Applications
Methods
Theory,$\mathrm{M}.\mathrm{T}$[82] Watanabe, Y., Yamamoto, N. &Nakao, M.T., A numerical verification method ofsolutions for the Navier-Stokes equations, Reliable Computing 5 (1999), 347-357.

[83] Yamamoto, N. &Nakao, M.T., Numerical verifications of solutions for elliptic equations in nonconvex polygonal domains, Numerische Mathematik 65 (1993), 503-521.
[84] Yamamoto, N. &Nakao, M.T., Numerical verifications for solutions to elliptic equations using residual iterations with higher order finite element, Journal of Computational and Applied Mathematics 60 (1995), 271-279.
[85] Yamamoto, N., A numerical verification method for solutions of boundary value
problems with local uniqueness by Banach’s fixed point theorem, SIAM J. Numer.
Anal. 35 (1998), 2004-2013.
[86] Yamamoto, T., Kanno, S. &Atanassova, , Validated computation of polynomial
zeros by Durand-Kerner method, in Topics in Validated Computations, Herzberger,
J.(ed.), North-Holland, Amsterdam (1994), 27-62.
$\mathrm{L}$[87] Zeidler, E., Nonlinear functional analysis and its applications I, Springer-Verlag,
New York, 1986.