Computing with Intervals

The ultimate symbolic computation?

Yves A. Papegay and David Daney

INRIA Sophia Antipolis, COPRIN
(contact Yves.Papegay@sophia.inria.fr)

$NumberMarks=True;
<<UnCertainties`;

Abstract

Two years ago, we presented at the Wolfram Technology Conference 2005, the ALIASTICA package, interfacing Mathematica with ALIAS, an extensive C++ library of interval analysis algorithms. This interface is bringing a flavour of symbolic computation skills on top of the numerical library. But we are now convinced, that symbolic computation should interact with interval computation at every points of the interval analysis algorithms to make them even more efficients. That leads us to start the development of an interval analysis package completely implemented in Mathematica.

Interval Arithmetic and Interval Analysis

It is important to distinguish the very old concept of intervals, from the Interval Arithmetic - formally defined in the 60's, and from the scientific activity of Interval Analysis with its own community and its own journal: Reliable Computing.

Interval Arithmetic and Interval Analysis

Interval Arithmetic is (well) implemented in Mathematica. A few features of Interval Analysis are in the NumericalMath`IntervalRoots legacy package that is now available on MathSource at
http://library.wolfram.com/infocenter/MathSource/6779

Motivations for using intervals are the certification of the results of computations and the ability of considering set-valued parameters.

Perform Certified Computations

The Rump function is just another example of numerical instability.

RumpFunc[x_,y_]:=(1335/4-x^2)y^6+x^2(11 x^2 y^2-121 y^4-2)+(11/2) y^8+x/(2y)

RumpFunc[77617,33096]//N

RumpFunc[77617.,33096.]

RumpFuncN[x_,y_]:=(1335/4-x^2)y^6+x^2(11 x^2 y^2-121 y^4-2)+(5.5) y^8+x/(2y)

RumpFuncN[77617,33096]

RumpFuncN[77617.,33096.]

Overestimation

Interval Arithmetic certifies results, but only produces poor results, overestimated.

RumpFuncN[Interval[77617.],Interval[33096.]]

This overestimation is intrinsic : it comes from decorrelation of variables and the so-called wrapping effect.

Decorrelation

x=Interval[{-1.,2.}];

x*x
x-x

Overestimation

Wrapping Effect

Bisection and Filtering

Bisection and filtering are the numerical approaches to reduce the intrinsic overestimation.

Bisection

Bisection is used at evaluation time.

x1=Interval[{-1.,0.}];x2=Interval[{0.,2.}];x=IntervalUnion[x1,x2]

x-x
IntervalUnion[x1-x1,x2-x2]

x*x
IntervalUnion[x1*x1,x2*x2]

Bisection and Filtering

Bisection and filtering are the numerical approaches to reduce the intrinsic overestimation.

Filtering

Filtering is a way to cope with overestimation when solving (eg. Newton process)

NewtonFilter[Sin[x]+x-3,{x,Interval[{1.,5.}]}]

Symbolic Preprocessing

Symbolic manipulation allows to minimize overestimation when evaluating and could be done completely off-line.

Manipulation of Expressions

p[x_]:=x^2+x-6

p[Interval[{-1.,1.}]]

Factor[p[t]]
%/.t->Interval[{-1.,1.}]

HornerForm[p[t]]
%/.t->Interval[{-1.,1.}]

TrueQ[Expand[(t+1/2)^2-25/4]==p[t]]

(t+1/2)^2-25/4/.t->Interval[{-1.,1.}]//N

Symbolic Preprocessing

Manipulation of Expressions

Number of Occurences of Variables

x=.

p1=-6x^3+4x^3z^3+15x y^3z^3-3y^5 z^3-12x^2y z^3-3x y z^3+y^3 z^3;

p1/.{x->Interval[{-1,2}],y->Interval[{-2,1}],z->Interval[{-1.5,1.5}]}

HornerForm[p1]
%/.{x->Interval[{-1,2}],y->Interval[{-2,1}],z->Interval[{-1.5,1.5}]}

NbOcc[x_]:=Length[Cases[x,_Symbol,Infinity]]
OccSimplify[expr_]:=FullSimplify[expr,ComplexityFunction->NbOcc[#]&]

OccSimplify[p1]
OccSimplify[p1]/.{x->Interval[{-1,2}],y->Interval[{-2,1}],z->Interval[{-1.5,1.5}]}

Symbolic Preprocessing

Analysis of expressions during the evaluation also recquires symbolic computations that can be done off-line.

Analysis of Expressions

To be done during evaluation.

p2=-6x^3-4x^3z^3+15x y^3z^3+3y^5 z^3+12x^2y z^3+3x y z^3+y^3 z^3;

OccSimplify[p2]/.{x->Interval[{0.6,1}],y->Interval[{0.6,1.}],z->Interval[{0.6,1.}]}

D[p2,y]/.{x->Interval[{0.5,1}],y->Interval[{0.5,1}],z->Interval[{0.5,1}]}

p2m=p2/.y->0.6
p2M=p2/.y->1.

IntervalHull[OccSimplify[p2m]/.{x->Interval[{0.6,1}],z->Interval[{0.6,1.}]},
OccSimplify[p2M]/.{x->Interval[{0.6,1}],z->Interval[{0.6,1.}]}]

D[p2M,z]/.{x->Interval[{0.6,1}],z->Interval[{0.6,1.}]}

IntervalHull[OccSimplify[p2m]/.{x->Interval[{0.6,1}],z->Interval[{0.6,1.}]},OccSimplify[p2M/.z->0.6]/.x->Interval[{0.6,1}],
OccSimplify[p2M/.z->1]/.x->Interval[{0.6,1}]]

Symbolic Filtering

Symbolic filtering is an adaptation of Constraints Programming methods.

2-B Filter

f=Sin[x]+x-1;

f1=Replace[ f,x->y,{1}]

sol=Solve[f1==0,y]

filter[i_]:=N[y/.sol[[1]]/.x->i]

filter[Interval[{-10,10}]]

NestList[filter,Interval[{-10.,10.}],8]

NestWhile[filter,Interval[{-10.,10.}],Width[#]>0.00000001&]

Symbolic Execution of Algorithms

Newton Exemple

f=.

f[x_]:=Sin[x]+x-3

m=MidPoint[Interval[{1.,5.}]]

IntervalIntersection[Interval[{1.,5.}],m-f[m]/f'[Interval[{1.,5.}]]]

Symbolic Execution of Algorithms

Newton Exemple

IntervalIntersection[Interval[{1.,5.}],m-f[m]/f'[Interval[{1.,5.}]]]

2D exemple

g1[{x_,y_}]:=x^2+y^2-1
g2[{x_,y_}]:=(x-1)^2+y^2-1
m={0.1,0.1};

D[{g1[{x,y}],g2[{x,y}]},{{x,y}}]
%/.{x->Interval[{0.,1.}],y->Interval[{0.,1.}]}

m-IntervalLinearSolve[(({{Interval[{0.,2.}], Interval[{0.,2.}]}, {Interval[{-2.,0.}], Interval[{0.,2.}]}})),(({{-0.98}, {-0.18}}))]

m-IntervalLinearSolve[(({{2 x, 2 y}, {2 (-1+x), 2 y}})),(({{-0.98}, {-0.18}})),{x,Interval[{0.,1.}]},{y,Interval[{0.,1.}]}]

Map[IntervalIntersection[Interval[{0.05,0.9}],#]&,m-LinearSolve[(({{2 x, 2 y}, {2 (-1+x), 2 y}})),(({{-0.98}, {-0.18}}))]/.{x->Interval[{0.05,0.9}],y->Interval[{0.05,0.9}]}]//MatrixForm