\documentstyle[IEEEtran]{article}
\begin{document}
\tolerance 10000
\title{
A Posteriori Error Estimation and Corrected  Solution
of Partial Differential Equation }
\author{Boris S. Dobronets}

\pagestyle{myheadings}

\markright{APIC'95, El Paso,
Extended Abstracts,
A Supplement to the international journal of {\rm Reliable
Computing}\ \ \ \ \ \ \ \ \ \ \  \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \ \ \ \ \ \ \ \ \ \ \ \ }

\maketitle

\begin{abstract}
In this paper, we show that a posteriori error estimate technique
can be used to produce 
corrected solutions of boundary value problem for partial differential
equations.
\end{abstract}
  
\auffil{The author is with Computer Center of Siberian Department of 
Russian Academy of Sciences,
Akademgorodok,
660036, Krasnoyarsk,
Russia,
E-mail: {\tt dobr@cckr.krasnoyarsk.su}.
This work was partially supported
        by Krasnoyarsk Regional Science Foundation, Krasnoyarsk,
        Russia.}

\section{Formulation Of The Boundary Value Problem}
We will illustrate our idea on the following 
elliptic boundary value
problem
\begin{eqnarray}
\label{ell}
&& Lu = f(x,u), \quad  x\in\Omega ,\\
&& u(x) = 0, \quad x\in\partial\overline\Omega,
\end{eqnarray}
where $\Omega$ is a bounded open convex domain in $R^2$, with piecewise
smooth boundary $\partial\overline\Omega$, and
$$
Lu = \sum_{i=1}^{2} \partial^2 u/\partial x_i^2.
$$
In this paper,
we assume that
$$
\partial f(x,u) / \partial u \ \geq \ q(x) \geq 0,
$$
and that a positive constant $K$ exists such that
$$
|f(x,\eta)| \leq K(1+| \eta |),$$ for all $x \in \Omega$ ,
and for all $\eta \in [ \min_\Omega  u, \max_\Omega u].$

The solution of (\ref{ell}) is understood in the following weak sense:
find a function $u\in W_2^1 (\Omega)$ such that
\begin{equation}
L(u,v) = (f,v), \quad\forall v\in W_2^1 (\Omega),
\end{equation}
where $(\cdot,\cdot)$ is the inner product in $L_2$:
$$
L(u,v) = \int_\Omega \sum_{i=1}^2
\partial_i u\,\partial_i v\,\mbox{d}\Omega.
$$

\section{A Numerical Solution Of The Boundary Value Problem}
To describe the solution, we must first define 
the finite element space $S^n_l$. We fix a mesh size $h$, and 
divide the domain $\Omega$ into
subdomains of size $h$. The set of these subdomains will be denoted by
${\cal T}_h$. The finite-element space is defined as a set of
piece-wise polynomial functions (i.e., function that are polynomial on
each of the subdomains $T\in {\cal T}_h$):
\begin{equation}
S^n_l = \{\, s(x) \mid s\in W_2^l (\overline\Omega),
s|_T \in{\cal P}^n, T\in{\cal T}_h \},
\end{equation}
where ${\cal P}^n$ is the set of polynomials of degree $n$.

We define the finite element solution $u^h$ of the problem (\ref{ell}) as
a function from $S^1_1$ that satisfies the conditions
\begin{eqnarray}
&& L(u^h,v^h) = (f(x,u^h(x)),v^h), \forall v^h \ \mbox{in} \ S^1_1, \\[2mm]
&& u_h = 0,\quad \mbox{on} \ \partial\overline\Omega. \nonumber
\end{eqnarray}
Since the function $u^h$ is uniquely determined by its (finitely
many) polynomial coefficients, we have a (non-linear) system of
equations, that can be solved, e.g., if we use iterative techniques. 

\section{Estimating Error Of This Numerical Solution}
We would like to use the {\it defect} $$Lu^h-f(x,u^h)$$ of the approximate
solution
$u^h$ to 
estimate its accuracy, i.e., the difference between this approximate
solution and the (unknown) precise one. However, we cannot compute the
defect of the original solution, because by definition, the function
$u^h$ is differentiable only once, so, we cannot compute its second
derivatives that form an operator $L$. 
  
So, to apply this idea, we use 
the finite element solution $u^h$ to construct
a new smooth approximate solution $s\in S_2^3 (\Omega)$.
Informally speaking, this solution $s$ must satisfy two conditions:
\begin{itemize}
\item First, it must be close to the original solution $u^h$; i.e.,
for every mesh point $v_i$, we must have $s(v_i)\approx u^h(v_i)$.
\item Second, this function $s$ must be an (approximate) solution of
the equations (1). This means, in particular, that for every mesh
point $v_i$, we must have $Ls(v_i)\approx f(x,s(v_i))$. The function
$f$ is non-linear, so, this condition is non-linear in $s$. However,
we can linearize it if we take into consideration that $s(v_i)\approx
u^h(v_i)$. Then, this condition can be reformulated as follows:
$Ls(v_i)\approx f(v_i,u^h(v_i))$.
\end{itemize}
We cannot exactly satisfy both sets of conditions, so we must apply a
techniques similar to least squares method to combine them. 
As a result, we arrive at the following method of determining $s$:
For an arbitrary mesh point 
$x_0 = (x^0_1, x^0_2)$, we define the (local) polynomial solution 
$p$ as a polynomial $$p =
\sum_{i,j: 0\le i+j\le n}a_{i,j} (x_1 -x^0_1)^i (x_2 -
x^0_2)^j$$ such that
$$\sum_{i=1}^N \alpha_i| p(v_i) - u^h (v_i) |^2 +$$
\begin{equation}
\label{p}
\sum_{i=1}^M \beta_i| Lp(w_i) - f(w_i,u^h(w_i)) |^2 \rightarrow \min ,
\end{equation}
where $\{v_i\}_{i=1}^{N}$ 
and $\{w_i \}_{i=1}^{M}$ are nodes of the mesh,
$$
\alpha_i = a_1 / (\rho(v_i,x_0) + a_2),
$$
$$
\beta_i = b_1 / (\rho(w_i,x_0) + b_2),
$$
$a_i$ and $b_i$ are positive constants, and 
$\rho(x,y)$ denotes the distance between the points $x,y$. 

The problem
(\ref{p}) is quadratic in the unknown coefficients $a_{ij}$. Thus, 
if we differentiate with respect to $a_{ij}$, and equate the
derivatives to 0, we thus reduce this problem to 
the solution of the following system of linear algebraic
equations:
$$
Ba = d,
$$
where
$$
a = (a_{0,0}, a_{1,0}, a_{0,1}, \ldots,a_{n,0},a_{n-1,1},\ldots,a_{0,n} ),$$
$$d=(u^h_1,\ldots,u^h_N,$$
$$f(w_1,u^h(w_1)),\ldots,f(w_M,u^h(w_M))).$$

After this smooth approximate solution $s$ is computed, 
we can use the following defect \cite{DoS}
\begin{equation}
\phi (x,s) = Ls - f(x,s), \quad x\in\Omega.
\end{equation}
To get the desired interval estimate, we apply the finite element
method (FEM) to solve
numerically the following additional problem:
\begin{eqnarray}
\label{el1}
&& L_1 u_1 = 1,\quad x\ \mbox{in} \ \Omega,\\
&& u_1(x) = 0,\quad x \ \mbox{on} \ \partial\overline\Omega,
\end{eqnarray}
where $L_1 = Lu - qu$.
We solve this auxiliary problem by constructing a special spline $s_1$.

Then, the desired interval solution has the form
$$
{\bf u} = s + [\underline\alpha , \overline\alpha] s_1 +
[\underline\beta , \overline\beta],
$$
where
$$ {\bf \overline\alpha} =\max_{\overline\Omega}
 (\phi / (L_1 s_1), 0),\ 
 {\bf \underline\alpha}=\min_{\overline\Omega}
 (\phi / (L_1 s_1), 0),
$$
$$
{\bf \overline\beta} = \max_{\partial\Omega}
 (-\overline\alpha s_1 - s, 0),
\ {\bf \underline\beta}=\min_{\partial\Omega}
 (-\underline\alpha s_1 - s, 0),
$$

\section{How To Correct The Solution Of The Boundary Value Problem}
To correct the solution, we must somehow find a good approximation for
the difference $u-s$ between the exact (unknown) and the approximate
solutions. To estimate this difference, let us find a differential
equation that this function $\varepsilon=u-s$ must satisfy. 

If we substitute $u=s+\varepsilon$ into the equation (1), we get the
following equation:
$$Ls+L\varepsilon=f(x,s+\varepsilon),$$ or
$$L\varepsilon=f(x,s+\varepsilon)-Ls.$$ To estimate $\varepsilon$, we
can esither solve this non-linear equation (using the same techniques
as above), or simplify it.
We can simplify this 
non-linear equation if 
we take into consideration the fact that 
$\varepsilon$ is small and therefore, we can expand $f$ into a power
series in $\varepsilon$ and keep only linear terms. As a result, we
get the following equation:
$$L\varepsilon=f(x,s)+r(x)\varepsilon-Ls,$$ where we denoted
$$r(x)={\partial f\over\partial u}(x,s(x)).$$
If we move all the terms that contain $u$ into the left-hand side, 
and all the other terms into the right-hand side, we 
get the following equation:
$$L\varepsilon-r(x)\varepsilon=\psi(x),$$
where $\psi(x)=f(x,s(x))-(Ls)(x).$ 

Let $\varepsilon^h$ be a solution of this approximate equation, obtained by a
similar finite-element technique. Then, for the corrected solution
$$s_{cor}=s+\varepsilon^h,$$ the following result is true:
\newpage

\noindent
{\bf THEOREM.} {\it If $u\in W_2^4 (\Omega)$, then
$$
\| u - s_{cor}\|\leq C\cdot h^4 \| u \|_{W_2^4(\Omega)} ,
$$
where $h$ is the mesh size, and $C$ is a constant independent of $h$.}
\smallskip

To get an even better solution, we can apply the correction procedure
to the corrected solution; then, if necessary, we can apply the
correction procedure again and again, and get better and better
solutions. 

In principle, instead of starting with an approximate solution, we 
can start with the 
$u^h=0$ and $s=0$. Thus, we will get a new iterative method of
numerical solution for such equations.

\begin{thebibliography}{99}
\bibitem{DoDC} B. S. Dobronets, ``Difference correction  by  finite  elements
         of higher order'', In: {\bf Mathematical  models and
methods of solving the problems of
          continuous media}, Krasnoyarsk, 1986, pp. 80--88 (in Russian).
\bibitem{DoS} B. S. 
Dobronets, V. V. Shaydurov, {\bf Two-sided Numerical Methods},
         Nauka Publ. (Siberian Branch), Novosibirsk, 1990 (in Russian).
\bibitem{Streng} G. Streng  and 
G. F. Fix, {\bf An Analysis of the Finite Element
         Method}, Englewood Cliffs, Prentice-Hall, 1973.
\end{thebibliography}
\end{document}