Some recent advances in validated methods for IVPs for ODEs.
Compared to standard numerical methods for initial value problems (IVPs) for ordinary differential equations (ODEs), validated methods (often called interval methods) for IVPs for ODEs have two important advantages: if they return a solution to a problem, then (1) the problem is guaranteed to have a unique solution, and (2) an enclosure of the true solution is produced.We present a brief overview of interval Taylor series (ITS) methods for IVPs for ODEs and discuss some recent advances in the theory of validated methods for IVPs for ODEs. In particular, we discuss an interval Hermite-Obreschkoff (IHO) scheme for computing rigorous bounds on the solution of an IVP for an ODE, the stability of ITS and IHO methods, and a new perspective on the wrapping effect, where we interpret the problem of reducing the wrapping effect as one of finding a more stable scheme for advancing the solution.
Preprint submitted to Elsevier Preprint 20 November 2000
is D. The condition (2) permits the initial value
to be in an interval, rather than specifying a particular value. We assume
that the representation of f contains a nite number of constants, variables,
elementary operations, and standard functions. Since we assume f 2 C k 1 (D),
we exclude functions that contain, for example, branches, abs, or min. For
expositional convenience, we consider only autonomous systems. This is not
a restriction of consequence, since a nonautonomous system of ODEs can be
converted into an autonomous one. Moreover, the methods discussed here can
be extended easily to nonautonomous systems.
We consider a grid t 0 < which is not necessarily equally spaced,
and denote the stepsize from t j 1 to t j by h . The step from
is referred to as the jth step. We denote the solution of (1) with an
initial condition For an interval, or
an interval vector [y j 1 ], we denote by y(t; t the set of solutions
Our goal is to compute interval vectors [y j m, that are guaranteed
to contain the solution of (1{2) at
Standard numerical methods for IVPs for ODEs attempt to compute an approximate
solution that satises a user-specied tolerance. These methods are
usually robust and reliable for most applications, but it is possible to nd
examples for which they return inaccurate results. On the other hand, if a validated
method for IVPs for ODEs returns successfully, it not only produces a
guaranteed bound for the true solution, but also veries that a unique solution
to the problem y exists for all y
There are situations when guaranteed bounds are desired or needed. For ex-
ample, a guaranteed bound on the solution could be used to prove a theorem
[28]. Also, some calculations may be critical to the safety or reliability of a
system. Therefore, it may be necessary or desirable to ensure that the true
solution is within the computed bounds.
One reason validated solutions to IVPs for ODEs have not been popular in
the past is that their computation typically requires considerably more time
and memory than does that of standard methods. However, now that \chips
are cheap", it seems natural to shift the burden of determining the reliability
of a numerical solution from the user to the computer | at least for standard
problems that do not require a very large amount of computer resources.
In addition, there are situations where interval methods for IVPs for ODEs
may not be computationally more expensive than standard methods. For ex-
ample, many ODEs arising in practical applications contain parameters. Often
these parameters cannot be measured exactly, but are known to lie in certain
intervals, as, for example, in economic models or in control problems. In these
situations, a user might want to compute solutions for ranges of parameters. If
a standard numerical method is used, it has to be executed many times with
dierent parameters, while an interval method can \capture" all the solutions
at essentially no extra cost.
Important developments in the area of validated solutions of IVPs for ODEs
are the interval methods of Moore [19], Kruckeberg [13], Eijgenraam [8], and
Lohner [17]. All these methods are based on Taylor series. One reason for the
popularity of this approach is the simple form of the error term. In addition,
the Taylor series coe-cients can be readily generated by automatic dierenti-
ation, and both the order of the method and its stepsize can be changed easily
from step to step.
Usually, validated methods for IVPs for ODEs are one-step methods, where
each step consists of two phases:
Algorithm I: validate existence and uniqueness of the solution with some
stepsize, and
Algorithm II: compute a tight enclosure for the solution.
The main di-culty in Algorithm I is how to validate existence and uniqueness
of the solution with a stepsize that corresponds to the order of Algorithm
II. The main obstacle in Algorithm II is how to reduce the so-called
wrapping eect, which arises when the solution set, which is not generally a
box, is enclosed in | or wrapped by | a box on each integration step, thus
introducing overestimations in the computed bounds. Currently, Lohner's QR-
factorization method is the most-eective standard scheme for reducing the
wrapping eect.
The methods considered in this paper are based on high-order Taylor series
expansions with respect to time. Recently, Berz and Makino [4] proposed a
method for reducing the wrapping eect that employs high-order Taylor series
expansions in both time and the initial conditions. Berz and Makino's scheme
validates existence and uniqueness and also computes tight bounds on the
solution at each grid point t j in one phase, while the methods considered here
separate these tasks into two phases, as noted above.
The purpose of this paper is to present a brief overview of interval Taylor
series (ITS) methods for IVPs for ODEs and to discuss some recent advances
in the area of validated ODE solving. In particular, we discuss
an interval Hermite-Obreschko (IHO) scheme for computing tight enclosures
on the solution [20, 21];
instability in interval methods for IVPs for ODEs due to the associated
formula for the truncation error [20, 21], which appears to make it di-cult
to derive eective validated methods for sti problems; and
a new perspective on the wrapping eect [22], where we view the problem
of reducing the wrapping eect as one of nding a more stable scheme for
advancing the solution.
Section 2 introduces interval-arithmetic operations, enclosing ranges of func-
tions, and automatic generation of Taylor series coe-cients. In Section 3, we
brie
y discuss methods for validating existence and uniqueness of the solution,
explain the wrapping eect, and describe Lohner's QR-factorization method
for reducing the wrapping eect.
In Section 4, we outline the IHO method. Section 5 shows that the stability of
the ITS and the IHO methods depends not only on the formula for advancing
the solution, as in point methods for IVPs for ODEs, but also on the associated
formula for the truncation error. Section 6 explains how the wrapping eect
can be viewed as a source of instability in interval methods for IVPs for ODEs
and how Lohner's QR-factorization scheme improves the stability of an interval
method.
The set of intervals on the real line R is dened by
a
If a
a then [a] is a point interval; if a
0 then [a] is nonnegative ([a]  0);
and if a
then [a] is symmetric. Two intervals [a] and [b] are equal if a
and  a =  b.
Let [a] and [b] 2 IR, and - 2 =g. The interval-arithmetic operations
are dened [19] by
which can be written in the following equivalent form (we omit  in the nota-
a
a
b;
a b
min
a
b;  ab
a
b;  ab
1=  b; 1=b
The strength of interval arithmetic when implemented on a computer is in
computing rigorous enclosures of real operations by including rounding errors
in the computed bounds. To include such errors, we round the real intervals
in (4{7) outwards, thus obtaining machine interval arithmetic. For example,
when adding intervals, we round a
+b  down and round  a+  b up. For simplicity
of the discussion, we assume real interval arithmetic in this paper. Because
of the outward roundings, intervals computed in machine interval arithmetic
always contain the corresponding real intervals.
Denition (3) and formulas (4{6) can be extended to vectors and matrices. If
the components of a vector or matrix are intervals, we have interval vectors
and matrices, respectively. The arithmetic operations involving interval vectors
and matrices are dened by the standard formulas, except that reals are
replaced by intervals and real arithmetic is replaced by interval arithmetic in
the associated computations.
We have an inclusion of intervals
[a]  [b] () a
and  a   b:
The interval-arithmetic operations are inclusion monotone. That is, for real
intervals [a], [a 1 ], [b], and [b 1 ] such that [a]  [a 1 ] and [b]  [b 1 ],
Although interval addition and multiplication are associative, the distributive
law does not hold in general [1]. That is, we can easily nd three intervals [a],
[b], and [c], for which
However, for any three intervals [a], [b], and [c], the subdistributive law
does hold. Moreover, there are important cases in which the distributive law
does hold. For example, it holds if [b] [c]  0, if [a] is a point interval, or if [b]
and [c] are symmetric. In particular, for  2 R, which can be interpreted as
the point interval [; ], and intervals [b] and [c], we have
For an interval [a], we dene width and midpoint of [a] as
a a
and
respectively [19]. Width and midpoint are dened component-wise for interval
vectors and matrices.
Using (4{6) and (9), one can easily show that
for any [a], [b] 2 IR and  2 R.
The equalities (11) and (12) also hold when [a] and [b] are interval vectors.
If A is an n  n real matrix and [a] is an n dimensional interval vector, then
from (11{12),
where jAj is obtained by taking absolute values on each component of A.
2.1 Ranges of Functions
R be a function on D  R n . The interval-arithmetic evaluation
of f on [a]  D, which we denote by f ([a]), is obtained by replacing each
occurrence of a real variable with a corresponding interval, by replacing the
standard functions with enclosures of their ranges, and by performing interval-
arithmetic operations instead of the real operations. From (8), the range of f ,
g, is always contained in f ([a]). Although f ([a]) is not unique,
since expressions that are mathematically equivalent for scalars, such as x(y
z) and xy may have dierent values if x, y, and z are intervals, a value
for f ([a]) can be determined from the code list, or computational graph, for
f . This implies that, if we rearrange an interval expression, we may obtain
tighter bounds.
continuously dierentiable on D  R n and [a]  D, then,
for any y and b 2 [a], some  2 [a] by the
mean-value theorem. Thus,
[19]. The mean-value form, f M ([a]; b), is popular in interval methods since
it often gives tighter enclosures for the range of f than the straightforward
interval-arithmetic evaluation of f itself. Like f , f M is not uniquely dened,
but a value for it can be determined from the code lists for f and f 0 .
2.2 Automatic Generation of Taylor Coe-cients
Since the interval methods considered here use Taylor series coe-cients, we
outline the scheme for their generation. More details can be found in [3] or
[19], for example.
If we know the Taylor coe-cients (u) i  u (i) =i! and (v) i  v (i) =i! for
for two functions u and v, then we can compute the pth Taylor
coe-cient of u  v, uv, and u=v by standard calculus formulas [19].
We introduce the sequence of functions

@y
f
For the IVP y the ith Taylor coe-cient of y(t) at t j ,
where
is the (i 1)st Taylor coe-cient of f evaluated at y j . Using
and formulas for the Taylor coe-cients of sums, products, quotients, and
the standard functions, we can recursively evaluate (y j  1. It can be
shown that, to generate k Taylor coe-cients, we need at most O(k 2 ) times as
much computational work as we require to evaluate f(y) [19]. Note that this is
far more e-cient than the standard symbolic generation of Taylor coe-cients.
If we have a procedure to compute the point Taylor coe-cients of y(t) and
perform the computations in interval arithmetic with [y j ] instead of y j , we
obtain a procedure to compute the interval Taylor coe-cients of y(t) at t j .
3 Overview of Interval Methods for IVPs for ODEs
3.1 Algorithm I: Validating Existence and Uniqueness of the Solution
Using the Picard-Lindelof operator and the Banach xed-point theorem, one
can show that if h j 1 and [~y 0
then (1) with has a unique solution y(t;
We refer to a method based on (16) as a rst-order enclosure method. Such a
method can be easily implemented, but a serious disadvantage of this approach
is that it often restricts the stepsize Algorithm II could take.
One can obtain methods that enable larger stepsizes by using polynomial
enclosures [18] or more Taylor series terms in the sum in (16), thus obtaining
a high-order Taylor series enclosure method [6, 24]. In the latter, we determine
holds for all
has a unique solution y(t;
for all
In [24], we show that, for many problems, an interval method based on (17)
is more e-cient than one based on (16).
3.2 Algorithm II: Computing a Tight Enclosure of the Solution
Consider the Taylor series expansion
its lth component
n) evaluated at Using the a
priori bounds [~y we can enclose the local truncation
error of the Taylor series (18) on
We can also bound the ith Taylor coe-cient f [i] (y
Therefore,
contains However, as explained below, (19) illustrates how replacing
real numbers in an algorithm by intervals often leads to large overestimations

Taking widths on both sides of (19) and using (11), we obtain
w([y
which implies that the width of [y j ] almost always increases 2 with j, even if
the true solution contracts.
A better approach is to apply the mean-value theorem to f [i] in (18), obtaining
for any ^
I
An equality is possible only in the trivial cases h
for all
where J
is the Jacobian of f [i] with its lth row evaluated at
This formula is
the basis for the ITS methods of Moore [19], Eigenraam [8], Lohner [17], and
Rihm [27]; see also [23].
Let
and
(The ith Jacobian can be computed by generating the Taylor coe-cient
dierentiating it [2, 3]. Alternatively, these Jacobians
can be computed by generating the Taylor coe-cients for the associated variational
equation [17].)
Using (21), we can rewrite (20) as
We refer to a method implementing (22) as the direct method.
If we compute enclosures with (22), the widths of the computed intervals may
decrease. However, this approach frequently works poorly, because the interval
vector
signicantly overestimate the set
Such overestimations accumulate as the integration proceeds. This is often
called the wrapping eect.
In

Figure

1, we illustrate the wrapping of the parallelepiped f Ax j x 2 [x] g
by the box A[x], where
A =B @ 1 2
Fig. 1. The wrapping of the parallelepiped f Ax j x 2 [x] g by the box A[x].
3.2.1 The Wrapping Eect
The wrapping eect is clearly illustrated by Moore's example [19],
The interval vector [y 0 ] can be viewed as a box in the (y 1
the true solution of (23) is the rotated box shown in Figure 2. If we want to
Fig. 2. The rotated box is wrapped at
enclose this box in an interval vector, we have to wrap it by another box with
sides parallel to the y 1 and y 2 axes. On the next step, this larger box is rotated
and so must be enclosed in a still larger one. Thus, at each step, the enclosing
boxes become larger and larger, but the true solution set is a rotated box of
the same size as [y 0 ]. Moore showed that at 2, the interval inclusion is
in
ated by a factor of e 2  535 as the stepsize approaches zero [19].
3.2.2 Lohner's Method
Here, we describe Lohner's QR-factorization method [16, 17], which is one
of the most successful, general-purpose methods for reducing the wrapping
eect.
Let
(j  1), where I is the identity matrix. From (20{21) and (24{25), we compute
and propagate for the next step the interval vector
where A j 2 R nn is nonsingular for
In Lohner's QR-factorization method,
where Q j is orthogonal and R j is upper triangular. Other choices for A j are
discussed in [16, 17, 20, 22].
One explanation why this method is successful at reducing the wrapping eect
is that we enclose the solution on each step in a moving orthogonal coordinate
system that \matches" the solution set; for more details, see [16, 17, 20, 23].
In section 6, we show that this choice for A j ensures better stability of the
QR-factorization method, compared to the direct method.
4 An Interval-Hermite Obreschko Method
For more than Taylor series has been the only eective approach for
computing rigorous bounds on the solution of an IVP for an ODE. Recently,
we developed a new scheme, an interval Hermite-Obreschko (IHO) method
[20, 21]. Here, we outline the method and its potential.
Let
c q;p
(q; p, and i  0), y It can be shown
that
c p;q
[20, 31]. The Hermite-Obreschko formula (29) is the basis of our IHO meth-
od. As in an ITS method, we can easily bound the local truncation error in
the IHO method, since we can readily generate the interval Taylor coe-cient
The method we propose in [20] consists of two phases, which can be considered
as a predictor and a corrector. The predictor computes an enclosure [y (0)
of
the solution at t j , and using this enclosure, the corrector computes a tighter
enclosure
j ] at t j . If q > 0, (29) is an implicit scheme. The corrector
applies a Newton-like step to tighten [y (0)
We have shown in [20, 21] that for the same order and stepsize, our IHO
method has smaller local error, better stability, and requires fewer Jacobian
evaluations than an ITS method. The extra cost of the Newton step is one
matrix inversion and a few matrix multiplications.
5 Instability from the Formula for the Truncation Error
In this section, we investigate the stability of the ITS and IHO methods, when
applied with a constant stepsize and order to the test problem
where  and y 0 2 R, and  < 0. Since we have not dened complex interval
arithmetic, we do not consider problems with  complex. Note also that the
wrapping eect does not occur in one-dimensional problems.
For the remainder of this paper, we consider ITS methods with a constant
stepsize and order k truncation error and in this section, we consider the IHO
scheme with a constant stepsize and order
5.1 The Interval Taylor Series Method
Suppose that at t j 1 > 0, we have computed a tight enclosure [y ITS
of the
solution of (30) with an ITS method, and [~y ITS
is an a priori enclosure of the
solution on
Let
r
Using (31), an ITS method for computing tight enclosures of the solution to
(30) can be written as
[y ITS
[~y ITS
cf. (20). Since [~y ITS
we obtain from (31{32)
that
w([y ITS

w([y ITS
(Note that w([y ITS
Therefore, if
the ITS method given by (32) is asymptotically unstable, in the sense that
lim j!1 w([y ITS
This result implies that we have restrictions on the
stepsize not only from the function T k 1 (h), as in point methods for IVPs
for ODEs, but also from the factor jhj k =k! in the remainder term.
5.2 The Interval Hermite-Obreschko Method
we assume that [y IHO
computed with an IHO
method, and [y IHO
. The formula (29) reduces to
c p;q
are dened in (28). Dene
R p;q
and Q p;q
c q;p
Also let [~y IHO
be an a priori enclosure of the solution on
From (34{35), we compute an enclosure [y IHO
[y IHO
[~y IHO
where
w([y IHO

w([y IHO
Therefore, the IHO method is asymptotically unstable in the sense that
lim j!1 w([y IHO
In (33) and (36),
are approximations to e z of the same order. In particular, R p;q (z) is the Pade
rational approximation to e z (see for example [26]). For the ITS method,
only if h is in the nite stability region of T k 1 (z).
However, for the IHO method with 0 >  2 R, jR p;q (h)j < 1 for any h > 0 if
Roughly speaking, the
stepsize in the ITS method is restricted by both
while in the IHO method, the stepsize is limited mainly by
In the latter case,
p;q =Q p;q (h) is usually much smaller than one; thus, the
stepsize limit for the IHO method is usually much larger than for the ITS
method.
An important point to note here is that an interval version of a standard
numerical method, such as the Hermite-Obreschko formula (29), that is suitable
for sti problems may still have a restriction on the stepsize. To obtain an
interval method without a stepsize restriction, we must nd a stable formula
not only for the propagated error, but also for the associated truncation error.
6 A New Perspective on the Wrapping Eect
The problem of reducing the wrapping eect has usually been studied from
a geometric perspective as nding an enclosing set that introduces as little
overestimation of the enclosed set as possible. For example, parallelepipeds
[8, 13, 17, 19], ellipsoids [11, 12, 25], convex polygons [29], and zonotopes [14]
have been employed to reduce the wrapping eect.
In [22], we linked the wrapping eect to the stability of an ITS method for
IVPs for ODEs and interpreted the problem of reducing the wrapping eect as
one of nding a more stable scheme for advancing the solution. This allowed
us to study the stability of several ITS methods (and thereby the wrapping
eect) by employing eigenvalue techniques, which have proven so useful in
studying the stability of point methods.
In this section, we explain rst how the wrapping eect can cause instability
in the direct method and then show that Lohner's QR-factorization method
provides a more stable scheme for propagating the error.
6.1 The Wrapping Eect as a Source of Instability
Consider the IVP
2.
applied with a constant stepsize h and order k to
(37), the direct method (22) reduces to
Taking widths on both sides of (38) and using (11) and (13), we obtain
w([y
We can interpret w([y j ]) as the size of the bound on the global error and
w([z j ]) as the size of the bound on the local error. We can keep w([z j ]) small
by decreasing the stepsize or increasing the order. However, T [y can be
a large overestimation of the set f Ty g. The reason for this
is that f Ty is not generally a box, but it is enclosed, or
wrapped, by the box T [y j 1 ]; cf. Figure 1. Note that we wish to compute
g, but this is generally infeasible. So we compute T [y
instead. That is, we normally compute on each step products of the form
consequently, we may incur a wrapping on each step, resulting
in unacceptably wide bounds for the solution of (37).
Consider now the point Taylor series (PTS) method given by
If we denote by - the global error of this method at t j , then
where
z j is the local truncation error of this method.
An important observation is that the global error in the PTS method propagates
with T , while the global error in the direct method propagates with
jT j. Denote by (A) the spectral radius of A 2 R nn . It is well-known that
Assuming that T and jT j can be diagonalized, we showed in
[22] that
then the bounds on the global error in these two methods
are almost the same, and the wrapping eect is not a serious di-culty for
the direct method; and
then the global error in the direct method can be much
larger than the global error in the PTS method, and the direct method may
suer signicantly from wrapping eect.
Thus, we can associate the wrapping eect in an ITS method with its global
error being signicantly larger than the global error of the corresponding PTS
method. In particular, if (T ) < 1 and (jT the PTS method is stable,
but the associated ITS method is likely asymptotically unstable in the sense
that lim j!1 kw([y j To reduce the wrapping eect, or improve the
stability of an ITS method, we must nd a more stable scheme for advancing
the solution.
6.2 How Lohner's QR-factorization Method Improves Stability
When applied with a constant stepsize and order to (37), an ITS method
incorporating Lohner's QR-factorization scheme, which together we refer to
as Lohner's method, can be written as
The transformation matrices A j 2 R nn satisfy
where A is orthogonal, and R j is upper triangular.
Using (42{43), we write (40{41) in the following equivalent form:
The interval vector [r j ] in (45) can be interpreted as an enclosure of the global
error (at that is propagated to the next step. Obviously, we are interested
in keeping the overestimation in [r j ] as small as possible.
is not much bigger than
w([z j ]), which we keep su-ciently small (by reducing the stepsize or increasing
the order). Since
we can consider jR as the matrix for propagating the global error in the QR
method. As we discuss below, the nature of the R j matrices ensures that this
method is normally more stable than the direct method.
A key observation here is that (42{43) is the simultaneous iteration (see for
example [32]) for computing the eigenvalues of T . This iteration is closely
related to Francis' QR algorithm [9] for nding the eigenvalues of T . To see
this, let
(j  1): (47)
Note that S
Now observe that the iteration (48{49) is just the unshifted QR algorithm for
nding the eigenvalues of T .
6.2.1 Eigenvalues of Dierent Magnitudes
Let T be nonsingular with eigenvalues f
of distinct mag-
nitudes. Then, from Theorem 3 in [9], under the QR iteration (48{49), the
elements below the principal diagonal of S j tend to zero, the moduli of those
above the diagonal tend to xed values, and the elements on the principal
diagonal tend to the eigenvalues of T .
Using this theorem, we showed in [22] that lim j!1 jR
matrix R with eigenvalues f j
g. Thus,
lim
Using (50), we derived in [22] an upper bound for the global error of Lohner's
method and showed that this bound is not much bigger than the bound for
the global error of the corresponding PTS method (39), for the same stepsize
and order. We also showed on several examples, that, at each t j , the global
error of Lohner's method is not much bigger than the global error of the PTS
method.
6.2.2 Eigenvalues of Equal Magnitude
If T is nonsingular with p eigenvalues of equal magnitude, then the matrices S j ,
dened in (47), tend to a block upper-triangular form with a pp block on its
main diagonal [10]. The eigenvalues of this block tend to these p eigenvalues.
If T is nonsingular with at least one complex conjugate pair of eigenvalues and
at most two eigenvalues of the same magnitude, then S j tends to a block upper-triangular
form with 1  1 and 2  2 blocks on its diagonal. The eigenvalues
of these blocks tend to the eigenvalues of T . Note that S j does not converge
to a xed matrix, but the eigenvalues of its diagonal blocks converge to the
eigenvalues of T .
For this case, we showed in [22] that, as j !1,
if T has a dominant complex conjugate pair of eigenvalues, then
and (R j ) oscillates; and
if T has a unique real eigenvalue of maximal modulus, then we should generally
expect that
Since the matrices jR j j do not approach a xed matrix, as in the case of
eigenvalues of distinct magnitudes, deriving a bound for the global error in
this case is more di-cult. However, we illustrated with four examples in [22]
that the global error in Lohner's method is not normally much bigger than
the corresponding global error of the PTS method.



--R

Introduction to Interval Computations.
FADBAD, a exible C
TADIFF, a exible C


Validating an a priori enclosure using high-order Taylor series

The Solution of Initial Value Problems Using Interval Arithmetic.
The QR transformation: A unitary analogue to the LR transformation
The Algebraic Eigenvalue Problem.
Automatic error analysis for the solution of ordinary di
A computable error bound for systems of ordinary di


Computational Methods in Ordinary Di
Enclosing the solutions of ordinary initial and boundary value problems.

Step size and order control in the veri
Interval Analysis.
Computing Rigorous Bounds on the Solution of an Initial Value Problem for an Ordinary Di
An interval Hermite-Obreschko method for computing rigorous bounds on the solution of an initial value problem for an ordinary dierential equation
A new perspective on the wrapping e
Validated solutions of initial value problems for ordinary di

Global, rigorous and realistic bounds for the solution of dissipative di
A First Course in Numerical Analysis.
On a class of enclosure methods for initial value problems.
On the existence and the veri
A heuristic to reduce the wrapping e
On higher order stable implicit methods for solving parabolic di
On the integration of sti
Understanding the QR algorithm.
--TR
The algebraic eigenvalue problem
Rigorously computed orbits of dynamical systems without the wrapping effect
Computing rigorous bounds on the solution of an initial value problem for an ordinary differential equation

--CTR
Hidde De Jong, Qualitative simulation and related approaches for the analysis of dynamic systems, The Knowledge Engineering Review, v.19 n.2, p.93-132, June 2004
wrapping effect;interval methods;validated methods;QR algorithm;ordinary differential equations;simultaneous iteration;taylor series;initial value problems;stability