\documentstyle[epsfig,IEEEtran]{article}
\tolerance 10000

\title{Enclosing Solutions of Integral Equations}

\author{Ole Caprani$^1$, Kaj
Madsen$^2$, and Ole Stauning$^2$}

\begin{document}

\pagestyle{myheadings}

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


\maketitle

\auffil{The authors are with the $^1$University of {\AA}rhus
and with the $^2$Institute of Mathematical Modelling,
Technical University of Denmark}
 
\newcommand{\oh}{\mbox{$\mbox{O}$}}
\newcommand{\R}{\mbox{
\begin{picture}(0,0)(0,0)
\put(-4,-0){\makebox(0,0)[lb]{\rm\smash{I}}}
\put(-2,-0){\makebox(0,0)[lb]{\rm\smash{R}}}
\end{picture}$\;\;$}} 
\newcommand{\I}{\mbox{
\begin{picture}(0,0)(0,0)
\put(-4,-0){\makebox(0,0)[lb]{\rm\smash{I}}}
\put(-2,-0){\makebox(0,0)[lb]{\rm\smash{I}}}
\end{picture}$\;$}} 
\newcommand{\N}{\mbox{
\begin{picture}(0,0)(0,0)
\put(-4,-.5){\makebox(0,0)[lb]{\rm\smash{N}}}
\put(-3,-.5){\makebox(0,0)[lb]{\rm\smash{N}}}
\end{picture}$\;\;$}} 
\newcommand{\Z}{\mbox{
\begin{picture}(0,0)(0,0)
\put(-4,-.5){\makebox(0,0)[lb]{\rm\smash{Z}}}
\put(-3,-.5){\makebox(0,0)[lb]{\rm\smash{Z}}}
\end{picture}$\;\;$}} 
\newcommand{\C}{\mbox{
\begin{picture}(0,0)(0,0)
\put(-4,-.5){\makebox(0,0)[lb]{\rm\smash{C}}}
\put(-1.5,-.5){\makebox(0,0)[lb]{\rm\smash{|}}}
\end{picture}$\;\;$}} 
\newcommand{\qed}{ 
 \begin{flushright}
 $\vcenter{\hrule\hbox{\vrule height.6em\kern.6em\vrule}\hrule}$\\ 
 \end{flushright}
} 
\newcommand{\spn}{\mbox{span}} 
\newcommand{\wid}{\mbox{w}} 
\newcommand{\dis}{\mbox{d}} 
\newcommand{\iwid}{\overline{\underline{\mbox{w}}}} 
\newcommand{\Inf}[1]{\mbox{$\underline{#1}$}}       % inf
\newcommand{\Sup}[1]{\mbox{$\overline{#1}$}}        % sup
\newcommand{\Int}[1]{\mbox{$[\Inf{#1},\Sup{#1}]$}}
\newcommand{\ve}[1]{\mbox{$\nabla #1$}}
\newcommand{\fix}[1]{\mbox{$#1^\star$}}
\newcommand{\Bigone}[1]{\mbox{\Large#1\normalsize}}
\newcommand{\bigone}[1]{\mbox{\large#1\normalsize}}


\section{Formulation of the Problem}
In this paper, we will consider the problem of solving 
the integral equation
\begin{equation}
\label{igraleq}
Ty=y,
\end{equation}
where $T$ is an integral operator of the type
\begin{equation}
\label{igral}
Ty(x)=f(x)+\int_a^b g(x,t,y(x),y(t))dt,
\end{equation}
$f$ and $g$ are given smooth funtions, and 
$x$ takes values from a given interval $[a,b]$.  
As a domain of the operator $T$, we will consider the 
set of all functions $y$ 
for which the integral from (\ref{igral}) exists. 
This class will be denoted by ${\cal F}([a,b])$.

We want to provide 
{\it guaranteed bounds} for a solution. In other words, we want to find a 
narrow function interval $\fix{Y}$ which is guaranteed to contain a solution 
$\fix{y}$ of (\ref{igral}).
By a {\em function interval} $Y$, we mean the following set of
functions:
$$Y=\Int{y}=$$
$$\{ y\in{\cal F}([a,b]) |\forall t \in [a,b] 
    : \Inf{y}(t) \leq y(t) \leq \Sup{y}(t)\},$$
where $\Inf{y},\Sup{y}\in{\cal F}([a,b])$.

\section{The Basic Method}
In order to find $\fix{Y}$,
we use the following iterative procedure:
Suppose that we know how to, given a 
function interval $Y$, calculate an enclosure $ITY$ of the set
$$TY=\{Ty\;|\;y\in Y\}.$$
Then, if we start with an abritrary function interval $Y^{(0)}$, then
the iterations
\begin{equation}
\label{iteration}
Y^{(k+1)}=Y^{(k)} \cap ITY^{(k)} \;\;\; k=0,1,...
\end{equation}
define a nested sequence of function intervals with the following 
properties (see \cite{ca:inter}):
\begin{eqnarray}
\label{seqprob1}
\fix{y} \in Y^{(k)} & \Rightarrow & \fix{y} \in Y^{(k+1)} \\
\label{seqprob2}
Y^{(k+1)}=\emptyset & \Rightarrow & \mbox{no solution in } Y^{(k)}\\
\label{seqprob3}
ITY^{(k)} \subseteq Y^{(k)} & \Rightarrow & \mbox{solution exists in }
ITY^{(k)}
\end{eqnarray}

\section{Mean Value Forms}
We base our implementation of the iterative process (\ref{iteration})
on the mean value form of $T$ \cite{ca:mean}:
\begin{equation}
\label{mv}
T_mY=Ty_m+(T'Y)(Y-y_m),
\end{equation}
where as a (real-valued) function $y_m$, we take the lower bound of $Y$:
\begin{equation}
\label{defym}
y_m(x)=\Inf{y}(x).
\end{equation}
As an enclosure $(T'Y)(Y-y_m)$ of the 
first derivative, we take the following expression 
\cite{ca:mean,st}:
\begin{equation}
\label{imean1}
T_mY(x) =
Ty_m(x)\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;
\end{equation}
$$+(Y(x)-y_m(x))\int_a^b G'_3(x,t,Y(x),Y(t))\,dt$$
$$+\int_a^bG'_4(x,t,Y(x),Y(t))(Y(t)-y_m(t))\,dt,$$
where $G'_3$ and $G'_4$ are interval extensions of the partial derivatives
of $g$ with respect to its third and fourth variable, and integrals of
interval-valued functions are defined as in \cite{ca:integ}. 

For every $Y$, $TY\subseteq T_mY$. So, if we find an enclosure $IT_mY$
for $T_mY$, then this enclosure will automatically serve as an
enclosure for $TY$ as well. 

To obtain the enclosure $IT_mY$ for $T_mY$, we enclose 
the integrals by Riemann sums, that correspond to an 
appropriately chosen subdivision 
\begin{equation}
\label{subdiv}
X_k=[x_k,x_{k+1}],\;\;\;k=0,...,n-1,
\end{equation}
where $x_0=a$ and $x_n=b$.
Then, for $x\in X_i$, the value $IT_mY(x)$ of this enclosure takes the
following form:
\begin{equation}
\label{ITmYx}
IT_mY(x) =
Ty_m(x)\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;
\end{equation}
$$+(Y(x)-y_m(x))\sum_{j=0}^{n-1} G'_3(X_i,X_j,Y_i,Y_j)\wid(X_j)$$
$$+\sum_{j=0}^{n-1}G'_4(X_i,X_j,Y_i,Y_j)\int_{X_j}(Y(t)-y_m(t))dt,$$
where $Y_k=Y(X_k)$ is the smallest
(numerical) interval that contains the range $\{Y(x)\;|\;x\in X_k\}$
of the function interval $Y$ over $X_k$.

The only practical problem left here is calculation of the
non-interval term $Ty_m(x)$.

It can be shown that under reasonable assumptions, the distance between $IT_mY$
and $T_mY$ is of the order $\oh(\wid(Y)\max(\Delta,\wid(Y)))$, where 
$\Delta=\max(\wid(X_i))$.
Furthermore, it can be shown
that a contraction property of $T$ is preserved if we extend $T$
to its mean value form, and then enclose the 
interval integrals by Riemann sums.
We can therefore still obtain $Y^{(k)} \rightarrow [\fix{y},\fix{y}]$ for 
$\Delta\rightarrow 0$. Hence, the only source of error will 
be the non-interval term in (\ref{ITmYx}): $Ty_m(x)$ (and rounding errors 
when it comes to real calculations).

Finally, we describe how $Ty_m$ is enclosed. We give two different 
enclosions, using approximations of order zero and of first order.

\section{Enclosing $Ty_m$: Approximation Of Order Zero}
If $x\in X_i$,
then $Ty_m \subseteq I_0T\;y_m$, where
\begin{equation}
\label{zeroigral}
I_0T\;y_m(x)=F(X_i)\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;
\end{equation}
$$+\sum_{j=0}^{n-1} G(X_i,X_j,Y_{m,i},Y_{m,j})\wid(X_j),$$
$$Y_{m,k}=[\inf_{t\in X_k}(y_m(t)),\sup_{t\in X_k}(y_m(t))],$$ 
and $F$ and $G$ are interval extensions of $f$ and $g$.

Thus, the final enclosure of $TY$ is
(\ref{ITmYx}) with the enclosure (\ref{zeroigral})
of $Ty_m$ inserted. 

Notice that if $Y$ is bounded by piecewise constant functions, then the same 
holds for the enclosure.

\section{Enclosing $Ty_m$: First Order Approximation}
Since $y_m$ is a real-valued function (not a variable), when $Ty_m$ is
calculated, we can consider 
$g$ as a function of $x$ and $t$ only: 
\begin{equation}
h(x,t)=g(x,t,y_m(x),y_m(t)).
\end{equation}
If $x\in X_i$,
then $Ty_m \subseteq I_1T\;y_m$, where
\begin{equation}
\label{firstigral}
 I_1T\;y_m(x) = f(x_i) +
F'(X_i)(x-x_i)\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;
\end{equation}
$$+\sum_{j=0}^{n-1} \Bigone{\{} (h(x_i,x_j)+ H'_1(X_i,X_j)(x-x_i))\wid(X_j)$$
$$+ H'_2(X_i,X_j)\cdot \frac{1}{2} \wid(X_j)^2 \Bigone{\}}$$
and $F'$, $H'_1$ and $H'_2$ are interval extensions of the gradients
$f'$, $h'_1$ and $h'_2$.

Thus, the final enclosure of $TY$ is
(\ref{ITmYx}) with the enclosure (\ref{firstigral})
of $Ty_m$ inserted. Notice that if $Y$ is bounded by piecewise linear
functions 
then the same holds for the enclosure.

\section{Convergence Results}
Assume that $T$ is an interval contraction. 
Under reasonable assumptions, we then
obtain:
\begin{itemize}
\item  convergence to a function interval 
$\fix{Y}$ of the width $\oh(\Delta)$, where $\Delta$ is $\max(\wid(X_i))$,
for the approximation of
zero order, and 
\item  convergence to a function interval 
$\fix{Y}$ of the width $\oh(\Delta^2)$
for the first order approximation. 
\end{itemize}

\section{A Numerical Example}
For the numerical calculations, we have used the C++ interval package PROFIL 
\cite{kn:bia}.

Consider the integral equation 
\cite{ra:num}:
\begin{equation}
\label{ambartsumian}
y(x) = 1 + \frac{\lambda}{2} \int_0^1 \frac{x}{x+t}y(x)y(t)dt,
\end{equation}
where $\lambda \in [0,1]$. 

This equation comes from {\it astrophysics}: 
it has been derived by Chandrasekhar
\cite{c} to describe the isotropic scattering of radiation in
plane-parallel atmosphere. 

The results shown in the following tables
correspond 
to $\lambda=\frac{1}{2}$; for this $\lambda$, $T$ is a contraction. 

Notice that our method
also allows an interval to be inserted for $\lambda$.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%AMBARTSUMIAN ZERO
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%DEGREE%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{table}[htb]
\begin{center}
\begin{tabular}{|c||c|c|c|}
\hline\hline
$i$ & $n=15$ & $n=30$ & $n=60$  \\
\hline
0 & 0.41000 & 0.41000 & 0.41000 \\
1 & 0.21141 & 0.16021 & 0.13589 \\
2 & 0.14997 & 0.08580 & 0.05631 \\
3 & 0.13061 & 0.06556 & 0.03552 \\
4 & 0.12437 & 0.06012 & 0.03040 \\
5 & 0.12232 & 0.05864 & 0.02912 \\
6 & 0.12164 & 0.05824 & 0.02880 \\
7 & 0.12141 & 0.05813 & 0.02872 \\
8 & 0.12133 & 0.05810 & 0.02870 \\
9 & 0.12130 & 0.05809 & 0.02869 \\
$\cdots$ & $\cdots$ & $\cdots$ & $\cdots$ \\
$\infty$ & 0.12129 (34) & 0.05809 (30) & 0.02869 (29)\\
\hline\hline
\end{tabular}
\end{center}
\caption{Zero order approximation: The width of $Y^{(i)}$ as a
function of the iteration number $i$ using 
$n=15,30,60$ discretization intervals. The number in parenthesis is
the total number
of iterations before the algorithm stops (when $Y^{(i)}$ becomes
invariant to the iterations).\label{ambar0degtab}}
\end{table}

\begin{table}[htb]
\begin{center}
\begin{tabular}{|c||c|c|c|}
\hline\hline
$i$  & $n=15$ & $n=30$ & $n=60$ \\
\hline
0 & 0.41000 & 0.41000 & 0.41000  \\
1 & 0.16132 & 0.13337 & 0.12205\\
2 & 0.07382 & 0.04560 & 0.03573\\
3 & 0.05370 & 0.02113 & 0.01266\\
4 & 0.05333 & 0.01628 & 0.00725\\
5 & 0.05323 & 0.01569 & 0.00659\\
6 & 0.05321 & 0.01556 & 0.00648\\
7 & 0.05320 & 0.01553 & 0.00646\\
8 & 0.05320 & 0.01552 & 0.00645\\
9 & 0.05319 & 0.01552 & 0.00645\\
$\cdots$ & $\cdots$ & $\cdots$ & $\cdots$ \\
$\infty$  & 0.05319 (33) & 0.01551 (30) & 0.00645 (28) \\
\hline\hline
\end{tabular}
\end{center}
\caption{First order approximation: The width of $Y^{(i)}$ as a
function of the iteration number $i$ using 
$n=15,30,60$ discretization intervals. The number in parenthesis is
the total number
of iterations before the algorithm stops (when $Y^{(i)}$ becomes
invariant to the iterations).\label{ambar1degtab}}
\end{table}

Since the function $g$ has a singularity near $x=t=0$, we use a
non-equidistant subdivision $$x_i = \frac{i^2}{n^2},\;i=0,...,n.$$ 
As the initial function interval, we took 
$Y^{(0)}(x)\equiv [.99,1.4]$. 

Table \ref{ambar0degtab} shows the width of $Y^{(i)}$ as a function of 
the iteration number $i$ using the zero order approximation (\ref{zeroigral}) 
with $n=15, 30, 60$ discretization intervals. 
Table \ref{ambar1degtab} shows the width of $Y^{(i)}$ as a function of 
the iteration number $i$ using the first order approximation 
(\ref{firstigral}) 
with $n=15, 30, 60$ discretization intervals. 

The four figures following this text show the results of the
corresponding iteration processes for $n=15$ and $n=60$: 
\begin{itemize}
\item the first two
figures show the results of using approximation of zero order:
\begin{itemize}
\item[$\bullet$] for $n=15$ (the first figure);
\item[$\bullet$] for $n=60$ (the second figure);
\end{itemize} 
\item the
third and the fourth figures show the results of using first order 
approximation:
\begin{itemize}
\item[$\bullet$] for $n=15$ (the third figure);
\item[$\bullet$] for $n=60$ (the fourth figure).
\end{itemize}
\end{itemize}

All solutions are self-validated, i.e. the final function 
interval is proven to contain a solution to
(\ref{ambartsumian}) because of the inclusion $IT_mY^{(i)} \subseteq Y^{(i)}$.
 
\begin{thebibliography}{99}
  
\bibitem{ca:iter} O. Caprani, K. Madsen, 
``Iterative methods for interval inclusion of fixed points'',
 {\it BIT}, 1978, Vol. 18, pp. 42--51.

\bibitem{ca:inter} O. Caprani, K. Madsen, 
``Interval contractions for the solution of integral equations'',
In: K. Nickel (ed.),  {\bf Interval Mathematics 1980}, Academic Press,
1980, pp. 281--290.  
 
\bibitem{ca:mean} O. Caprani, K. Madsen, 
``Mean value forms in interval analysis'',
{\it  Computing}, 1980, Vol. 25, pp. 147--154.

\bibitem{ca:integ} O. Caprani, K. Madsen, and L. B. Rall,
``Integration of interval functions'',
 {\it SIAM J. Math. Anal}, 1981, Vol. 12, No. 3, pp. 321--341.

\bibitem{c}
S. Chandrasekhar, {\bf Radiative transfer}, Dover, N.Y., 1960.

\bibitem{kn:bia} O. Kn\"{u}ppel,
{\bf PROFIL -- Programmer's Runtime Optimized Fast Interval Library,}
 Technische Universit\"{a}t, Hamburg - Harburg, 1993.

\bibitem{mo:int} R. E. Moore,
{\bf Interval analysis.}
 Prentice-Hall, Englewood Cliffs, NJ, 1966.

\bibitem{mo:met} R. E. Moore,
 {\bf Methods and application of interval analysis,}
 SIAM Stud. App. Math., 1978.  

\bibitem{ra:num} L. B. Rall, 
 ``Numerical integration and the solution of integral 
 equations by the use of Riemann sums,''
 {\it SIAM Review}, 1965, Vol. 7, pp. 55--64. 

\bibitem{st} O. Stauning, 
 {\bf Solving integral equations using interval analysis.}
 Institute of Mathematical Modelling, Technical University of Denmark,
Internal Report, 1994.

\end{thebibliography}

\end{document}
 


