Volume 4 No. 4 2003, 379-456
In tern ational Journal of Pure and A pplied Mathem atics
TAYLOR MODELS AND OTHER VA LIDATED
FUNCTIONA L INCLUSION METHODS
Kyoko Makino
1
,MartinBerz
2 §
1
Department of Physics
Universit y of Illinois at Urbana-Cham paign
1110 W. Green Street, Urbana, IL 61801-3080, U SA
2
Departmen t of Physics and A stronomy
Mich igan State Univ ersity
East Lansing, M I 48824, USA
e-mail:
Abstract: A detailed co mp arison bet ween Taylor model methods an d
other tools for validated com pu ta tion s is pro v ided . Basic elements of the
Tay lor model (TM ) methods are r ev iewed, beginning with the ar ithm etic
for elementary operations and in trinsic functions. We discuss some of
the fundamental properties, i ncluding high approximation order and the
abilit y t o control the dependency problem, a nd poin ters to man y of the
more advanced TM tools are pro vid ed. Aspects of the curren t imple-
mentation, and i n particular the issue of floating point e rror con trol, a re
discussed.
For t he purpose of p ro viding range enclosures, w e compare with
modern v ersions of centered form s and mean value f orms, as well as t he
direct co mputation o f remainder bounds by high-order interval auto-
matic differentiation and sho w th e adva ntages of the TM methods.
We also compare with the so-called boundary arithmetic (BA) of

terms and usually hav e a large number of local extrema; to mak e
matters worse, t hey exhibit a very significant cancellation problem. The
normal form defect functions themselv es are obtained from the high-
order dependence of solutions of ODEs on initial condition s. In various
meetings an d a large number of private discussions, the authors posed
this combined range bounding and integration problem to the interval
community as an interesting pr oject. How e ver, it w as uniform ly believe d
that because of dependency pr oblem in the nor m al for m defect functions,
the dimensionalit y, and the need to determine high-order d ependencies
on initial conditions in the OD E integration, the pro blem is intractable
through any of the tools known in the community. A nd indeed, the
attempt to apply var ious state of the a rt packages was not successful.
As a remedy to this problem, we developed the Ta ylor model ap-
proac h as an augmentation to earlier work on high-order m ultivariate
autom atic differentiation a n d the d ifferen tial algebraic methods to solv e
ODEs. Specifically, final variables in a code list are expressed in terms of
a high-order m ultivariate floating point Ta ylo r polyno m ia l of initial vari-
ables, plus a remainder bound a ccounting for the approximation error.
TAYLOR MODELS AND OTHER VALIDATED 381
Over suitably small d om ains, the polynomial represen tation is n atura lly
free of most of the dependency problem that the underlying function
ma y hav e had. At each node of the code list, the remainder bound is
calculated in parallel to the floating point coefficients; s in ce this only re-
quires information about the cu rrent Taylor coefficients, its calculation
itself is also free of much of t h e dependency problem of the original code
list; details will become clear in t he definition of the a rithm etic and in
the various examples that will be pro vided.
For the purpose of motivation, consider the problem of studying the
behavior of the polynomial function
f(x)=−371.9362500 − 791.2465656 · x +4044.944143 · x

17
+0.0002274682229 · x
18
(1.1)
in a validated w ay ove r a sufficiently small ra nge including the point x =
2. Becauseofthelargecoefficients and the alternating signs, a treatment
with in terval arithmetic, or mor e advan ced tools like cen tered forms,
will suffer from significant o verestimation because of th e c ancellation o f
terms. Ho wever, if before evaluation of th e function, the fu nction is first
re-expanded in po wers of (x − 2), it assu m es the following form
f(x)=−.1181179453 − 4.339394861 · (x −2) − 23.05727974 · (x − 2)
2
+14.04340823 · (x − 2)
3
+316.6727626 · (x − 2)
4
+583.1235424 · (x − 2)
5
− 157.0468495 · (x − 2)
6
− 1261.784612 · (x − 2)
7
− 858 .7604751 · (x − 2)
8
+271.5211596 · (x − 2)
9
+454.2310790 · (x − 2)
10
+107.4309653 · (x − 2)
11

(x − 2) is th e polynomial comprised of orders 0 th ro ugh
12 of f an d w e are in terested in studying o ver the domain [1.9, 2.1], then
o ver this d om a in w e can a sser t f(x) ∈ P
12
(x −2) + [−2 · 10
−12
, 2 ·10
−12
].
Ev en in this t runcated form, we can study m uch o f the beha vior of the
function; f or example, range bounding will only incur an additional over-
estimation of about 10
−12
, and integration can be done t o that a ccur acy
as well. S o we observ e that the simple trick of re-expanding around a
suitable point g reat ly simplified th e f un ctional behavior for the purpose
of using validated m eth ods.
Apparently the idea applies to any polynomial function, a lso in more
than one variables. It also easily gen eraliz es to ra tional fun ctions, since
these can be written as o rdered pairs (P, Q) of polynomials that can be
studied separately. The ordered pairs can be added and m ultiplied in
the ob v ious way.
The Taylor m odel m ethods in troduced in [ 112], [113] an d discussed
below capitalize on this observation b y represen ting an y functional de-
pendency in terms of a (Taylor) polynomial of sufficiently high order,
plus a small interval bound c ap turing the parts of the f unctio n th at d e-
viate from the polynomial. As such i t is m erely a v alidated extension of
autom atic differentiation methods[63], [ 2 0], n a mely tho se of high order
in many variables [11], [14], [61]; or in a mor e gen e ral context, the fact
known to scien tists of all back grou nd s that locally, smooth functions can

computer c ode list by a Taylor polynomial and a remainder bound
with a sharpness that scales with order (n +1) of t he width of the
domain.
2. The abilit y to alleviate the dependency problem in the calculation.
3. T h e ability to scale fav ora ble to hig h er dimensional problems.
We begin with a review of the d efinitions of the basic operations.
Definition 1. (Ta ylo r Model) Let f : D ⊂ R
v
→ R be a function
that is (n +1) tim es c o ntin u ously partially differen tiab le on an open set
384 K.Makino,M.Berz
containing the domain D. Let x
0
be a point in D and P the n-th order
Taylor polynomial of f around x
0
. Let I be an interval such that
f(x) ∈ P(x −x
0
)+I for all x ∈ D. (2.1)
Then we call the pair (P, I) an n-th order Taylor m odel of f around x
0
on D.
Apparently P + I encloses f between tw o hypersurfaces on D. As a
first step, we dev elop methods to calculate Taylor models from those of
smaller pieces.
Definition 2. (Ad dition and Multip lication of Taylo r Models)
Let T
1,2
=(P

1·2
is the part of the polynomi al P
1
· P
2
up to order n and
I
1·2
= B(P
e
)+B(P
1
) · I
2
+ B(P
2
) · I
1
+ I
1
· I
2
where P
e
is t he part of the polynomial P
1
·P
2
of orders (n+1) to 2n,and
B(P ) denotes a bound of P on th e domain D. We demand that B(P )

P (x −x
0
)=P (x −x
0
) −c
f
, so th at (
¯
P,I) is a
Taylor m odel for
¯
f. For the various intrinsics, we proceed as follow s.
TAYLOR MODELS AND OTHER VALIDATED 385
Exponential. We first write
exp(f(x)) = exp
¡
c
f
+
¯
f(x)
¢
=exp(c
f
) · exp
¡
¯
f(x)
¢
=exp(c

f(x)
¢
¾
, (2.2)
where 0 <θ<1.Takingk ≥ n,thepart
exp(c
f
) ·
½
1+
¯
f(x)+
1
2!
(
¯
f(x))
2
+ ···+
1
n!
(
¯
f(x))
n
¾
is merely a polynom ial of
¯
f, of which w e can obtain the Ta ylor m odel
via Taylor model addition and m ultiplication. The remainder part of

(k +1)-st powers of the Taylor model (
¯
P,I) of
¯
f will have vanishing
polynomial part, and thus so does the e n tire remainder part (2.3). The
remainder bound interval for the Lagrange remainder term
exp(c
f
)
1
(k +1)!
(
¯
f(x))
k+1
exp
¡
θ ·
¯
f(x)
¢
can be estim ated because, for a ny x ∈ D,
¯
P (x − x
0
) ∈ B(
¯
P ),and
0 <θ<1,andso

)+I) ⊂
(0, ∞), we first write as follow s
log(f(x)) = log c
f
+
¯
f(x)
c
f
−
1
2
(
¯
f(x))
2
c
2
f
+ ···+(−1)
k+1
1
k
(
¯
f(x))
k
c
k
f

0
)+I), wewriteasfollows:
1
f(x)
=
1
c
f
·
(
1 −
¯
f(x)
c
f
+
(
¯
f(x))
2
c
2
f
− ···+(−1)
k
(
¯
f(x))
k
c

c
f
·
(
1+
1
2
¯
f(x)
c
f
−
1
2!2
2
(
¯
f(x))
2
c
2
f
+ ···+(−1)
k−1
(2k − 3)!!
k!2
k
(
¯
f(x))

tribution from the rema inder term .
TAYLOR MODELS AND OTHER VALIDATED 387
M u ltip lic a tive inverse of squar e root. Under the condition ∀x ∈
D, B(P (x −x
0
)+I) ⊂ (0, ∞), we rew rite the expression
1
p
f(x)
=
1
√
c
f
·
(
1 −
1
2
¯
f(x)
c
f
+
3!!
2!2
2
(
¯
f(x))

k+1
c
k+1
f
1
¡
1+θ ·
¯
f(x)/c
f
¢
k+3/2
and evaluate in Taylor model arithmetic, o btaining a pure interval con-
tribution from the remain der term .
Sine. We use the addition theorem and power series expansion of
the sine function and obt ain
sin(f(x)) = sin(c
f
)+cos(c
f
) ·
¯
f(x) −
1
2!
sin(c
f
) · (
¯
f(x))

½
cos(c
f
+ θ ·
¯
f(x)) if k is even,
sin(c
f
+ θ ·
¯
f(x)) else,
and evaluate in Taylor model arithmetic; t he last term generates merely
an in terval contribution.
Cosine. Si m ilarly, we have
cos(f(x)) = cos(c
f
) −sin(c
f
) ·
¯
f(x) −
1
2!
cos(c
f
) · (
¯
f(x))
2
+

f
+ θ ·
¯
f(x)) if k is ev en ,
cos(c
f
+ θ ·
¯
f(x)) else.
388 K.Makino,M.Berz
Hy perbolic sine. In a similar v ein, we ha ve
sinh(f(x)) = sinh(c
f
)+cosh(c
f
) ·
¯
f(x)+
1
2!
sinh(c
f
) · (
¯
f(x))
2
+
1
3!
cosh(c

f
) ·
¯
f(x)+
1
2!
cosh(c
f
) · (
¯
f(x))
2
+
1
3!
sinh(c
f
) · (
¯
f(x))
3
+ ···+
1
(k +1)!
(
¯
f(x))
k+1
· J,
where

2
´
.
Utilizing th at
g(x) ≡ f(x) ·
q
1 −c
2
f
− c
f
·
p
1 −(f(x))
2
does not h av e a constan t part, w e ha ve
arcsin(g(x)) = g(x)+
1
3!
(g(x))
3
+
3
2
5!
(g(x))
5
+
3
2

(a)=(1+2a
2
)/(1 − a
2
)
5/2
,
A recursive form u la for the higher order derivatives of arcsin
arcsin
(k+2)
(a)=
1
1 −a
2
{(2k +1)a arcsin
(k+1)
(a)+k
2
arcsin
(k)
(a)}
is useful [13 2]. Then, eva lu at in g i n Taylor model arithmetic yields the
desired r esult, where aga in the terms involving θ only produce in terval
contributions.
Arccosine. Use arccos(f(x)) = π/2 −arcsin(f(x)).
Arctangen t. Using an addition form u la for the arctangent, w e ha ve
arctan(f(x)) = arctan(c
f
) + arctan
µ

5
−
1
7
(g(x))
7
+ ···+
1
k +1
(g(x))
k+1
· cos
k+1
(arctan(θ · g(x)) · sin
³
(k +1)·
³
arctan(θ · g(x)) +
π
2
´´
and proceed as usual.
An tideriv ation. We note that a Taylor model for the integ ral with
respect to variable i of a function f ca n be obtained from the Taylor
model (P, I) of t he function b y merely integrating t he part P
n−1
of order
up to n −1 of the polynomial, and bounding the n-th order into the new
remainder bound. Specifically, w e have
∂

,I
f,h
) and (P
g,h
,I
g,h
) be n-th order Taylor mode ls for f and
g around x
h
on x
h
+[−h, h]
v
⊂ D. Let the remainder bounds I
f,h
and
I
g,h
satisfy I
f,h
=O(h
n+1
) and I
g,h
=O(h
n+1
). Then the Taylor mod els
(P
f+g
,I

Pro of. The proof for the binary o perations follo ws directly from the
definition of the remainder bounds for the binaries. Similarly, the p roof
for the intrinsics follow s because all in trin sics are composed of binary
operations as well as an additional i nterval, the width o f which scales at
least with th e (n +1)-st power of a bound B of a function that scales
at least linearly with h. ¤
Remark 1. (High O rder S calin g P roperty) T he h igh ord er scaling
property of Taylor model arithmetic states that a giv en function f can
be appro xim ated by another function P (apolynomial)withanerror
that scales with high order as the domain decreases. This approximation
statement follows sta n dard mathem atica l p r actice. How ever, in the in-
terval comm unity it is customary to study another related but different
meanin g of scaling: namely the behavior of the o verestima tion of a given
method to determine the range of a function. In the con ventional in-
terval comm un ity, this scalin g propert y is importan t because in tervals,
including range in tervals, pla y a leading role. In the w orld of Taylor
model a lg orithm s, the u se of intervals them selves is m uch reduced, s ince
as a gener al rule, expression s are kept in Taylor model form as muc h a s
possible, for e xa m ple to retain the a b ility to suppress dependency. Thus
in general, the high order scaling propert y as s tated in t h e previous the-
orem is the r elevant one. This, ho wever, a pp lies only in a lim ited sense
to the question of range bounding; more about this matter below an d
in [120].
TAYLOR MODELS AND OTHER VALIDATED 391
Having defined the in trinsics of Taylor model arithmetic as abo v e,
we can summarize the main p ropert y of Taylor model arithme tic in the
following th e o r em:
Theorem 2. ( F TT MA , Fundamental Theorem of Taylor Model
Arithme tic) Let the f un ctio n f : R
v

sho w some resu lts of a computation of the function sin
2
(exp(x +1))+
cos
2
(exp(x+1)), executed with an impleme ntation of Taylor model arith-
metic as discussed in the next section. Of course the function is identica l
to 1, but the validated methods cannot capitalize on this information;
so this fu nction can serv e as a good example t o a ssess th e t ig htness of
various enclosure schemes. The left picture in Fig ure 1 sho w s th e result
of the enclosure of the function by inter va ls, m ean value form, cen ter ed
form, and t he r esult of the Taylor model range bounding algorithm for
the dom ain s [−2
−j
, 2
−j
] for j =1, , 7; more comparison s about these
met hods a n d Taylor models follow below. Also shown in the rig ht pic -
ture are empirically compu ted approxim ation orders as a function of j.
Indeed it ca n be seen th at the width of the com pu ted h igher order re-
maind er intervals scale with o rd er (n +1)forTaylormodelsofordern,
un til near the floor of m ac hine precision, a t which point rounding e ffects
dominate.
Asasidenotewealsoobservethatintherepresentationofafunction
through its Taylor model, it is apparent t hat some functions that can be
represented exactly by intervals cannot be represented exactly by Tay-
lor models; a situation that also occurs with other advanced inclusion
tools like cen tered forms. As a n exam ple of this effect, we consider t he
func tion f(x)=1/x. Figure 2 shows the beha vior of the TM method
392 K.Makino,M.Berz

6
9
9
9
9
9
9
9
12
12
12
12
12
12
12
j
l
o
g
1
0
q
INTERVAL
CENTERED
MEAN VALUE
1ST ORDER TM
3RD ORDER TM
6TH ORDER TM
9TH ORDER TM
12TH ORDER TM

3
3
3
6
6
6
6
6
6
9
9
9
9
12
12
j
E
A
O
INTERVAL
CENTERED
MEAN VALUE
1ST ORDER TM
3RD ORDER TM
6TH ORDER TM
9TH ORDER TM
12TH ORDER TM
1
3
6

Since in the Taylor model approach, the coefficients are floating point
(FP) numbers, care must be taken t hat t he i naccuracies of conven tional
FP ar ithm etic are properly accou nted for. Algorithmically the methods
are rather straigh tforward; however for practical use of the meth ods,
the more importan t question i s that o f the sound n ess of the actual im-
plementation. Besides the tests performed in the develop m ent of the
TAYLOR MODELS AND OTHER VALIDATED 393
1 2 3 4 5 6 7
-16
-14
-12
-10
-8
-6
-4
-2
0
1
1
1
1
1
1
1
3
3
3
3
3
3

CENTERED
MEAN VALUE
1ST ORDER TM LDB
3RD ORDER TM LDB
5TH ORDER TM LDB
7TH ORDER TM LDB
9TH ORDER TM LDB
1
3
5
7
9
1 2 3 4 5 6
0
1
2
3
4
5
6
7
8
9
10
11
12
13
1
1
1

3
5
7
9
Figure 2: Relativ e overestim ation q (left) and em pirica l app roxima -
tion order (righ t) for the function 1/x with LDB range bounder in
2+[−2
−j
, 2
−j
].
program , various other tests have been performed. Corliss and Yu per-
formed extensive t ests of the COSY in terval tools by porting of COSY
interval results to Maple in binary format and comparison with Maple
computations with nearly 1000 digits of accuracy. Sev era l thousand
cases that are to be considered particularly difficult as w ell as around
10
6
random tests spanning all orders of magnitude of allowed domains
of the in trin sics were performed[36]. Independently, Rev ol performed
around 10
8
random tests of the in terval arith m e tic b y com pariso n with
a guaranteed precision library for elementary operations and intrinsic
functions[156]. I n addition, Revol prove d the soundness of the algo-
rithms in the floating point coefficient treatment of the Taylor model
implementation and chec k ed the actual coding [157].
Definition 4. (Adm issib le FP Arithmetic) We assume com pu ta-
tion is performed in a floa ting point environment supporting the four
elementary operations ⊕, ⊗, Ä, ®. We call the arithmetic admissible if

,b
2
] of floating point num bers and any
° ∈ {⊕, ⊗, Ä, ®} and c orresponding real operation ◦∈{+, ×, −,/},
we have
a[a
1
,b
1
] ° [a
2
,b
2
] ⊃ {x ◦ y|x ∈ [a
1
,b
1
],y∈ [a
2
,b
2
]}, (3.1)
and fu rthermore, for any interval intrinsic s ∈ S representin g the real
func tion s, we have
s([a, b]) ⊃ {s(x)|x ∈ [a, b]}. (3.2)
For the specific purposes of Taylor model arithmetic, some addi-
tional considerations are necessary. First we note that combinatorial
argumen ts sho w [17] that the number of nonzero c oefficients in a poly-
nomial of order n in v va riables c an not exceed(n+v)!/(n! · v!). Further-
more, as a lso sho w n in [17], the n u mber of multiplications necessary to

threshold ε
c
can be ch osen over a wide possible range, but since it is later
used to control the number of coefficients actively retained in the Tay lor
model arithmetic, a value n ot too far below ε
m
,suchasε
c
=10
−20
, is
a good cho ice for man y cases. Furthermore, for essentially all practically
conceivab le ca ses o f n and v,thechoicee =2is satisfactory, and this is
thenumberusedinourimplementation.
Under t he assum p tion o f the above properties of the floating poin t
arithmetic, inter va l arithm etic, and the Ta ylor model arithmetic con-
stan ts, w e no w describe th e algorithms for Taylor model arith metic,
whichwillleadtothedefinition of admissible FP Taylor model arith-
metic.
Storage. In the COSY implementation, a Taylor model T of order
n and dim ension v is represented b y a collection of nonzero floating point
coefficients a
i
, as w ell a s two coding integers n
i,1
and n
i,2
that contain
uniqu e information allo w ing to iden tify the term to which the coefficient
a

| >ε
c
, are retained. In man y p ractical ca ses, this en-
tails significant savings in space and execution time; more on how the
non-retained terms are treated is described belo w . Since b y require-
ment, ε
2
c
>ε
u
, th e m ultiplica tion of two retained coefficien ts can never
lead to underflo w . B esides the coefficients and coding integ ers, eac h
TM also conta ins an interval I composed of tw o floating poin t n umbers
representing rigorous enclosures of the remainder bound.
Err o r collection. In the elementary operations of Taylor models,
396 K.Makino,M.Berz
the errors due to floa ting poin t arithm etic are accum u lated in a float-
ing point “tallyin g variable” t whichintheendisusedtoincreasethe
remainder bound in terval I by an interva l of the form e ⊗ ε
m
⊗ [−t, t].
The factor e assures a safe upper bound of all floating poin t errors of
adding up the (positive) con tributions to t. Accounting for t he error
through a single floating point v ariable t with the factor e· ε
m
“factored
out” notably increases computational efficiency. In a ddition , there is a
“sweeping variable” s that w ill be used to ab sorb terms that fall below
the cutoff threshold ε
c

ing polynom ial, and k will be incremen ted . If |b
k
| <ε
c
, the sw eeping
variable s is increment ed by |b
k
|. After all terms have been treated, the
total remainder bound o f t he result of the s calar m ultiplication is set to
be [c, c] ⊗I ⊕e ⊗ε
m
⊗[−t, t] ⊕e ⊗[−s, s], which i s perform ed in out ward
rounded in terval arithmetic.
Addition. Add ition of two Taylor m odels T
(1)
and T
(2)
with coef-
ficients a
(1)
i
and a
(2)
j
, coding integers (n
(1)
i,1
,n
(1)
i,2

(1)
i,2
) 6=(n
(2)
j,1
,n
(2)
j,2
), the t erm
that should come first according to the ordering is merely copied, and
its pointer as well as k ar e incremen ted. In c ase (n
(1)
i,1
,n
(1)
i,2
)=(n
(2)
j,1
,n
(2)
j,2
),
w e proceed as follo ws. We determ ine the floating poin t coefficient b
k
=
a
(1)
i
⊕a

⊕ e ⊗ε
m
⊗ [−t, t] ⊕ e ⊗ [−s, s], which
TAYLOR MODELS AND OTHER VALIDATED 397
is performed in outward rounded in terval a rithm etic.
M u ltiplication . T he multiplication of t wo Ta y lor models T
(1)
and
T
(2)
of order n with coefficien ts a
(1)
i
and a
(2)
j
and coding in tegers (n
(1)
i,1
,
n
(1)
i,2
) and (n
(2)
j,1
,n
(2)
j,2
), r espectively, is performed as follo ws. T he con-

(1)
i
⊗ a
(2)
j
of the coefficients. To account for the error, w e incre-
ment t by |p|. We add the term p to the coefficient b
l
. To acco unt for
the error, we incremen t t by max(|p|, |b
l
|).
After all monomial multiplications have been executed, all result-
ing total coefficients b
l
of the product polynomial will be studied for
sw eeping. If |b
l
| ≥ ε
c
, the term will be included in the resulting poly-
nomial, and l will be increm ented. If |b
l
| <ε
c
, the sweeping variable
s is incremented by |b
l
|, but l will not be incremented, i.e. the term
is not retained. In the end, the remainder bound I is increm ented by

Definition 7. (Admissible F P Taylor Model Ar ithm etic) We call
a Ta ylor model arithmetic ad m issible if it is ba sed o n an admissible FP
and interva l arithm etic and it adheres to th e algorithms for storage,
scalar multiplication, ad dition , mu ltiplication , an d intrinsic functions
described above.
Remark 3. (FP Ta ylor Model Arithmetic in COSY INFINITY )
The code COSY INFINITY con tains an admissible Taylor model arith-
metic in arbitrary order and in arbitrarily many variables. The code
consists of around 50, 000 lines of FOR TRAN’ 77 source that also cross-
compiles to s tand ard C. I t can be used in the env iron m ent of the COSY
language,aswellasinF77 and C. It is a lso a vailable as classes in F90
and C++. T h e code is highly optimized for performance in that any
o verhead for addressing of polynom ial coefficients amounts t o less than
30 percen t o f t he floating poin t arithmetic n ecessary for the coefficient
arithme tic [11]. It also h as full spar sity support i n that coefficients be-
low the cutoff threshold do not con trib ute to execution time and storage.
Remark 4. (Verificatio n and Validation of the COSY FP Taylor
Model Arithmetic) The FP TM arithmetic implem e nted in COSY is
currently being verified and validated b y two outside groups [36], [156]
with a s uite of c hallenging test problems. Independen t ly, the va lidity of
thealgorithmsformingthecoreoftheCOSY Tay lor model FP algorithm
have been verified by Rev ol [157].
4. Taylor Model Algorithms
The above algorithms for Taylor model arithmetic assure that also in
a computer en v ironment subject to floating point errors, any computa-
tions using Taylor models lead to rigorous enclosures, an d we obtain t h e
following result.
TAYLOR MODELS AND OTHER VALIDATED 399
Theorem 3. (Taylor M od el En closure T h eorem ) Let the function
f : R

properties of the remainder bound in a rigo rou s sense is lost. However,
these pr operties of Ta ylo r models are r etained in an app roximate fash ion .
Remark 5. (Influence of Floating Point Ar ithmetic) In th e p r es-
ence of floating point errors, the polynomial P will be a floating point
appro ximation of the Taylor polynomial of g ◦ f if P
f
w as an approxi-
mate Taylor polynomial for f. Fu rth er mor e, any (n +1)-st o rder scaling
property for the remainder interval will prevail approxim ately until near
the floor of machine precision.
As an imm ed iate conseq u ence, we obtain the following:
Algorithm 1. (Range Bounding with Taylor M odels)
Input: a finite code list in volving elementary operations and intrin-
sics describing the function f ov er the multivariate dom ain box D
Output: an enclosure of f in a Taylor model P
f
+I
f
, and an interval
bound B(f) for the ran ge of f over D
1. Set up a Taylor model T
I
enclosing the identity function. This is
comp rised of the linear m ultiva ria te polynomial P (x)=x plus the
remainder bound [0, 0].
2. Evaluate the code list for f in Taylor model arith m etic. As a
result, obtain P
f
+ I
f

f the sharpn ess of whic h scales with order
(n +1)with D
1. Set up a Ta ylor model T
I
enclosing the identity function. T his is
comp rised of the linear m u ltiva ria te polynomial P (x)=x plus the
remainder bound [0, 0].
2. Evalu ate the code list for f in Taylor model arithmetic. As a
result, obtain P + I.
3. Integrate the polynomial b y manipulation of coefficients to obtain
aprimitiveP
I
for P, and insert the endpoints of D into P
I
to
obtain the in teg ral
R
D
P.
4. Obtain an enclosu re for
R
D
f as
R
D
f ⊂
R
D
P + |D|·I
Various applications of the method are described in detail in [25].

It has recently been suggested that it would be useful to have a de-
tailed comparison between Taylor models and the cen tered form (CF)
and mean value form (MF) [127], [100] , [155], [98], [2], [1], [131] for range
boundin g. Since t he latter two usually pro v id e sharper enclosures than
intervals and earlier com parison s of Taylor models were m o stly with in-
tervals, it was suspected that for mere range bounding, t h e performan ce
of Ta ylor models would be rather similar to CF a nd MF, whic h are
kno w n to h ave t h e quadratic approximation property. In this section w e
attempt a comparison based on what w e believe to be a limited collec-
tion of meaningful examples. We comp are with Ta ylor model m ethods
of various orders, and sub sequent bounding sch em es based on either
naive interva l evalua tion of the Taylor polynom ial, or based on th e lin-
ear dominated bounder LDB [120]. To increase the demand on the LDB
method, in all examples sho w n no domain subdivisions as utilized in the
vario us Bernste in-b ase d sc h em e s [133], [134] are allo wed. Ap p are ntly
allowing subdivision before applying LDB wo u ld increase the applica-
bilit y of LDB to larger doma ins. We observe that overall, Tay lor m odels
suppress dependency mu ch better t han centered forms and mean value
402 K.Makino,M.Berz
forms, resulting in frequently muc h sharper inclusions. Furthermor e, in
many cases the LDB method leads to higher order enclosures of esti-
mated ranges.
All com putations a re performed using COSY for the Taylor models,
while interva ls, centered f orm s, a nd slopes were eva luated using the i m -
plementation in the INTLA B toolbo x for Matlab [165]. Specifically, we
used IN T LAB Version 3.1 under Matlab Version 6 . We believ e we ha ve
used the code in the proper way, although documentation is somewhat
terse; as the author puts it, “To be frankly, there is not m uch other
documentation about INTLA B . In every routine, of course, the func-
tionalit y is documented. Otherwise, w e think INTLA B code is m uch

the denom inator th en scales with the second (or a higher) po wer of the
domain wi dth. We usually list the E AO only until the floor of m ac hine
precision is reac h ed. We frequently also list the aver age empiric al ap-
proximation order (AEAO) for various methods, which is obtain ed by
TAYLOR MODELS AND OTHER VALIDATED 403
a v eraging the EAO data for the given method o v er all choices of the
domain width.
For n otation a l simplicit y, in the fo llowing pictures, resu lts obtained
using inter va l evaluation will be denoted by the sym bol ¡, remin iscent
of an interval bo x, while those obtained by the mean value form and
cen tered form will be denoted by the sym bols ∇ an d 4, reminiscent of
a gradient and a differ en ce quotient, respectiv e ly. Taylor models will be
identified b y numbers corresponding to their orders.
We begin our discussion with the study of a simple three dimensional
examp le function w ith m odest dependency but ov erall rather in nocen t
behavior studied in [112]. The function has the form
f
1
(x, y, z)=
4tan(3y)
3x + x
q
6x
−7(x−8)
− 120 − 2x − 7z(1 + 2y)
− sinh
µ
0.5+
6y
8y +7

asMF,asisfrequentlyobserved. Thefirst o rder Taylor model m ethod
behaves sim ilar to CF , and is in fact slightly superior. The h igh er order
Ta ylor models, while still sho w in g order 2 scaling, prov id e enclosures
that is about 1 order of m agnitude sharper t han those of CF.
With LDB, the approximation order of the Taylor model o f order n
increases t o (n +1), until the floor of mac hine precision is reached. At
the most favorable point, the sharpness o f the 9-th o rd er Taylo r model
method is about 11 orde rs of m agnitude higher th an that of CF.