\documentstyle[IEEEtran]{article}
\begin{document}
\tolerance 10000
\title{               A Method of Estimating the Solution
                of a System of Ordinary Differential Equations
by Using Interval Taylor Series}
\author{ V. S. Zyuzin and E. A. Novikova}

\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 Department of Mechanics and Mathematics,
 Saratov State University,
 Astrakhanskaya 83,
 Saratov, 410071, Russia,
 e-mail zyuzin@scnit.saratov.su.}

\section{Ellipsoid Method for Estimating the Accuracy of the Solution
of a System of ODE}
     In \cite{1,2},  Filippov  proposed the following {\it ellipsoid
method} for
estimating the error
of the approximate solution of a system of ordinary differential
equations (ODE): In this method, 
at any moment of time $t$, we have an ellipsoid $E(t)$ that is guaranteed to 
contain the actual value of the solution $x(t)=(x_1(t),...,x_n(t))$.
The initial ellipsoid represents the initial uncertainty. 
The equations that describe how this ellipsoid changes in time are
obtained from the original system of ODE, if we take into account 
rounding errors, and additional errors caused by the method itself.

The evolution equations can be then solved using 
an appropriate iterative technique. 
In particular, 
Filippov uses Adams' method.  

\section{Taylor Series Techniques Leads to Better Estimates for an
Ellipsoid Method}
Instead of Adams' method, we used Taylor series techniques to solve the
evolution equations for the ellipsoid. As a result, for several
reasonable systems of ODE, we got better error estimates.

\section{We Can Get Even Better Error Estimates If We Use
Parallelepipeds Instead of Ellipsoids}
We implemented this method using 
a PASCAL-SC compiler \cite{3} and a PASCAL-XSC compiler
\cite{4}. It turned out 
after processing our program, these compilers actually ended up
with algorithm in which the parameters of the ellipsoid $E(t)$ were
computed as follows:
\begin{itemize}
\item first, we compute a changing {\it
parallelepiped} that contains the actual solution of the system, and
then 
\item compute an ellipsoid that contains this parallelepiped.
\end{itemize}
Because of that, we get a better 
estimate if we do not use ellipsoids at all, and instead, use 
a changing {\it parallelepiped} to estimate the errors of the approximate
solution.
Thus, we arrive at the following method:

\section{Mathematical Formulation of the Problem}
     Let the system of ordinary differential equations
$$x'(t)=f(t,x(t)) \eqno{(1)}$$
with the initial conditions
                        $$x(t_0)=x_0 \eqno{(2)}$$
be given, where 
$x(t)$ and $f(t,x)$ are functions whose values are 
vectors in $n-$dimensional Euclidean space.

Let us assume that for every initial value $x(t_0)$, the system (1) 
has a unique solution $x(t)$.

\section{The Idea Behind Our Method}
     Let $\tilde x(t)$  be an approximate solution to the problem
(1)--(2). Usually, it is  a vector function, that is obtained by some simple
interpolation procedure (piecewise constant or piecewise linear)
from its values for $t_i=t_0+ih$, $i=1,..,n$ that are computed during
the numerical solution (e.g., if we apply the Runge-Kutta method). 

To describe how the parallelepiped evolves in time, i.e., how much error
this method can add on the interval between two consequent points, we
we will define the notion of a {\it local error}. Suppose that at the moment
$t_{i-1}$, we know the value of $x(t)$ precisely, i.e, that
$x(t_{i-1})=\tilde x(t_{i-1})$. Since we assumed that the solution of 
our system (1) is uniquely determined by the initial conditions, the
condition $x(t_{i-1})=\tilde x(t_{i-1})$ uniquely determines a solution; we
will denote this solution by $x_{i-1}(t)$. As $t$ goes from $t_{i-1}$
to $t_i$, 
the value $x_{i-1}(t)$ of this solution may differ
from $\tilde x(t)$. 
This difference $l(t)=\tilde x(t_i)-x_{i-1}(t_i)$ will be called the
{\it local error}.

 The local error can be estimated as 
a sum of the error of the approximation method
(that is usually known) and the rounding error.  

Our goal is to estimate the difference 
between the exact solution
 and the approximate solution of ODE system (1):
$$e(t)=x(t)-\tilde x(t) \eqno{(3)}.$$
For $t$ from $t_{i-1}$ to $t_i$, the desired difference $e(t)$
 can be  represented  as  a difference $$e(t)=y(t)-l(t)$$ 
 between a function $$y(t)=x(t)-x_{i-1}(t)\eqno{(4)}$$
and the local error $l(t)$.
Since we already know the bounds for local error, to estimate $e(t)$,
it is sufficient to be able to estimate $y(t)$. 

Since both $x(t)$ and $x_{i-1}(t)$ are solutions of the equation (1),
after differentiating both sides of the equation (4), we get 
the following equation:
               $$y'(t)=x'(t)-x_{i-1}'(t)=f(t,x(t))-f(t,x_{i-1}(t)).$$
Substituting $x(t)=x_{i-1}(t)+y(t)$ into this equation, we 
obtain the first order differential equation for $y$:
$$y'(t)=f(t,x_{i-1}(t)+y(t))-f(t,x_{i-1}(t)).\eqno{(5)}$$

The initial value for $y$ is $y(t_{i-1})=x(t_{i-1})-x_{i-1}(t_{i-1})$.
By definition of $x_{i-1}$, this difference is equal to 
$y(t_{i-1})=x(t_{i-1})-\tilde x(t_{i-1})$.
So, for $i=0$, we have $y(t_0)=0$, and for further $i$, as an estimate
for this initial value, we can take the error bound that was obtained on the
previous step of the error estimation.

To solve the system (5), we use the method of Taylor series \cite{5}
(we use Lohner's method to take the remainder into account).

\section{Testing}
     The method of  parallelepipeds  was  successfully 
tested  on  two  simple  ODE
 systems:
$$x_1'(t)=x_1-x_2,\ x_2'(t)=2x_1-x_2^3,$$ $$x_1(0)=0.25,\ x_2(0)=0,$$
and
$$x_1'(t)=-x_2+x_1\cdot(x_1^2+x_2^2+1),$$ 
$$x_2'(t)=x_1+x_2\cdot (x_1^2+x_2^2+1),$$ 
$$x_1(0)=0.5,\ x_2(0)=0.$$

\section{Potential Applications}
We are currently applying 
our method to estimating
the phase state of dynamical systems \cite{7}, in particular, to
intertial navigation problems.

\begin{thebibliography}{99}

\bibitem{1} A. F. Filippov, 
``Ellipsoidal estimates for a solution of a system of
    differential equations'', {\it Interval Computations}, 1992, No.
2(4), pp. 6--16.

\bibitem{2} A. F. Filippov,  
``Ellipsoidal  error  estimates  for Adams method'',
    {\it Interval Computations}, 1992, No. 3(5), pp. 75--79.

\bibitem{3} U. Kulisch (ed.),  {\bf 
PASCAL-SC:  A PASCAL Extension for  Scientific
    Computation,  Information  Manual  and  Floppy Disks,  Version IBM
    PC/AT (DOS)}, Teubner, Stuttgart, 1987.

\bibitem{4} R. Klatte, U. Kulisch, M. Neaga,  D. Ratz, and  Ch. Ullrich,
{\bf PASCAL-XSC,  Sprachbeschrebung mit Biespielen},  Springer  Verlag,
    1991.

\bibitem{5} R. E. Moore, {\bf 
Methods and applications of interval analysis}, SIAM,
    Philadelphia, 1979.

\bibitem{6} V. S. Zyuzin  and I. V. Leskina, ``Ellipsoidal  
Estimates  via Taylor's
    Series for an Approximate Solution of  a  System  of  Differential
    Equations'', {\it SCAN-93}, Vienna, Austria, 1993.

\bibitem{7} 
F. L. Chernousko, {\bf Estimating the  phase  state  dynamical  system},
    Nauka, Moscow, 1988 (in Russian).

\end{thebibliography}


\end{document}