intpakX.html

intpakX v1.0

written by Grimmer, Markus, Department of Mathematics, University of Wuppertal, Germany, http://www.math.uni-wuppertal.de/wrswt
<© 1999-2002 Scientific Computing/Software Engineering Research Group, University of Wuppertal, Germany>
intpakX
by Krämer, Walter, Geulig, Ilse and Grimmer, Markus, Department of Mathematics, University of Wuppertal, Germany, http://www.math.uni-wuppertal.de/wrswt
<© 1999-2002 Scientific Computing/Software Engineering Research Group, University of Wuppertal, Germany>
intpak
by Corless, R. and Connell, A., University of Western Ontario, Canada
<© 1992-1993 R. Corless and A. Connell>

> restart;

###########################################################################

intpakX v1.0 - Interval Arithmetic in Maple

This is a new, integrated version of the Maple package intpakX for interval computations.
It integrates and updates the package

- intpak by R. Corless and A. Connell
and its extension
- intpakX by I. Geulig and W. Kraemer

intpakX v1.0 was written/updated by M. Grimmer.

intpakX v1.0 provides basic data types and operations for interval arithmetic as well as additional
features for further interval computation.
Through this package, the features of the above mentioned packages are available in a
combined version as a Maple 7 module.
A couple of updates have been made to provide a version working with Maple 7 properly.

For information how to create a Maple package from this document, see the end of this worksheet.

intpakX v1.0 (c) Scientific Computing/Software Engineering Research Group,
Department of Mathematics, University of Wuppertal, Germany

Contact:
kraemer@math.uni-wuppertal.de
markus.grimmer@math.uni-wuppertal.de
http://www.math.uni-wuppertal.de/wrswt

Disclaimer: see end of module

###########################################################################

###########################################################################
Module header
###########################################################################

>    intpakX:=module()
export init,Evalf,min,intpakX_min,max,intpakX_max,ilog10,intpakX_ilog10,
  Interval_ulp,ulp,
  Interval_Round_Up,Interval_Round_Down,ru,rd,
  intpakX_greater,is_in,construct,expinfinity,
  Interval_exp,infinityln,`&exp`,`&ln`,
  Interval_ln,addinfinity,Interval_add,`&+`,subtractinfinity,
  Interval_subtract,`&-`,timesinfinity,Interval_times,`&*`,
  Interval_reciprocal,inv,Interval_divide,`&/`,sqrtinfinity,Interval_sqrt,`&sqrt`,
  sqrinfinity,Interval_sqr,`&sqr`,
  Interval_option_zero,powerinfinity,Interval_Integerpower,`&intpower`,
  Interval_power,`&**`,Interval_midpoint,midpoint,
  Interval_width,width,Interval_intersect,`&intersect`,
  Interval_union,`&union`,Interval_trig_ru,Interval_trig_rd,
  Interval_scale,Interval_range_values,Interval_sin,Interval_cos,Interval_tan,
  Interval_arcsin,Interval_arccos,Interval_arctan,Interval_hyp_rd,Interval_hyp_ru,
  Interval_cosh,Interval_sinh,Interval_tanh,`&sin`,`&cos`,`&tan`,`&arcsin`,
  `&arccos`,`&arctan`,`&sinh`,`&cosh`,`&tanh`,
  ext_int_div,rel_diam,mid,compute_naive_interval_range,compute_mean_value_range,
  compute_monotonic_range,compute_combined_range,compute_taylor_form_range,`&Convex_Hull`,
  subdivide_equidistant,subdivide_adaptive,interval_list_plot,compute_range,
  max_abs_error,cni_range3d,
  subdivide_equi3d,interval_list_plot3d,compute_range3d,newton,compute_all_zeros,newton_plot,
  newton_with_plot,compute_all_zeros_with_plot,
  `&cabs`,complex_disc_plot,complex_polynom_plot,`&cadd`,`&csub`,`&cmult`,`&cdiv`,x0_start,
  `&cmult_opt`,`&cdiv_opt`,horner_eval_cent,horner_eval_opt,centred_form_eval,cexp;
option package;

###########################################################################

Name definitions (see below for body of corresponding procedures)

###########################################################################

> min:=intpakX_min;
max:=intpakX_max;
ilog10:=intpakX_ilog10;
width:=Interval_width;
ulp:=Interval_ulp;
ru:=Interval_Round_Up;
rd:=Interval_Round_Down;
`&exp`:=Interval_exp;
`&ln`:=Interval_ln;
`&+`:=Interval_add;
`&-`:=Interval_subtract;
`&*`:=Interval_times;
inv:=Interval_reciprocal;
`&/`:=Interval_divide;
`&sqrt`:=Interval_sqrt;
`&sqr`:=Interval_sqr;
`&intpower`:=Interval_Integerpower;
`&**`:=Interval_power;
midpoint:=Interval_midpoint;
`&intersect`:=Interval_intersect;
`&union`:=Interval_union;
`&sin`:=Interval_sin;
`&cos`:=Interval_cos;
`&tan`:=Interval_tan;
`&arcsin`:=Interval_arcsin;
`&arccos`:=Interval_arccos;
`&arctan`:=Interval_arctan;
`&cosh`:=Interval_cosh;
`&sinh`:=Interval_sinh;
`&tanh`:=Interval_tanh;

###########################################################################

Initialization

The init procedure contains type definitions for the various (procedural) interval types.

In Maple, type defintions must be global; they seemingly cannot be module exports.

###########################################################################

> init:=proc()

> global Interval_fnlist,`type/num_or_FAIL`,`type/interval_comp`,`type/interval`,
`type/complex_disc`,`type/complex_interval`,`convert/interval`,inapply;

>

> Interval_fnlist:=[sqrt=`&sqrt`,ln=`&ln`,exp=`&exp`,sin=`&sin`,cos=`&cos`,

> tan=`&tan`,arcsin=`&arcsin`,arccos=`&arccos`,arctan=`&arctan`,sinh=`&sinh`,

> cosh=`&cosh`];

- - in proc init - - - - - - - - - - - - - - -

This type is included for ease of writing and clarity of code.

> `type/num_or_FAIL`:=proc(a)

> local bool,Constants:

> Constants:={constants}:

> bool:=type(a,numeric) or a=-infinity or a=infinity or a=FAIL

> or member(a,Constants):

> bool:

> end:

- - in proc init - - - - - - - - - - - - - - -

Simplified version of type interval_comp (cf. intpak/intpakX).
Pi, gamma, Catalan, false, true are not allowed as interval borders.

> `type/interval_comp`:=proc(x)

> local bool;

> bool:=type(x,float) or x=FAIL or x=-infinity or x=infinity or x=0;

> bool;

> end:

- - in proc init - - - - - - - - - - - - - - -
PROCEDURE - TYPE/INTERVAL

This procedure defines the procedural type interval.
It tests an argument to see if it is an interval. An interval

is defined here to be a list with either zero elements,+/- infinity,FAIL

or a list with two floating point members.

`type/interval`:=proc(x)

> local bool;

> bool:=false;

> if type(x,list) then

> if nops(x) = 2 then

> if type(x[1],interval_comp) and type(x[2],interval_comp) then

> if intpakX_greater(x[2],x[1]) then bool:=true

> elif x[1] = FAIL or x[2] = FAIL then bool:=true

> else ERROR(`enter the lowest endpoint first`)

> fi

> elif nops(x)=0 then bool:=true

> fi;

> bool;

> end:

   - - in proc init - - - - - - - - - - - - - - -
   PROCEDURE - TYPE/COMPLEX_DISC

   This procedure implements (i.e. tests an argument to be of) the procedural type complex disc
   for use in complex interval arithmetic.

> `type/complex_disc`:=proc(z)

> local bool;

> bool:=false;

> if type(z,list) then

> if nops(z)=3 then

> if type(z[3],numeric) then

> if z[3] >= 0 then

> bool:= type(z[1],numeric) and

> type(z[2],numeric);

> fi;

> bool;

> end:

   - - in proc init - - - - - - - - - - - - - - -
   PROCEDURE - TYPE/COMPLEX_INTERVAL

   This procedure implements (i.e. tests an argument to be of) the procedural type complex interval
   (~rectangular interval)
   for use in complex interval arithmetic.

> `type/complex_interval`:=proc(z)

> local bool;

> bool:=false;

> if type(z,list(interval)) then

> if nops(z)=2 then

> if z[1]=[] then

> ERROR(`First Element of List must be a real interval.`);

> elif z[2]=[] then

> ERROR(`Second Element of List must be a real interval.`);

> else

> bool:=type(z[1][1],numeric) and type(z[1][2],numeric) and

> type(z[2][1],numeric) and type(z[2][2],numeric);

> fi;

> bool;

> end:

- - in proc init - - - - - - - - - - - - - - -
PROCEDURE `CONVERT/INTERVAL`

A utility program to convert Maple expressions to interval arithmetic.
convert(1+x + x^2,'interval') returns (1 &+ x) &+ (x &^ 2),

whereas inapply(1+x+x^2,x) yields the operator x -> (1 &+ x) etc.

> `convert/interval`:=proc(e)

> local ope, mope, fn;

> # option system:

>

> if type(e,interval) or type(e,interval_comp) then e;

> elif type(e,`+`) then

> if type(op(2,e),`*`) and op(1,op(2,e)) = (-1) then

> `convert/interval`(op(1,e)) &- `convert/interval`(op(1,e)-e);

> else

> `convert/interval`(op(1,e)) &+ `convert/interval`(e-op(1,e));

> fi;

> elif type(e,`*`) then

> if op(1,e)=1.0 then

> `convert/interval`(subs(op(1,e)=1,e));

> else

> `convert/interval`(op(1,e)) &*

> (`convert/interval`(e/op(1,e)));

> fi;

> ### WARNING: note that `I` is no longer of type `^`

> elif type(e,`^`) then

> if type(op(2,e),posint) then

> `convert/interval`(op(1,e)) &intpower op(2,e);

> elif type(op(2,e),integer) then

> inv(`convert/interval`(op(1,e)) &intpower (-op(2,e)));

> elif op(2,e)=(1/2) then

> &sqrt(`convert/interval`(op(1,e)));

> else (`convert/interval`(op(1,e)) &** `convert/interval`(op(2,e)))

> fi;

> elif type(e,function) then

> ope:=[op(e)];

> mope:=op(map(`convert/interval`,ope));

> fn:=subs(Interval_fnlist,op(0,e));

> fn(mope):

> else e

>

> fi:

> end:

   - - in proc init - - - - - - - - - - - - - - -
   PROCEDURE INAPPLY

  Converts the argument to an interval function.

> inapply := proc(); unapply(convert(args[1],'interval'),args[2..nargs]); end:

- - end proc init global type definitions - - - - - - - - - - - - - - -

> with(plots):

> with(geometry):

>
macro(evalf=Evalf);

> RETURN();

> end:

> init();

#### end initialization

###########################################################################

> Evalf := proc(x)

> local e;

> e := evalf(x);

>                 if e = -1.*infinity then -infinity
             elif e = 1.*infinity then infinity
             else e fi

> end:

>

###########################################################################

BASIC INTERVAL DEFINITIONS

-----------------------------------------------------------------------------

Definition of `type/num_or_FAIL`:

## see init procedure!

-----------------------------------------------------------------------------

Definition of `type/interval_comp`:

## see init procedure!

-----------------------------------------------------------------------------

> Interval_ulp := proc(x):

> if x=0 then (0) else

> Float(1,length(op(1,x))+op(2,x)-Digits)

> fi:

> end:

-----------------------------------------------------------------------------

> Interval_Round_Up := proc(x):

> if x=-infinity then x

> elif x=infinity then infinity

> elif x=FAIL then FAIL

> else x + ulp(x)

> fi:

> end:

-----------------------------------------------------------------------------

> Interval_Round_Down := proc(x):

> if x=-infinity then x

> elif x=infinity then infinity

> elif x=FAIL then FAIL

> else x - ulp(x)

> fi:

> end:

-----------------------------------------------------------------------------

Definition of `type/interval`:

## moved to init procedure!

-----------------------------------------------------------------------------

OPERATOR &GREATER/INTERVAL_GREATER

Checks if (Arg1>=Arg2) applies.

> intpakX_greater:=proc(a,b)

> local a_rational, b_rational;

>

> if (a=FAIL) or (b=FAIL) then FAIL;

> elif not(type(a,numeric) or member(a,{infinity,-infinity})) then

> ERROR(`first argument must be a numeric, FAIL, infinity or

> -infinity`);

> elif not(type(b,numeric) or member(b,{infinity,-infinity})) then

> ERROR(`second argument must be a numeric, FAIL, infinity or

> -infinity`);

> elif (a=infinity) then true;

> elif (b=-infinity) then true;

> elif (a = -infinity) then false;

> elif (b=infinity) then false;

> elif a > b then true

> elif type(a,rational) and type(b,rational) then

> if a >= b then true else false fi;

> elif type(a,rational) and type(b,float) then

> b_rational:=op(1,b)*(10^(op(2,b)));

> if a >= b_rational then true else false fi;

> elif type(a,float) and type(b,rational) then

> a_rational:=op(1,a)*(10^(op(2,a)));

> if a_rational >= b then true else false fi;

> elif type(a,float) and type(b,float) then

> a_rational:=op(1,a)*(10^(op(2,a)));

> b_rational:=op(1,b)*(10^(op(2,b)));

> if a_rational >= b_rational then true else false fi;

> fi;

> end:

-----------------------------------------------------------------------------

PROCEDURE - IS_IN

This procedure takes two parameters. It tests to see whether the interval,

or numeric x is contained in the interval a. If a is a float then

the procedure constructs an interval to test.

The operator &greater/Interval_greater is used here.

> is_in:=proc(x,a)

>

> if type(a,interval) then

> if x=FAIL then FAIL

> elif x=infinity then member(infinity,{a[1],a[2]})

> elif x=-infinity then member(-infinity,{a[1],a[2]})

> elif not(type(x,{numeric,interval})) then

> ERROR(`first argument must be a numeric or an interval`)

> elif a=[] and x<>[] then false

> elif a[1]=FAIL or a[2]=FAIL then FAIL

> elif type(x,numeric) then

> intpakX_greater(x,a[1]) and intpakX_greater(a[2],x)

> elif type(x,interval) then

> if x=[] then false #????????????

> elif x[1]=FAIL or x[2]=FAIL then FAIL

> else intpakX_greater(x[1],a[1]) and intpakX_greater(a[2],x[2])

> fi;

> elif a=infinity then evalb(x=infinity)

> elif a=-infinity then evalb(x=-infinity)

> elif a=FAIL then FAIL

> elif type(a,numeric) then

> if x=FAIL then FAIL

> elif x=infinity or x=-infinity then false

> elif type(x,numeric) then

> intpakX_greater(a,x) and intpakX_greater(x,a);

> elif type(x,interval) then

> if x=[] then false

> else intpakX_greater(x[1],a) and intpakX_greater(a,x[2])

> fi;

> else ERROR(`first argument must be a numeric or an interval`)

> fi;

> else ERROR(`second argument must be an interval or a numeric`)

> fi;

> end:

-----------------------------------------------------------------------------

PROCEDURE - CONSTRUCT

This procedure can accept a single argument, to construct a degenerate

interval, two arguments, to construct an interval form the low/high

endpoints. The option `rounded` can be entered as the last argument in

each case to construct a rounded interval.

Special type checking is performed to find (Evalf(-infinity)). Where it

occurs as an argument, Evalf(-infinity) is always rounded to -infinity, to

eliminate Evalf(-infinity) from being an interval endpoint.

Suitable arguments are numeric, FAIL, +-infinity and constants.

> construct:=proc()

> local p,q,Constants;

> if member(false,{args}) then

> ERROR(`false and true ar not accepted as arguments`);

> fi;

>

> if member(FAIL,{args}) then RETURN([FAIL,FAIL]) fi;

> Constants:={Pi,gamma,Catalan};

>

> if nargs = 3 then

> if args[3] = 'rounded' and

> type(args[1],num_or_FAIL) and

> type(args[2],num_or_FAIL) then

> p:=Evalf(args[1]);

> q:=Evalf(args[2]);

> [rd(intpakX_min(p,q)),ru(intpakX_max(p,q))]

> else ERROR(`first,second args must be numeric,third arg must be "rounded"`)

> fi

>

> elif nargs=2 then

> if args[2]=rounded and

> type(args[1],num_or_FAIL) then

> [rd(Evalf(args[1])),ru(Evalf(args[1]))]

> elif type(args[1],num_or_FAIL) and

> type(args[2],num_or_FAIL) then

>

> p:=(args[1]);

> q:=(args[2]);

> if member(p,Constants) and member(q,Constants)

> then construct(p,q,rounded);

> elif member(p,Constants) then

> if Evalf(p) < q then

> if intpakX_greater(Evalf(q),q) then

> [rd(Evalf(p)),Evalf(q)]

> else [rd(Evalf(p)),ru(Evalf(q))]

> fi;

> elif Evalf(p) > q then

> if intpakX_greater(q,Evalf(q)) then

> [Evalf(q),ru(Evalf(p))]

> else [rd(Evalf(q)),ru(Evalf(p))]

> fi;

> else construct(p,q,rounded);

> fi;

> elif member(q,Constants) then

> if Evalf(q) < p then

> if intpakX_greater(Evalf(p),p) then

> [rd(Evalf(q)),Evalf(p)]

> else [rd(Evalf(q)),ru(Evalf(p))]

> fi;

> elif Evalf(q) > p then

> if intpakX_greater(p,Evalf(p)) then

> [Evalf(p),ru(Evalf(q))]

> else [rd(Evalf(p)),ru(Evalf(q))]

> fi;

> else construct(p,q,rounded);

> fi;

> elif intpakX_greater(q,p) then

> if intpakX_greater(Evalf(q),q) and

> intpakX_greater(p,Evalf(p)) then

> [Evalf(p),Evalf(q)];

> elif intpakX_greater(Evalf(q),q) then

> [rd(Evalf(p)),Evalf(q)];

> elif intpakX_greater(p,Evalf(p)) then

> [Evalf(p),ru(Evalf(q))];

> else [rd(Evalf(p)),ru(Evalf(q))];

> fi;

> else

> if intpakX_greater(Evalf(p),p) and

> intpakX_greater(q,Evalf(q)) then

> [Evalf(q),Evalf(p)];

> elif intpakX_greater(Evalf(p),p) then

> [rd(Evalf(q)),Evalf(p)];

> elif intpakX_greater(q,Evalf(q)) then

> [Evalf(q),ru(Evalf(p))];

> else [rd(Evalf(q)),ru(Evalf(p))];

> fi;

>

> fi

> else ERROR(`incorrect arguments entered`)

> fi

>

> elif nargs=1 then

> if type(args[1],num_or_FAIL) then

> p:=args[1];

> if member(p,Constants) then

> [rd(Evalf(p)),ru(Evalf(p))];

> elif is_in(p,Evalf(p)) then

> [Evalf(p),Evalf(p)];

> elif intpakX_greater(p,Evalf(p)) then

> [Evalf(p),ru(Evalf(p))];

> else [rd(Evalf(p)),Evalf(p)];

> fi;

> else ERROR(`incorrect arguments entered`)

> fi

> else ERROR(`maximum of three arguments accepted`)

> fi

> end:

###########################################################################

Logarithm and Exponential Functions:

PROCEDURE &EXP/INTERVAL_EXP

expinfinity is only called from &exp. It deals with FAIL and +/- infinity.

Like most of the other subroutines &exp takes floating point intervals

or numerics (which are converted into intervals).

> expinfinity:=proc(x):

> if x=FAIL then FAIL

> elif x=infinity then infinity

> elif x=-infinity then 0

> else Evalf(exp(x))

> fi:

> end:

> Interval_exp:=proc(x):

> if type(x,'interval') then

> if x=[] then []

> elif x[1]=FAIL or x[2]=FAIL then [FAIL,FAIL]

> else [rd(expinfinity(x[1])),ru(expinfinity(x[2]))]

> fi:

> elif type(x,'num_or_FAIL') then

> Interval_exp(construct(x))

> else

Return unevaluated.

> 'Interval_exp(x)'

ERROR(`floating point interval or scalar arguments are required`)

> fi: end:

---------------------------------------------------------------------------

PROCEDURE &LN

> infinityln:=proc(x)

> if x=FAIL then FAIL

> elif x=infinity then infinity

> elif x=0 then -infinity

> elif (min(x,0)=x ) then FAIL

The above line returns FAIL, as opposed to Maple`s ln function which

returns an ERROR message stating that a singularity has been encountered.

> else ln(x)

> fi:

> end:

> Interval_ln:=proc(x)

> if type(x,'interval') then

> if x=[] then []

> else [rd(infinityln(x[1])),ru(infinityln(x[2]))]

> fi:

> elif type(x,'num_or_FAIL') then Interval_ln(construct(x))

> else

Return unevaluated

> 'Interval_ln(x)'

ERROR (`floating point intervals or scalars required`)

> fi:

> end:

---------------------------------------------------------------------------

###########################################################################

BASIC INTERVAL ARITHMETIC OPERATIONS

All the interval arithmetic subroutines perform type checking.

They accept scalars (of type numeric) or intervals (floating point).

Floating point and integer scalars are made into intervals.

-------------------------------------------------------------------------------

The ***infinity subroutines correct any problems that may occur with

infinite and FAIL results

-------------------------------------------------------------------------------

PROCEDURE &+

> addinfinity:=proc(x,y)

> if (x=infinity and y=(-infinity)) or (x=(-infinity) and y=infinity)

> then FAIL

> elif x=FAIL or y=FAIL then FAIL

> elif x=infinity or y=infinity then infinity

> elif x=-infinity or y=-infinity then -infinity

> else x+y

> fi:

> end:

> Interval_add:=proc(a,b)

> if type(a,'interval') and type(b,'interval') then

> if a=[] or b=[] then []

> else [rd(addinfinity(a[1],b[1])),ru(addinfinity(a[2],b[2]))]

> fi:

> elif type(a,'interval') and type(b,'num_or_FAIL') then Interval_add(a,construct(b))

> elif type(a,'num_or_FAIL') and type(b,'interval') then Interval_add(construct(a),b)

> elif type(a,'num_or_FAIL') and type(b,'num_or_FAIL')

> then Interval_add(construct(a),construct(b))

> else

Return unevaluated

> 'Interval_add'(a,b)

ERROR(`floating point interval and scalar arguments required`)

> fi:

> end:

-------------------------------------------------------------------------------

PROCEDURE &-

> subtractinfinity:=proc(x,y)

> if (x=infinity and y=(infinity)) or (x=(-infinity) and y=(-infinity))

> then FAIL

> elif x=FAIL or y=FAIL then FAIL

> elif x=infinity then infinity

> elif y=infinity then -infinity

> elif x=-infinity then -infinity

> elif y=-infinity then infinity

> else x-y

> fi:

> end:

> Interval_subtract := proc(a,b)

> if type(a,'interval') and type(b,'interval') then

> if a=[] or b=[] then [] else

> [rd(subtractinfinity(a[1],b[2])),ru(subtractinfinity(a[2],b[1]))]

> fi:

> elif type(a,'interval') and type(b,'num_or_FAIL')

> then Interval_subtract(a,construct(b))

> elif type(a,'num_or_FAIL') and type(b,'interval')

> then Interval_subtract(construct(a),b)

> elif type(a,'num_or_FAIL') and type(b,'num_or_FAIL')

> then Interval_subtract(construct(a),construct(b))

> else

Return unevaluated

> 'Interval_subtract(a,b)'

ERROR(`floating point interval and scalar arguments are required`)

> fi:

> end:

-------------------------------------------------------------------------------

PROCEDURE &*

> timesinfinity:=proc(a,b)

> if a=0 or b=0 then 0

> elif a=FAIL or b=FAIL then FAIL

> elif (a=-infinity and min(b,0)=0) or (b=-infinity and min(a,0)=0)then -infinity

> elif (a=-infinity and min(b,0)=b) or (b=-infinity and min(a,0)=a) then infinity

> elif (a=infinity and min(b,0)=0) or (b=infinity and min(a,0)=0) then infinity

> elif (a=infinity and min(b,0)=b) or (b=infinity and min(a,0)=a)then -infinity

the min function is called so that infinity*(-infinity) will give -infinity

if either of the arguments in Interval_times are FAIL the result will

be [FAIL,FAIL].

The above code is long but takes into account every possibility.

> else a*b

> fi:

> end:

> Interval_times := proc(a,b)

> local xy,xY,Xy,XY:

> if type(a,'interval') and type(b,'interval') then

> if a=[] or b=[] then [] else

> xy := timesinfinity(a[1],b[1]):

> xY := timesinfinity(a[1],b[2]):

> Xy := timesinfinity(a[2],b[1]):

> XY := timesinfinity(a[2],b[2]):

> [ rd(min(xy,xY,Xy,XY)),ru(max(xy,xY,Xy,XY))]

> fi:

> elif type(a,'interval') and type(b,'num_or_FAIL') then

> Interval_times(a,construct(b))

> elif type(a,'num_or_FAIL') and type(b,'interval') then

> Interval_times(construct(a),b)

> elif type(a,'num_or_FAIL') and type(b,'num_or_FAIL') then

> Interval_times(construct(a),construct(b))

> else

Return unevaluated

> 'Interval_times(a,b)'

ERROR(`floating point interval and scalar arguments are required`)

> fi:

> end:

---------------------------------------------------------------------

PROCEDURE INV/INTERVAL_RECIPROCAL

This procedure returns 1/infinity =0.

If zero is contained in the denominator(interval) the procedure returns

[-infinity,infinity]

> Interval_reciprocal := proc(x):

> if type(x,'interval') then

> if x=[] then []

> elif x[1]=FAIL or x[2]=FAIL then [FAIL,FAIL]

> elif is_in(0.0,x) then [-infinity,infinity]

> elif abs(x[1])=infinity and abs(x[2])=infinity then [0,0]

> elif abs(x[1])=infinity then [rd(1./x[2]),0]

> elif abs(x[2])=infinity then [0,ru(1./x[1])]

> else [rd(1./x[2]),ru(1./x[1])]

> fi:

> elif type(x,'num_or_FAIL') then

> Interval_reciprocal(construct(x))

> else

Return unevaluated

> 'Interval_reciprocal(x)'

ERROR(`a floating point interval or scalar argument is required`)

> fi:

> end:

----------------------------------------------------------------------

PROCEDURE `&/`/INTERVAL_DIVIDE

This procedure also performs type checking. A check is also done

to see if 0.0 is contained in the denominator. [-infinity,infinity] is

returned.

> Interval_divide := proc(a,b):

> if type(a,'interval') and type(b,'interval') then

> if a=[] or b=[] then []

> elif is_in(0.0,b) then [-infinity,infinity]

> elif (abs(b[1])=infinity or abs(b[2])=infinity) and (abs(a[1])=infinity

> or abs(a[2])=infinity) then [FAIL,FAIL]

> else a &* inv(b)

> fi:

> elif type(a,'interval') and type(b,'num_or_FAIL') then

> Interval_divide(a,construct(b))

> elif type(a,'num_or_FAIL') and type(b,'interval') then

> Interval_divide(construct(a),b)

> elif type(a,'num_or_FAIL') and type(b,'num_or_FAIL') then

> Interval_divide(construct(a),construct(b))

> else

Return unevaluated

> 'Interval_divide(a,b)'

ERROR(`floating point interval and scalar arguments are required`)

> fi:

> end:

---------------------------------------------------------------------

PROCEDURE &SQRT

An error message is returned if a negative argument is entered.

> sqrtinfinity:=proc(x)

> if x=FAIL then FAIL

> elif x=0 then 0

> elif x=infinity then infinity

> elif min(x,0)=x then ERROR (`cannot compute the sqrt of a negative number`)

> else sqrt(x)

> fi:

> end:

> Interval_sqrt := proc(x)

> if type(x,'interval') then

> if x=[] then [] else

> [rd(sqrtinfinity(x[1])),ru(sqrtinfinity(x[2]))]

> fi:

> elif type(x,'num_or_FAIL') then Interval_sqrt(construct(x))

> else

Return unevaluated

> 'Interval_sqrt(x)'

ERROR(`floating point interval or scalar argument is required`)

> fi:

> end:

----------------------------------------------------------------------

PROCEDURE &SQR

> sqrinfinity:=proc(x)

> if x=FAIL then FAIL

> elif abs(x)=infinity then infinity

> else x**2

> fi:

> end:

> Interval_sqr:=proc(x)

> local a,b:

> if type(x,'interval') then

> if x=[] then []

> elif is_in(0,x) then

> [0,ru(max(sqrinfinity(x[1]),sqrinfinity(x[2])))]

> else

> a := min(abs(x[1]),abs(x[2])):

> b := max(abs(x[1]),abs(x[2])):

> [rd(sqrinfinity(a)),ru(sqrinfinity(b))]

> fi:

> elif type(x,'num_or_FAIL') then Interval_sqr(construct(x))

> else

Return unevaluated

> 'Interval_sqr(x)'

ERROR (`a floating point interval or scalar argument is required`)

> fi:

> end:

------------------------------------------------------------------------

> Interval_option_zero:=false:

If option zero is true then 0**0=1

If Interval_option_zero is false then 0**0=FAIL.

This global variable can be changed by the user depending on which they

prefer. The disadvantage with the false case is that even if only one of

the endpoints is zero raised to zero, the resulting interval will

be [FAIL,FAIL].

------------------------------------------------------------------------

PROCEDURE &INTPOWER / INTERVAL_INTEGERPOWER

Integerpower: This takes interval or num_or_FAIL arguments, x, and raises

them to an integer power.

Powerinfinity is a subroutine used to evaluate such cases as

infinity**2 etc.

> powerinfinity:=proc(x,n)

> if n=0 then 1.

> elif n=infinity and min(x,0)=x then FAIL

> elif n=infinity then infinity

> elif n=-infinity then 0

> elif x=0 then 0

This is included to prevent error messages for such cases as 0**(-3)

> elif x=infinity then if n>0 then infinity else 0 fi:

> elif x=-infinity then

> if type(n/2,integer) and n>0 then infinity else (-infinity) fi:

> elif x=FAIL then FAIL

> else x**n

> fi:

> end:

> Interval_Integerpower:=proc(x,n)

> local a,b:

> if type (x,'interval') and type(n,integer) then

> if x=[] then []

> else

> a:=powerinfinity(x[1],n):

> b:=powerinfinity(x[2],n):

The following is a check for monotonicity. If n is even and n>0 then

if zero is in the interval it represents the lowest endpoint in the

returned interval, and the the max endpoint is the max of x[1] and x[2]

to the power n.

Otherwise the function is monotonic, n>0 and the endpoints are evaluated

directly for the maximum, and minimum values the function takes on the

interval,x.

If n<0 and 0 is contained in the interval, x, then [-infinity,infinity]

is returned. Otherwise the value of the two endpoints raised to the negative

power n are returned.

> if n<0 and type(n/2,integer) and is_in(0,x) then [rd(min(a,b)),infinity]

> elif n<0 and (not type(n/2,integer)) and is_in(0,x) then [-infinity,infinity]

> elif (n<0) then [rd(min(a,b)),ru(max(a,b))]

> elif n>0 and type(n/2,integer) and is_in(0,x) then

> [0,ru(max(a,b))]

> else [rd(min(a,b)),ru(max(a,b))]

The else case covers n=0, n>0 and odd (ie monotonic).

> fi:

> elif type(x,'num_or_FAIL') then Interval_Integerpower(construct(x),n)

> else

Return unevaluated

> 'Interval_Integerpower(x,n)'

ERROR

(`arg[1] must be a float interval or numeric, arg[2] must be an integer`)

> fi:

> end:

-----------------------------------------------------------------------------

PROCEDURE &** / INTERVAL_POWER

This procedure calculates an interval raised to the power of another interval.

Digits is extended to reduce rounding error.

Note the conditions on the ops of x, that they be numeric, so that the

ilog10 function can be applied.

The procedure has been enhanced to accept numbers of type fraction as exponents.

> Interval_power:=proc(x,n)

> local logx,prod,result,oldDigits;

> if type(x,interval) and (type(n,interval)
or type(n,interval_comp) or type(n,fraction)) then

> oldDigits:=Digits;

> if x = [] then []

> elif n = [] then []

> elif not(type(x[1],numeric) and type(x[2],numeric)) then Digits:=Digits

> elif not type(x[1],numeric) then

> Digits:=Digits + 2 + intpakX_ilog10(x[2]);

> elif not type(x[2],numeric) then

> Digits:=Digits + 2 + intpakX_ilog10(x[1])

> else

> Digits:=intpakX_max(intpakX_ilog10(x[1]),intpakX_ilog10(x[2])) + Digits + 2

> fi;

> if type(n,interval) and (abs(n[1])=infinity or abs(n[2])=infinity)
and (x[1]=0 or x[2]=0) then [FAIL,FAIL]

> elif type(n,interval_comp) and abs(n)=infinity and (x[1]=0 or x[2]=0) then
[FAIL,FAIL]

> elif type(n,interval) and (n[1]=0 or n[2]=0) and (x[1]=0 or x[2]=0)
and not Interval_option_zero then [FAIL,FAIL]

> elif type(n,interval_comp) and n=0 and (x[1]=0 or x[2]=0)
and not Interval_option_zero then [FAIL,FAIL]

> elif type(n,fraction) then

> if n=(1/2) then &sqrt(x)

> else

> logx:=&ln(x);

> prod:=construct(n,rounded) &* logx;

> result:=&exp(prod);

> Digits:=oldDigits;

> construct(result[1],result[2],rounded)

> fi

> else

> logx:=&ln(x);

> prod:=n &* logx;

> result:= &exp(prod);

> Digits:=oldDigits;

> construct(result[1],result[2],rounded)

> fi

> elif type(x,num_or_FAIL) then construct(x) &** n

> elif type(n,integer) then x &intpower n

> else 'x &** n'

> fi

> end:

>

------------------------------------------------------------------------------

PROCEDURE MIDPOINT / INTERVAL_MIDPOINT

> Interval_midpoint:=proc(x):

> if type(x,'interval') then

> if x=[] then FAIL

> elif x[1]=FAIL or x[2]=FAIL then FAIL

> else (addinfinity(x[1],x[2]))&/2.:

> fi:

> elif type (x,'num_or_FAIL') then Interval_midpoint(construct(x))

> else

Return unevaluated

> 'Interval_midpoint(x)'

ERROR(`floating interval or scalar argument required`)

> fi:

> end:

------------------------------------------------------------------------

PROCEDURE WIDTH / INTERVAL_WIDTH

> Interval_width:=proc(x):

> if type(x,'interval') then

> if x=[] then []

> elif x[1]=FAIL or x[2]=FAIL then [FAIL,FAIL]

> else subtractinfinity(x[2],x[1])

> fi:

> elif type(x,'num_or_FAIL') then

> if not (x=infinity or x=-infinity) then 0

> else Interval_width(construct(x)) fi:

> else

Return unevaluated

> 'Interval_width(x)'

ERROR(`floating interval or scalar argument required`)

> fi:

> end:

###########################################################################

Union and Intersect

------------------------------------------------------------------------------

PROCEDURE &INTERSECT / INTERVAL_INTERSECT

> Interval_intersect:=proc(a,b)

> if type(a,'interval') and type(b,'interval') then

> if a=[] or b=[] then []

> elif a[1]=FAIL or a[2]=FAIL or b[1]=FAIL or b[2]=FAIL then [FAIL,FAIL]

> elif (max(a[1],b[2])=a[1] and not a[1]=b[2])

> or (max(b[1],a[2])=b[1] and not b[1]=a[2]) then []

> else [(max(a[1],b[1])),(min(a[2],b[2]))]

> fi:

> elif type(a,'num_or_FAIL') and type(b,'interval') then

> Interval_intersect(construct(a),b)

> elif type(a,'interval') and type(b,'num_or_FAIL') then

> Interval_intersect(a,construct(b))

> elif type(a,'num_or_FAIL') and type(b,'num_or_FAIL') then

> Interval_intersect(construct(a),construct(b))

> else

Return unevaluated

> 'Interval_intersect(a,b)'

ERROR(`floating point interval or scalar arguments required`)

> fi:

> end:

-------------------------------------------------------------------------

PROCEDURE &UNION / INTERVAL_UNION

This accepts floating point interval or scalar arguments. If the intersection

of the two arguments is the empty interval then the union of the two arguments

are the arguments themselves.

> Interval_union:=proc(a,b)

> if type(a,'interval') and type(b,'interval') then

> if a=[] then b

> elif b=[] then a

> elif a[1]=FAIL or a[2]=FAIL or b[1]=FAIL or b[2]=FAIL then [FAIL,FAIL]

> elif &intersect(a,b)=[] then RETURN(a,b)

> else [(min(a[1],b[1])),(max(a[2],b[2]))]

> fi:

> elif type(a,'num_or_FAIL') and type(b,'interval') then

> Interval_union(construct(a),b)

> elif type(a,'interval') and type(b,'num_or_FAIL') then

> Interval_union(a, construct(b))

> elif type(a, 'num_or_FAIL') and type(b,'num_or_FAIL') then

> Interval_union(construct(a), construct(b))

> else

Return unevaluated

> 'Interval_union(a,b)'

ERROR(`floating point interval or scalar arguments required`)

> fi:

> end:

###############################################################################

TRIGONOMETRIC FUNCTIONS
(sin, cos, tan and corresponding arc and hyperbolic procedures)

----------------------------------------------------------------------------

Auxiliary functions

The following two procedures are called from &cos and &sin to ensure

that the maximum and minimum the result can take are 1. and -1.

There is no rounding out of 1. and -1.

> Interval_trig_ru:=proc(x) if x= 1. then 1. else ru(x) fi; end:

> Interval_trig_rd:=proc(x) if x=-1. then -1. else rd(x) fi; end:

----------------------------------------------------------------------------

The following was written by Dr Dave Hare.

> intpakX_ilog10:=proc(x) if x=0 then 0 else length(op(1,x))+op(2,x)-1; fi; end:

----------------------------------------------------------------------------

The following code is similar to Dr Dave Hare's code for scaling in

Maple's sin function, except the arguments are scaled down

as follows:

k:=xx/2*Pi and x:=xx-2*Pi*k , where x is returned.

In order that the shape of the interval be preserved, the arguments

are divided by 2*Pi.

The two endpoints of the interval argument, xx, are scaled to the

same degree, to the minimum k value.

After scaling the interval 'result' is rounded out.

> Interval_scale:=proc(xx)

> local k,x,d,k1,k2, mag_x1,mag_x2,mag_x,result:

> mag_x1:=ilog10(xx[1]):

> mag_x2:=ilog10(xx[2]):

> mag_x:=max(mag_x1,mag_x2):

> if mag_x1 >=0 then

> d:=max(Digits, 2+mag_x):

> k1:=round(Evalf(xx[1]/(2*Pi),d)):

> fi:

> if mag_x2>=0 then

> d:=max(Digits,2+mag_x):

> k2:=round(Evalf(xx[2]/(2*Pi),d)):

> fi:

> if mag_x1>=0 and mag_x2>=0 then

> k:=min(k1,k2):

> x[1]:=rd(Evalf((xx[1])-2*Pi*k,Digits+d+mag_x)):

> x[2]:=ru(Evalf((xx[2])-2*Pi*k,Digits+d+mag_x)):

> result:= [x[1],x[2]]:

> else result:= [xx[1],xx[2]] fi:

> result

> end:

------------------------------------------------------------------

> Interval_range_values:=proc(x,y1,y2)

> local k1,k2,p1,p2,int:

The following code is used to test for maxima and minima in the &sin

and &cos subroutines.

(In order for this testing to work for &cos, Pi/2 must first be added to the

arguments. This is done in &cos.)

First the number of the 2Pi interval, k, in which the endpoint occurs

is assigned. The 2Pi intervals either side of x=zero have value k=0.

In the positive x axis k>=0 , in the negative x axis k<=0. It is then

possible to compare the k1 and k2 values. Testing is done for both the

positive and negative cases if abs(k2-k1)=1. If x[2]-x[1]>=2*Pi then

obviously there is a maximum and a minimum in the interval and [-1.,1.]

is returned.

> k1:=trunc(Evalf(x[1]/(2.*Pi))):

> k2:=trunc(Evalf(x[2]/(2.*Pi))):

The interval(x) is scaled down to `int`. Int is then tested for containment

of 1/4 and 3/4, those fractions of the 2*Pi interval at which the maxima and

minima occur.

> p1:=frac(Evalf(x[1]/(2.*Pi))):

> p2:=frac(Evalf(x[2]/(2.*Pi))):

> int:=construct(p1,p2):

> if abs(x[1]-x[2])>=Evalf(2*Pi) then [-1.,1.]

> elif abs(k2-k1)=1 then

> if Evalf(x[1])>=0 and Evalf(x[2])>=0 then

> if p1<=(1/4) then [-1.,1.]

> elif p1>=(1/4) and p1<=(3/4) and p2>=1/4 then [-1.,1. ]

> elif p1>=(1/4) and p1<=3/4 and p2<=(1/4) then [-1.,Interval_trig_ru(max(y1,y2))]

> elif p1>=3/4 and p2>=1/4 and p2<=3/4 then [Interval_trig_rd(min(y1,y2)),1.]

> elif p1>=3/4 and p2>=3/4 then [-1.,1.]

> elif p1>=3/4 and p2<=1/4 then [rd(min(y1,y2)),ru(max(y1,y2))]

> fi:

The case for all negative endpoints

> elif Evalf(x[1])<=0 and Evalf(x[2])<=0 then

> if p1<=-3/4 then [-1.,1.]

> elif p1>-3/4 and p1<=-1/4 and p2>=-3/4 then [-1.,1.]

> elif p1>-3/4 and p1<=-1/4 and p2<=-3/4 then [-1.,ru(max(y1,y2))]

> elif p1>-1/4 and p2>=-3/4 and p2<=-1/4 then [rd(min(y1,y2)),1.]

> elif p1>-1/4 and p2>=-1/4 then [-1.,1.]

> elif p1>-1/4 and p2<=-3/4 then [rd(min(y1,y2)),ru(max(y1,y2))]

> fi:

This is the case for k values the same, on the positive x axis.

> elif x[1]>=0 and x[2]>=0 then

> if is_in(1/4,int) and is_in(3/4,int) then [-1.,1.]

> elif is_in(1/4,int) then [Interval_trig_rd(min(y1,y2)),1.]

> elif is_in(3/4,int) then [-1.,Interval_trig_ru(max(y1,y2))]

> else [Interval_trig_rd(min(y1,y2)),Interval_trig_ru(max(y1,y2))]

> fi:

This is the case for k values the same, on the negative x axis.

> elif x[1]<=0 and x[2]<=0 then

> if is_in(-3/4,int) and is_in(-1/4,int) then [-1.,1.]

> elif is_in(-3/4,int) then [Interval_trig_rd(min(y1,y2)),1.]

> elif is_in(-1/4,int) then [-1.,Interval_trig_ru(max(y1,y2))]

> else [Interval_trig_rd(min(y1,y2)),Interval_trig_ru(max(y1,y2))]

> fi:

This is the case for k values the same, 0, either side of x=0.

> elif x[1]<=0 and x[2]>=0 then

> if is_in(-3/4,int) or is_in(3/4,int) then [-1.,1.]

> elif is_in(-1/4,int) and is_in(1/4,int) then [-1.,1.]

> elif is_in(-1/4,int) then [-1.,Interval_trig_ru(max(y1,y2))]

> elif is_in(1/4,int) then [Interval_trig_rd(min(y1,y2)),1.]

> else [Interval_trig_rd(min(y1,y2)),Interval_trig_ru(max(y1,y2))]

> fi:

> end:

--------------------------------------------------------------------------

PROCEDURE &SIN / INTERVAL_SIN

This function accepts floating-point interval arguments, or floating point

or integer scalar arguments.

> Interval_sin:=proc(s)

> local x,y1,y2:

> if type (s,'interval') then

> if s=[] then []

> elif s[1]=-infinity or s[2]=-infinity or s[1]=infinity

> or s[2]=infinity then [-1.,1.]

> elif s[1]=FAIL or s[2]=FAIL then [FAIL,FAIL]

The following test ensures that if the argument entered at either of the

endpoints has error in its last digit, then if the ulp of the last digit is

greater than 2*Pi, [-1.,1.] is returned.

> elif (ulp(s[1])/10.>=Evalf(2*Pi) or ulp(s[2])/10.>=Evalf(2*Pi)) and

> not s[1]=s[2] then [-1.,1.]

> else

> x:=Interval_scale(s):

> y1:=Evalf(sin(x[1])):

> y2:=Evalf(sin(x[2])):

> Interval_range_values(x,y1,y2):

> fi:

> elif type(s,'num_or_FAIL') then Interval_sin(construct(s))

> else

Return unevaluated

> 'Interval_sin(s)'

ERROR(` floating point intervals or scalars required`)

> fi:

> end:

---------------------------------------------------------------------------

PROCEDURE &COS / INTERVAL_COS

This function uses identical testing code from &sin to test for maximas

and minimas by adding Pi/2 to the argument of the cos function.

This saves repetition of code that has already been written and tested.

The actual values of cos at the endpoints are calculated before the Pi/2 is

added.

This function accepts floating point interval arguments, or floating point

or integer scalar arguments.

> Interval_cos:=proc(s)

> local y1,y2,r,t,x:

> if type (s,'interval') then

> if s=[] then []

> elif s[1]=-infinity or s[2]=-infinity or s[1]=infinity

> or s[2]=infinity then [-1.,1.]

> elif s[1]=FAIL or s[2]=FAIL then [FAIL,FAIL]

The same test as in &sin occurs here.

> elif (ulp(s[1])>=Evalf(2*Pi) or ulp(s[2])>=Evalf(2*Pi)) and

> not s[1]=s[2]

> then [-1.,1.]

> else

> r:=Interval_scale(s):

Digits are increased for accuracy.

> t:=[Evalf(s[1]+Evalf(.5*Pi,Digits+3),Digits+3),Evalf(s[2]

> +Evalf(.5*Pi,Digits+3),Digits+3)]:

> x:=Interval_scale(t):

> y1:=Evalf(cos(r[1])):

> y2:=Evalf(cos(r[2])):

> Interval_range_values(x,y1,y2):

> fi:

> elif type(s,'num_or_FAIL') then Interval_cos(construct(s))

> else

Return unevaluated

> 'Interval_cos(s)'

ERROR(` floating point intervals or scalars required`)

> fi:

> end:

--------------------------------------------------------------------------

PROCEDURE &TAN / INTERVAL_TAN

This function accepts floating-point interval arguments, or floating point

or integer scalar arguments.

For ease of testing, like &cos, Pi/2 is added, with increased Digits.

In the same way as for &cos and &sin, the interval is scaled down to a 2*Pi

interval and the k1 and k2 values are calculated.

The testing for this function is simpler than the &sin and &cos due

to the monotonicity of tan over the regions on which it is defined.

> Interval_tan:=proc(s)

> local int,k1,k2,p1,p2,y1,y2,r,t,x:

> if type (s,'interval') then

> if s=[] then []

> elif s[1]=-infinity or s[2]=-infinity or s[1]=infinity

> or s[2]=infinity then [-infinity,infinity]

> elif s[1]=FAIL or s[2]=FAIL then [FAIL,FAIL]

The following test ensures that if the argument entered at either of the

two endpoints (providing it's not a degenerate interval) has ulp greater than

Pi then [-infinity,infinity] is returned.

> elif (ulp(s[1])>=Evalf(Pi) or ulp(s[2])>=Evalf(Pi)) and

> not s[1]=s[2] then [-infinity,infinity]

> else

> r:=Interval_scale(s):

> t:=[Evalf(s[1]+Evalf(.5*Pi),Digits+3),Evalf(s[2]

> +Evalf(.5*Pi),Digits+3)]:

> x:=Interval_scale(t):

> y1:=Evalf(tan(r[1])):

> y2:=Evalf(tan(r[2])):

First the number of the 2Pi interval, k, in which

the endpoint occurs is assigned. The 2Pi intervals either side of x=zero

have value k=0. In the positive x axis k>=0 , in the negative x axis k<=0.

It is then possible to compare the k1 and k2 values.

> k1:=trunc(Evalf(x[1]/(2.*Pi))):

> k2:=trunc(Evalf(x[2]/(2.*Pi))):

The fraction of the 2*Pi interval is then assigned.

The interval(x) is scaled down to `int`. Int is then tested for containment

of 1/2, that fraction of the 2*Pi interval at which the singularity is

encountered.

> p1:=frac(Evalf(x[1]/(2.*Pi))):

> p2:=frac(Evalf(x[2]/(2.*Pi))):

> int:=construct(p1,p2):

If x[2]-x[1]>=Pi then obviously there is a singularity encountered in the

interval and [-infinity,infinity] is returned.

> if abs(x[1]-x[2])>=Evalf(Pi) then [-infinity,infinity]

> elif abs(k2-k1)=1 then [-infinity,infinity]

After scaling, x=0 is a singularity.

> elif k1=0 and k2=0 and x[1]<=0 and x[2]>=0 then [-infinity,infinity]

> elif k1=k2 and x[1]>=0 and x[2]>=0 then

If one of the endpoints is a singularity , for example [-Pi/2,0], then

[-infinity,infinity] is returned.

> if is_in(.5,int) or is_in (-.5,int) or is_in(0,int)

> then [-infinity,infinity]

> else [Interval_trig_rd(min(y1,y2)),Interval_trig_ru(max(y1,y2))]

> fi:

> elif k1=k2 and x[1]<=0 and x[2]<=0 then

> if is_in(.5,int) or is_in(-.5,int) or is_in (0,int)

> then [-infinity,infinity]

> else [Interval_trig_rd(min(y1,y2)), Interval_trig_ru(max(y1,y2))]

> fi:

> elif type(s,'num_or_FAIL') then Interval_tan(construct(s))

> else

Return unevaluated

> 'Interval_tan(s)'

ERROR(` floating point intervals or scalars required`)

> fi:

> end:

-------------------------------------------------------------------------------

PROCEDURE &ARCSIN / INTERVAL_ARCSIN

The following simple function returns the rounded interval of the

inverse sin function. The Maple V arcsin function carries one guard digit.

Here simply rounding the interval out ensures that the solution(s)

is/are always contained in the interval.

The results returned are in the interval [-Pi/2,Pi/2].

> Interval_arcsin:=proc(x):

> if type (x,'interval') then

> if x=[] then []

> elif x[1]=FAIL or [2]=FAIL then [FAIL,FAIL]

The testing for the correct argument range is complicated by the possiblity

of infinite arguments.

> elif (max(abs(x[1]),abs(1.))=abs(x[1]) and not (x[1]=1. or x[1]=-1.))

> or (max(abs(x[2]), abs(1.))=abs(x[2]) and not (x[2]=1. or x[2]=-1.)) then

> ERROR ( `the arguments must be in the range [-1.,1.]`)

> else [rd(Evalf(arcsin(x[1]))),ru(Evalf(arcsin(x[2])))]:

> fi:

> elif type(x,'num_or_FAIL') then

> Interval_arcsin(construct(x))

> else

Return unevaluated

> 'Interval_arcsin(x)'

ERROR(` floating point intervals or scalars required`)

> fi:

> end:

----------------------------------------------------------------------------

PROCEDURE &ARCCOS / INTERVAL_ARCCOS

The following simple function returns the rounded interval of the

arccos function. The Maple V arccos function carries one guard digit.

Here simply rounding the interval out ensures that the solution(s)

is/are always contained in the interval.

The results returned are in the interval [0, Pi].

> Interval_arccos:=proc(x):

> if type (x,'interval') then

> if x=[] then []

> elif x[1]=FAIL or [2]=FAIL then [FAIL,FAIL]

The testing for the correct argument range is complicated by the possiblity

of infinite arguments.

> elif (max(abs(x[1]),abs(1.))=abs(x[1]) and not (x[1]=1. or x[1]=-1.))

> or (max(abs(x[2]), abs(1.))=abs(x[2]) and not (x[2]=1. or x[2]=-1.)) then

> ERROR ( `the arguments must be in the range [-1.,1.]`)

> else [rd(Evalf(arccos(x[2]))),ru(Evalf(arccos(x[1])))]:

> fi:

> elif type(x,'num_or_FAIL') then

> Interval_arccos(construct(x))

> else

Return unevaluated

> 'Interval_arccos(x)'

ERROR(` floating point intervals or scalars required`)

> fi:

> end:

------------------------------------------------------------------------------

PROCEDURE &ARCTAN / INTERVAL_ARCTAN

The following procedure calculates the range of the arctan function,

given the x coordinate and the y coordinate of a point.

The default value for x is 1, and the answer returned is in the

[-Pi/2,Pi/2] interval. The arctan of the high and low endpoints of the

interval arguments are calculated. The rounded, widest possible interval

is returned.

> Interval_arctan:=proc(y,x)

> local a,b,c,d:

> if nargs=2 then

> if type(y,'interval') and (type(x,'interval') or x=[1,1]) then

The integer interval is needed for cases where only one argument is given

and [1,1] is the default value for the second argument. It is important for

cases involving infinity that the 1 be an integer as Maple's arctan function

returns cannot calculate such examples as arctan(infinity,1.0).

This may be amended in the next version of Maple, as a result of suggestions

to the writers of maple.

> if x[1]=FAIL or x[2]=FAIL or y[1]=FAIL or y[2]=FAIL then [FAIL,FAIL]

> else

> a:=Evalf(arctan(y[1],x[1])):

> b:=Evalf(arctan(y[1],x[2])):

> c:=Evalf(arctan(y[2],x[1])):

> d:=Evalf(arctan(y[2],x[2])):

> [rd(min(a,b,c,d)),ru(max(a,b,c,d))]

> fi:

> elif type(y,'interval') and type(x,'num_or_FAIL') then

> Interval_arctan(y,construct(x))

> elif (type(x,'interval') or x=[1,1]) and type(y,'num_or_FAIL') then

> Interval_arctan(construct(y),x)

> elif type(x,'num_or_FAIL') and type(y,'num_or_FAIL') then

> Interval_arctan(construct(y),construct(x))

> fi:

> elif nargs=1 then Interval_arctan(y,[1,1])

> else

Return unevaluated

> 'Interval_arctan(y,x)'

ERROR(`up to two floating point interval or scalar arguments accepted`)

> fi:

> end:

-------------------------------------------------------------------------------

Auxiliary procedures for sinh, cosh, tanh

The following rounding procedure ensure that the results of the hyperbolic

functions are always in the correct range. e.g. -1.<= tanh(x) <=1.

1.<=cosh(x).

> Interval_hyp_rd:=proc(x); if x=-1. then x else rd(x) fi; end:

> Interval_hyp_ru:=proc(x); if x= 1. then x else ru(x) fi; end:

------------------------------------------------------------------------------

PROCEDURE &COSH / INTERVAL_COSH

The following simple function returns the rounded interval result

of the cosh function over an interval range. Maple's cosh function is

called. A check is done for the inclusion of 0, a minimum.

> Interval_cosh:=proc(x):

> if type(x,'interval') then

> if x=[] then []

> elif x[1]=FAIL or x[2]=FAIL then [FAIL,FAIL]

> elif is_in(0,x) then [1.,ru(Evalf(max(cosh(x[1]),cosh(x[2]))))]

> else [Interval_hyp_rd(Evalf(cosh(x[1]))),ru(Evalf(cosh(x[2])))]

The absolute minimum of the cosh function is 1., so the lower endpoint is not

rounded below 1.

> fi:

> elif type (x,'num_or_FAIL') then Interval_cosh(construct(x))

> else

Return unevaluated

> 'Interval_cosh(x)'

ERROR (`floating point interval or scalar argument required`)

> fi:

> end:

------------------------------------------------------------------------------

PROCEDURE &SINH / INTERVAL_SINH

The following simple function returns the rounded interval result

of the sinh function over an interval range. Maple's sinh function

is called.

> Interval_sinh:=proc(x):

> if type(x,'interval') then

> if x=[] then []

> elif x[1]=FAIL or x[2]=FAIL then [FAIL,FAIL]

> else [rd(sinh(x[1])),ru(sinh(x[2]))]

> fi:

> elif type(x,'num_or_FAIL') then Interval_sinh(construct(x))

> else

Return unevaluated

> 'Interval_sinh(x)'

ERROR (`floating point interval or scalar argument required`)

> fi:

> end:

------------------------------------------------------------------------------

PROCEDURE &TANH / INTERVAL_TANH

The following simple function returns the rounded interval result

of the tanh function over an interval range. Maple's tanh function

is called. Results are in the range [-1.,1.].

> Interval_tanh:=proc(x):

> if type(x,'interval') then

> if x=[] then []

> elif x[1]=FAIL or x[2]=FAIL then [FAIL,FAIL]

> else [Interval_hyp_rd(Evalf(tanh(x[1]))),Interval_hyp_ru(Evalf(tanh(x[2])))]

> fi:

> elif type(x,'num_or_FAIL') then Interval_tanh(construct(x))

> else

Return unevaluated

> 'Interval_tanh(x)'

ERROR (`floating point interval or scalar argument required`)

> fi:

> end:

###############################################################################

PROCEDURES `CONVERT/INTERVAL` AND INAPPLY

-----------------------------------------------------------------------------

Definitions of `convert/interval` and inapply:

--> moved to init procedure!

-----------------------------------------------------------------------------

###############################################################################

Minimum and Maximum

-----------------------------------------------------------------------------

PROCEDURE MIN / INTPAKXI_MIN

This is called from &*, and numerous other procedures. It finds the minimum

of n arguments, of type(numeric) or infinity.

Note if one the n arguments given is FAIL then FAIL is returned

> intpakX_min:=proc() local a,i,result:

> a:={args} minus {infinity}:

> if nops(a)=0 then RETURN(infinity)

> elif member(FAIL,a) then RETURN (FAIL)

> elif member(-infinity,a) then RETURN (-infinity)

> else

> result :=a[1];

> for i from 2 to nops(a) do

> if a[i]<result then result := a[i] fi;

> od:

> result

> fi:

> end:

-----------------------------------------------------------------------------

PROCEDURE MAX / INTPAKXI_MAX

Note if FAIL is one of the n arguments given FAIL is returned.

> intpakX_max:=proc() local a,i,result:

> a:={args} minus {-infinity}:

> if nops(a)=0 then RETURN(-infinity)

> elif member (FAIL,a) then RETURN(FAIL)

> elif member (infinity,a) then RETURN(infinity)

> else

> result:=a[1];

> for i from 2 to nops(a) do

> if a[i]>result then result:=a[i] fi:

> od:

> result

> fi:

> end:

###########################################################################

EXTENSIONS (from intpakX by I. Geulig and W. Kraemer)

###########################################################################

03/2002: If there is any code or output in German left in this part, please contact me at

markus.grimmer@math.uni-wuppertal.de

-----------------------------------------------------------------------------

PROCEDURE EXT_INT_DIV

Extended Interval Division.

> ext_int_div:=proc(x,y)

>

> if (type(x,interval) and type(y,interval)) then

> if not( type(x[1],numeric) and type(x[2],numeric) ) then

> ERROR(`first arg must be a finite interval or a numeric`);

> fi;

>

> if not( type(y[1],numeric) and type(y[2],numeric) ) then

> ERROR(`second arg must be a finite interval or a numeric`);

> fi;

>

> if is_in(0,y) then

> if is_in(0,x) then

> [-infinity,infinity]

> elif (y[1]=0) and (y[2]=0) then []

> elif (x[2]<0) and (y[1]<y[2]) and (y[2]=0) then

> [rd(x[2]/y[1]),infinity]

> elif (x[2]<0) and (y[1]<0) and (0<y[2]) then

> [-infinity,ru(x[2]/y[2])] &union [rd(x[2]/y[1]),infinity]

> elif (x[2]<0) and (0=y[1]) and (y[1]<y[2]) then

> [-infinity,ru(x[2]/y[2])]

> elif (0<x[1]) and (y[1]<y[2]) and (y[2]=0) then

> [-infinity,ru(x[1]/y[1])]

> elif (0<x[1]) and (y[1]<0) and (0<y[2]) then

> [-infinity,ru(x[1]/y[1])] &union [rd(x[1]/y[2]),infinity]

> elif (0<x[1]) and (0=y[1]) and (y[1]<y[2]) then

> [rd(x[1]/y[2]),infinity];

> fi;

>

> else (x &/ y)

> fi;

>

> elif type(x,num_or_FAIL) and type(y,interval) then

> ext_int_div(construct(x),y);

> elif type(x,interval) and type(y,num_or_FAIL) then

> ext_int_div(x,construct(y));

> elif type(x,num_or_FAIL) and type(y,num_or_FAIL) then

> ext_int_div(construct(x),construct(y));

> else (x &/ y)

> fi;

> end:

-----------------------------------------------------------------------------

PROCEDURE REL_DIAM

Relative Diameter of an interval.

> rel_diam:=proc(x::interval)

> if x[1]=FAIL or x[2]=FAIL then

> ERROR(`Interval must be real.`);

> else

> if (-infinity < x[1]) and (x[2] < infinity) then

> if is_in(0,x) then

> width(x);

> else

> width(x) / (min(abs(x[1]),abs(x[2])));

> fi:

> else

> ERROR(`Interval must be real.`);

> fi:

> end:

-----------------------------------------------------------------------------

PROCEDURE MID

Midpoint of a real interval. The procedure returns a number which is guaranteed to be in the input interval.

> mid:=proc(z)

> local m;

> if type(z,interval) then

> if z=[] or z[1]=FAIL or z[2]=FAIL or

> is_in(infinity,z) or is_in(-infinity,z) then

> ERROR(` Parameter must be real Interval! `);

> else

> if width(z)=0 then

> m:=z[1];

> else

> m:=z[1] + 0.5 * width(z);

> if not( (z[1] < m) and (m < z[2])) then

> m:=z[2];

> fi;

> elif type(z,numeric) then

> m:=z;

> else

> RETURN('mid(z)');

> fi;

> RETURN(m);

> end:

###########################################################################

RANGE ENCLOSURE WITH GRAPHICAL REPRESENTATION

###########################################################################

FUNCTIONS OF ONE REAL VARIABLE

-----------------------------------------------------------------------------

PROCEDURE COMPUTE_NAIVE_INTERVAL_RANGE

Range enclosure using interval evaluation.

> compute_naive_interval_range:=proc(f,T)

> local F;

> F:=inapply(f(t),t);

> RETURN(F(T));

> end:

-----------------------------------------------------------------------------

PROCEDURE COMPUTE_MEAN_VALUE_RANGE

Range enclosure using mean value form.

compute_mean_value_range:=proc(f,t_0,T)

> local dF, F_Bound, MeanValueForm;

> dF:=inapply(diff(f(t),t),t);

> F_Bound:=dF(T);

> MeanValueForm:=inapply(f(t),t)(t_0) &+ (F_Bound &* (T &- t_0));

> RETURN(MeanValueForm);

> end:

-----------------------------------------------------------------------------

PROCEDURE COMPUTE_MONOTONIC_RANGE

Range Enclodure of a monotone twice continuously differentiable function.
If monotony cannot be proven, FAIL will be returned.

> compute_monotonic_range:=proc(f,T::interval)

> local df,dfBound,fLeft,fRight,c;

> df:=D(f);

> dfBound:=compute_naive_interval_range(df,T);

> fLeft:=compute_naive_interval_range(f,construct(T[1]));

> fRight:=compute_naive_interval_range(f,construct(T[2]));

>

> if dfBound[2] <= 0. then

> # Function is monotonously decreasing

> RETURN([fRight[1],fLeft[2]]);

> elif 0. <= dfBound[1] then

> # Function is monotonously increasing

> RETURN([fLeft[1],fRight[2]]);

> fi;

>

> c:=T[1] + 0.5 * (T[2]-T[1]);

> dfBound:=dfBound &intersect compute_mean_value_range(df,c,T);

>

> if dfBound[2] <= 0. then

> # Function is monotonously decreasing

> RETURN([fRight[1],fLeft[2]]);

> elif 0. <= dfBound[1] then

> # Function is monotonously increasing

> RETURN([fLeft[1],fRight[2]]);

> fi;

>

> # Monotony couldn't be proven

>

> RETURN(FAIL)

> end:

>

-----------------------------------------------------------------------------

PROCEDURE COMBINED_RANGE

Range enclosure combining the above three methods

> compute_combined_range:=proc(f,T)

> local fBound,c;

> fBound:=compute_monotonic_range(f,T);

>

> if not(fBound = FAIL) then

> RETURN(fBound)

> fi;

>

> c:=T[1] + 0.5*(T[2]-T[1]);

> fBound:=compute_naive_interval_range(f,T)

> &intersect compute_mean_value_range(f,c,T);

>

> RETURN(fBound);

> end:

>

-----------------------------------------------------------------------------

PROCEDURE COMPUTE_TAYLOR_FORM_RANGE

Range enclosure using second order Taylor form

> compute_taylor_form_range:=proc(f,t_0,T)

> local dF2, tf_bound, taylor_form;

> dF2:=inapply(diff(f(t),t$2),t);

> tf_bound:=dF2(T);

> taylor_form:=inapply(f(t),t)(t_0) &+

> (inapply(diff(f(t),t),t)(t_0) &* (T &- t_0)) &+

> ((tf_bound &/ 2) &* ((T &- t_0) &intpower 2));

> RETURN(taylor_form);

> end:

-----------------------------------------------------------------------------

PROCEDURES &CONVEX_HULL

Computes convex hull of two intervals.

> `&Convex_Hull`:=proc(x::interval,y::interval)

> RETURN([min(x[1],y[1]),max(x[2],y[2])]);

> end:

-----------------------------------------------------------------------------

PROCEDURE SUBDIVIDE_EQUIDISTANT, SUBDIVIDE_ADAPTIVE

Methods for interval subdivision during the computation of range enclosures.

> subdivide_equidistant:=proc(f::{function,procedure},

> list_of_intervals::evaln,

> list_of_ranges::evaln)

> local N, n, i, iter_counter, T, T1, T2, F1, F2, temp, t_mid,

> temp_list_of_intervals, temp_list_of_ranges;

>

> N:=3;

>

> if not(type(eval(list_of_intervals),listlist)) then

> lprint(`Second Parameter must be list of intervals!`);

> T:=readstat(`Enter Start interval: `);

> list_of_intervals:=[construct(T[1],T[2])];

> fi;

>

> n:=nops(eval(list_of_intervals));

>

> temp_list_of_intervals:=[];

> temp_list_of_ranges:=[];

>

> for i from 1 to n do

> T:=list_of_intervals[i];

>

> # Bisection of current interval

>

> t_mid:=T[1] + 0.5*(T[2]-T[1]);

> T1:=construct(T[1],t_mid);

> T2:=construct(t_mid,T[2]);

>

>
# Optional 4th parameter: Iteration counter iter_counter

> if (nargs >= 4) and (type(args[4],nonnegint)) then

> iter_counter:=args[4];

> else

> iter_counter:=1;

> fi;

>

> if iter_counter <= N then

> F1:=compute_naive_interval_range(f,T1);

> F2:=compute_naive_interval_range(f,T2);

> else

> F1:=compute_combined_range(f,T1);

> F2:=compute_combined_range(f,T2);

> fi;

>

> temp_list_of_intervals:=[op(temp_list_of_intervals),T1,T2];

> temp_list_of_ranges:=[op(temp_list_of_ranges),F1,F2];

> od:

>

> list_of_intervals:=[op(temp_list_of_intervals)];

> list_of_ranges:=[op(temp_list_of_ranges)];

>

> temp:=list_of_ranges[1];

> for i from 2 to nops(eval(list_of_ranges)) do

> temp:=&Convex_Hull(temp,list_of_ranges[i]);

> od;

>

> RETURN(temp);

> end:

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

> subdivide_adaptive:=proc(f::{function,procedure},

> list_of_intervals::evaln,

> list_of_ranges::evaln)

> local N, n, iter_counter, temp_max, temp_min,i, T, T1, T2,

> F1, F2, temp, t_mid, temp_list_of_intervals,

> temp_list_of_ranges, eps;

>

> N:=3;

>

> if not(type(eval(list_of_intervals),listlist)) then

> lprint(`Second parameter must be list of intervals!`);

> T:=readstat(`Enter start interval: `);

> T:=construct(T[1],T[2]);

> list_of_intervals:=[T];

> list_of_ranges:=[compute_naive_interval_range(f,T)];

> fi;

>

> n:=nops(eval(list_of_intervals));

>

> if nops(eval(list_of_ranges)) < n then

> ERROR(`Second list is too short!`);

> fi;

>

> temp_min:=min(seq(list_of_ranges[i][1],i=1..n));

> temp_max:=max(seq(list_of_ranges[i][2],i=1..n));

>

> eps:=0.1*width(construct(temp_min,temp_max));

>

> temp_list_of_intervals:=[];

> temp_list_of_ranges:=[];

>

> for i from 1 to n do

> T:=list_of_intervals[i];

>
# If T satisfies the selection criterion, T will be subdivided into halves

>

> if (list_of_ranges[i][1] <= temp_min + eps) or

> (list_of_ranges[i][2] >= temp_max -eps ) then

>

> # Bisection of current interval

>

> t_mid:=T[1] + 0.5*(T[2] - T[1]);

> T1:=construct(T[1],t_mid);

> T2:=construct(t_mid,T[2]);

>

> if (nargs >= 4) and (type(args[4],nonnegint)) then

> iter_counter:=args[4];

> else

> iter_counter:=1;

> fi;

>

> if iter_counter <= N then

> F1:=compute_naive_interval_range(f,T1);

> F2:=compute_naive_interval_range(f,T2);

> else

> F1:=compute_combined_range(f,T1);

> F2:=compute_combined_range(f,T2);

> fi;

>

> temp_list_of_intervals:=[op(temp_list_of_intervals),

> T1,T2];

> temp_list_of_ranges:=[op(temp_list_of_ranges),F1,F2];

> else

> temp_list_of_intervals:=[op(temp_list_of_intervals),T];

> temp_list_of_ranges:=[op(temp_list_of_ranges),

> list_of_ranges[i]];

>

> fi;

> od:

>

> list_of_intervals:=[op(temp_list_of_intervals)];

> list_of_ranges:=[op(temp_list_of_ranges)];

>

> temp:=list_of_ranges[1];

> for i from 2 to nops(eval(list_of_ranges)) do

> temp:=&Convex_Hull(temp,list_of_ranges[i]);

> od;

>

> RETURN(temp);

> end:

-----------------------------------------------------------------------------

PROCEDURE INTERVAL_LIST_PLOT

> interval_list_plot:=proc(interval_list::list(list(numeric)),

> range_list::list(list(numeric)))

> local n, i, T, F, c, opts, p;

>

> if nops(interval_list) > nops(range_list) then

> ERROR(`Second list is too short!`);

> fi;

>

> opts:=[args[3..nargs]];

> n:=nops(interval_list);

> for i from 1 to n do

> T:=interval_list[i];

> F:=range_list[i];

> c:=plottools[rectangle]([T[1],F[2]],[T[2],F[1]]);

> p[i]:=display(c,op(opts));

> od;

>

> RETURN(display([seq(p[i],i=1..n)]));

> end:

-----------------------------------------------------------------------------

PROCEDURE COMPUTE_RANGE

> compute_range:=proc(f::{function,procedure},xstart,iterationsteps::nonnegint)

> global list_of_intervals, list_of_ranges, max_iter_counter, q,

> r, iter_counter;

> local T, F, S, xrange, i, m, nx, h, arg_list, conv_linear,

> conv_quadratic, div_adaptive, c_list, j,x;

>

> # Initialization of global variables

>

> x:='x';

>

> if type(simplify(f(x)),numeric) or (simplify(diff(f(x),x)) = 0) then

> lprint(`Function is constant = `, simplify(f(x)));

> RETURN();

> fi;

>

> if type(xstart,{interval,numeric..numeric}) then

> T:=construct(op(1,xstart),op(2,xstart));

> elif type(xstart,{identical('X'),identical('x')} = {interval,numeric..numeric}) then

> T:=construct(op(1,rhs(xstart)),op(2,rhs(xstart)));

> else

> ERROR(`Second argument must be an interval or a range.`);

> fi;

>

> if not(type(compute_naive_interval_range(f,T),interval)) then

> ERROR(`First argument must be function in one real variable with real values.`);

> fi;

>

> list_of_intervals:=[T];

>

> max_iter_counter:=iterationsteps;

> nx:=1;

> conv_linear:=false;

> conv_quadratic:=false;

> div_adaptive:=false;

>

> arg_list:=[args[4..nargs]];

>

> for i from 1 to nops(arg_list) do

> if type(arg_list[i],identical('Nx') = nonnegint) then

> nx:=rhs(arg_list[i]);

> elif arg_list[i] = 'linear' then

> conv_linear:=true;

> elif (arg_list[i] = 'quadratic') or (arg_list[i] = 'quadratisch') then

> conv_quadratic:=true;

> elif (arg_list[i] = 'adaptive') or (arg_list[i] = 'adaptiv') then

> div_adaptive:=true;

> elif type(arg_list[i],identical('colorlist') = list) then

> c_list:=rhs(arg_list[i]);

> fi;

> od;

>

> if conv_quadratic then

> F:=compute_combined_range(f,T);

> lprint(`Start range enclosure = `,F);

> list_of_ranges:=[F];

> r[0]:=F;

> else

> F:=compute_naive_interval_range(f,T);

> lprint(`Start range enclosure = `,F);

> list_of_ranges:=[F];

> r[0]:=F;

> fi;

>

> # Subdivision of Start interval in Nx subparts,

> # given the optional parameter Nx

>

> if nx > 1 then

> h:=width(T)/nx;

> list_of_intervals:=[seq([T[1] + (i-1)*h, T[1] + i*h],i=1..nx)];

>

> if conv_quadratic then

> list_of_ranges:=[seq(compute_combined_range(f,S),S=list_of_intervals)];

> else

> list_of_ranges:=[seq(compute_naive_interval_range(f,S),S=list_of_intervals)];

> fi;

>

> r[0]:=[min(seq(list_of_ranges[i][1],i=1..nx)),

> max(seq(list_of_ranges[i][2],i=1..nx))];

> printf(`%s %d %s = \n[%.10g,%.10g]\n`,

> `Range enclosure after subdivision into `,nx,

> ` intervals`,r[0][1],r[0][2]);

>

> fi;

>

> xrange:=T[1]..T[2];

>

> # Colour list for graphical representation

>

> if c_list = 'c_list' then

> c_list:=[blue,red];

> fi;

>

> q[0]:=interval_list_plot(list_of_intervals,list_of_ranges,color=c_list[1],

> style=line,thickness=2);

> x:='x';

> q[max_iter_counter+1]:=plot(f(x),x=xrange,color=black,

> numpoints=500,linestyle=1);

>

> # Iteration

>

> iter_counter:=0;

>

> for i from 1 to max_iter_counter do

> iter_counter:=iter_counter + 1;

>

> if conv_linear then

> if div_adaptive then

> r[i]:=subdivide_adaptive(f,list_of_intervals,list_of_ranges);

> else

> r[i]:=subdivide_equidistant(f,list_of_intervals,list_of_ranges);

> fi;

> elif conv_quadratic then

> if div_adaptive then

> r[i]:=subdivide_adaptive(f,list_of_intervals,list_of_ranges,100);

> else

> r[i]:=subdivide_equidistant(f,list_of_intervals,list_of_ranges,100);

> fi;

> else

> if div_adaptive then

> r[i]:=subdivide_adaptive(f,list_of_intervals,list_of_ranges,i);

> else

> r[i]:=subdivide_equidistant(f,list_of_intervals,list_of_ranges,i);

> fi;

>

> j:= (i mod nops(c_list)) + 1;

> if i < max_iter_counter then

> q[i]:=interval_list_plot(list_of_intervals,list_of_ranges,color=c_list[j],

> style=line,thickness=2);

> else

> q[i]:=interval_list_plot(list_of_intervals,list_of_ranges,color=yellow,style=patch);

> fi;

> od;

>

> m:=max_iter_counter;

> printf(`%s %s %s %d = %s[%.10g,%.10g]\n`,

> Range_enclosure ,after,step,m,` `,r[m][1],r[m][2]);

>

> # Graphical representation of the last three steps of the iteration

>

> if m <= 1 then

> display([seq(q[i],i=0..(m+1))],axes=framed,
title=`Iterative range enclosure of f`,titlefont=[TIMES,BOLD,12]);

> else

> display([seq(q[i],i=(m-2)..m),q[m+1]],axes=framed,
title=`Iterative range enclosure of f`,titlefont=[TIMES,BOLD,12]);

> fi;

> end:

-----------------------------------------------------------------------------

PROCEDURE MAX_ABS_ERROR

> max_abs_error:=proc(list_of_ranges::list(interval))

> local index_of_max, index_of_min,temp_max,temp_min, n, i, new_list;

> index_of_max:=1;

> index_of_min:=1;

> n:=nops(list_of_ranges);

> if n=0 or (n=1 and list_of_ranges[1]=[]) then

> ERROR(`Input parameter must contain at least one real interval. `);

> fi;

>

> temp_max:=list_of_ranges[1][2];

> temp_min:=list_of_ranges[1][1];

>

> for i from 1 to n do

> if list_of_ranges[i][2] > temp_max then

> temp_max:=list_of_ranges[i][2];

> index_of_max:=i;

> fi;

> if list_of_ranges[i][1] < temp_min then

> temp_min:=list_of_ranges[i][1];

> index_of_min:=i;

> fi;

> od;

> temp_max:=list_of_ranges[index_of_max][1];

> temp_min:=list_of_ranges[index_of_min][2];

>

> if temp_min > temp_max then

> RETURN(max(seq(width(list_of_ranges[i]),i=1..n)));

> else

> new_list:=[];

> for i from 1 to n do

> if is_in(list_of_ranges[i],[temp_min,temp_max]) then

> next;

> else

> new_list:=[op(new_list),list_of_ranges[i]];

> fi;

> od;

> n:=nops(new_list);

> lprint(`Number of intervals: `,n);

> RETURN(max(seq(width(new_list[i]),i=1..nops(new_list))));

> fi;

> end:

>

###########################################################################

FUNCTIONS OF TWO REAL VARIABLES

-----------------------------------------------------------------------------

PROCEDURE CNI_RANGE3D

Interval evaluation of a function f: R^2 -> R.

> cni_range3d:=proc(f,T,S)

> local F;

> F:=inapply(f(t,s),t,s);

> RETURN(F(T,S));

> end:

The procedure subdivide_equi3d computes a range enclosure for a function of two real variables by
subdivision of the intervals in list_of_intervals2d in two intervals each. If k is odd, subdivision is
performed in x direction, if k is even, subdivision is performed in y direction.

An element of the global variable list_of_intervals2d should be a two dimensional interval, defined as
a list of two intervals, i.e. XY:=[[x1,x2],[y1,y2]].
The variable list_of_ranges contains the range enclosures
for the corresponding (sub-)intervals in list_of_intervals2d.

> subdivide_equi3d:=proc(f,

> list_of_intervals2d::evaln,

> list_of_ranges::evaln,

> k::nonnegint)

> local n, i, temp_loi, temp_lor, TS, T, S, T1, T2, S1, S2, F1, F2,

> F3, F4, t_mid, s_mid, temp;

>

> if not(type(eval(list_of_intervals2d),listlist)) then

> TS:=readstat(`Enter two-dimensional start interval [[x1,x2],[y1,y2]]:`);

> T:=construct(TS[1][1],TS[1][2]);

> S:=construct(TS[2][1],TS[2][2]);

> list_of_intervals2d:=[[T,S]];

> list_of_ranges:=[cni_range3d(f,T,S)];

> fi;

>

> n:=nops(eval(list_of_ranges));

> temp_loi:=[];

> temp_lor:=[];

> for i from 1 to n do

> TS:=list_of_intervals2d[i];

> T:=TS[1];

> S:=TS[2];

> if type(k,odd) then

> t_mid:=T[1] + 0.5*(T[2] - T[1]);

> T1:=construct(T[1],t_mid);

> T2:=construct(t_mid,T[2]);

> F1:=cni_range3d(f,T1,S);

> F2:=cni_range3d(f,T2,S);

> temp_loi:=[op(temp_loi),[T1,S],[T2,S]];

> temp_lor:=[op(temp_lor),F1,F2];

> else

> s_mid:=S[1] + 0.5*(S[2] - S[1]);

> S1:=construct(S[1],s_mid);

> S2:=construct(s_mid,S[2]);

> F1:=cni_range3d(f,T,S1);

> F2:=cni_range3d(f,T,S2);

> temp_loi:=[op(temp_loi),[T,S1],[T,S2]];

> temp_lor:=[op(temp_lor),F1,F2];

> fi;

> od;

>

> list_of_intervals2d:=[op(temp_loi)];

> list_of_ranges:=[op(temp_lor)];

>

> temp:=list_of_ranges[1];

> for i from 2 to nops(eval(list_of_ranges)) do

> temp:=&Convex_Hull(temp,list_of_ranges[i]);

> od;

>

> RETURN(temp);

> end:

>

-----------------------------------------------------------------------------

PROCEDURE INTERVAL_LIST_PLOT3D

Graphical representation of a single iteration step

> interval_list_plot3d:=proc(interval_list,range_list)

> local n, i, j, TS, T, S, F, c, c_list, p, opts, c_out;

>

> opts:=[args[3..nargs]];

>

> for i from 1 to nops(opts) do

> if type(opts[i],identical('colorlist') = list) then

> c_list:=rhs(opts[i]);

> elif type(opts[i],identical('cutout') = rational) then

> c_out:=rhs(opts[i]);

> fi;

> od;

>

> if c_list = 'c_list' then

> c_list:=[blue,red,green,magenta];

> fi;

>

> if c_out = 'c_out' then

> c_out := 8/10;

> fi;

>

> n:=nops(range_list);

>

> for i from 1 to n do

> TS:=interval_list[i];

> T:=TS[1];

> S:=TS[2];

>

> F:=range_list[i];

> c:=plottools[cuboid]([T[1],S[1],F[1]],[T[2],S[2],F[2]]);

>

> # Different colours are applied for individual intervals

> p[i]:=display(plottools[cutout](c,c_out),

> color=c_list[(i mod nops(c_list)) + 1]);

> # p[i]:=display([c],color=c_list[(i mod nops(c_list)) + 1]);

> od;

>

> RETURN(display([seq(p[i],i=1..n)]));

> end:

>

-----------------------------------------------------------------------------

PROCEDURE COMPUTE_RANGE3D

> compute_range3d:=proc(f::{function,procedure},X,Y,iterationsteps::nonnegint)

> global list_of_intervals2d,list_of_ranges,q, r;

> local TS, T, S, iter_counter, max_iter_counter, i, j, m, nx, ny, hx,

> hy, lx, ly, opts, c_list,c_out,n,opts_new,x,y;

>

> x:='x';

> y:='y';

>

> if type(simplify(f(x,y)),numeric) then

> lprint(`Funktion ist konstant = `, simplify(f(x,y)));

> RETURN();

> fi;

>

> if type(X,{interval,numeric..numeric}) then

> T:=construct(op(1,X),op(2,X));

> elif type(X,{identical('X'),identical('x')} =

> {interval,numeric..numeric}) then

> T:=construct(op(1,rhs(X)),op(2,rhs(X)));

> else ERROR(`Second argument must be an interval or a range`);

> fi;

>

> if type(Y,{interval,numeric..numeric}) then

> S:=construct(op(1,Y),op(2,Y));

> elif type(Y,{identical('Y'),identical('y')} =

> {interval,numeric..numeric}) then

> S:=construct(op(1,rhs(Y)),op(2,rhs(Y)));

> else ERROR(`Third argument must be an interval or a range`);

> fi;

>

> if not(type(cni_range3d(f,T,S),interval)) then

> ERROR(`First argument must be real function of two real variables!`);

> fi;

>

> lprint(`Start range enclosure = `,cni_range3d(f,T,S));

>

> # Initialization of global variables

>

> max_iter_counter:=iterationsteps;

> nx:=1;

> ny:=1;

>

> opts:=[args[5..nargs]];

>

> n:=nops(opts);

> opts_new:=[];

> for i from 1 to n do

> if type(opts[i],identical('Nx') = nonnegint) then

> nx:=rhs(opts[i]);

> elif type(opts[i],identical('Ny') = nonnegint) then

> ny:=rhs(opts[i]);

> elif type(opts[i],identical('colorlist') = list) then

> c_list:=rhs(opts[i]);

> elif type(opts[i],identical('cutout') = numeric) then

> c_out:=rhs(opts[i]);

> else opts_new:=[op(opts_new),opts[i]]

> fi;

> od;

>

> q[max_iter_counter+1]:=plot3d(f(x,y),x=T[1]..T[2],y=S[1]..S[2]);

>

> # Subdivision of the start interval into in nx*ny subparts,

> # given an optional parameter [nx, ny].

>

> if nx > 1 then

> hx:=width(T)/nx;

> lx:=[seq([T[1] + (i-1)*hx, T[1] + i*hx],i=1..nx)];

> else

> lx:=[T];

> fi;

> if ny > 1 then

> hy:=width(S)/ny;

> ly:=[seq([S[1] + (i-1)*hy, S[1] + i*hy],i=1..ny)];

> else

> ly:=[S];

> fi;

> list_of_intervals2d:=[seq(seq([lx[i],ly[j]],j=1..ny),i=1..nx)];

>

> list_of_ranges:=[seq(cni_range3d(f,TS[1],TS[2]),TS=list_of_intervals2d)];

>

> r[0]:=[min(seq(list_of_ranges[i][1],i=1..(nx*ny))),

> max(seq(list_of_ranges[i][2],i=1..(nx*ny)))];

> if nx >1 or ny > 1 then

> printf(`%s %d %s = \n[%.10g,%.10g]\n`, `Range enclosure after subdivision into`,nx*ny, `intervals`,r[0][1],r[0][2]);

> fi;

>

> if c_list = 'c_list' then

> c_list:=[blue,red,green,magenta];

> fi;

>

> if c_out = 'c_out' or not(type(c_out,rational) and (0 <= c_out) and

> (c_out <= 1)) then

> c_out:=8/10;

> fi;

>

> q[0]:=interval_list_plot3d(list_of_intervals2d,list_of_ranges,colorlist=c_list,cutout=c_out);

>

> iter_counter:=0;

>

> for i from 1 to max_iter_counter do

> iter_counter:=iter_counter + 1;

>

> r[i]:=subdivide_equi3d(f,list_of_intervals2d,list_of_ranges,iter_counter);

> # lprint(`Range enclosure after Iteration step`,i,` = `,r[i]);

printf(`%s %s %s %s %d = %s[%.10g,%.10g]\n`,

> Range,enclosure,after,step,i,` `,r[i][1],r[i][2]);

>
q[i]:=interval_list_plot3d(list_of_intervals2d,list_of_ranges,colorlist=c_list,cutout=c_out);

> od;

>

> m:=max_iter_counter;

> display([q[m],q[m+1]],op(opts_new),axes=framed);

>

> end:

>

###########################################################################

EXTENDED INTERVAL NEWTON METHOD

>

-----------------------------------------------------------------------------

PROCEDURE NEWTON

Newton-Iteration

> newton:=proc(f,df,xold::interval,eps::float,xunique::boolean)

> global zeros,infos,N,iter_counter;

> local c,ynew,xnew,i,unique,j,insert;

>

> if not(type(xold[1],numeric) and type(xold[2],numeric)) then

> RETURN(`Start_interval_must_be_real!`);

> fi:

>

> # Cancel Iteration if no zero exists in start interval

>

> if not(is_in(0.,f(xold))) then

> RETURN();

> fi:

>

> # Increase number of iteration steps by 1

>

> iter_counter:=iter_counter + 1;

>

> # Extended interval Newton step

>

> c:=midpoint(xold);

>

> # c:=mid(xold);

>

> ynew:=ext_int_div(f(c),df(xold));

>

> if (nops([ynew]) >= 2) then

> xnew[1]:=(c &- op(1,[ynew])) &intersect xold;

> xnew[2]:=(c &- op(2,[ynew])) &intersect xold;

> unique:=false;

> else

> ynew:=(c &- ynew);

> xnew[1]:=ynew &intersect xold;

>

> if xnew[1] = [] then RETURN() fi;

>

> # Strict inclusion implies existence and uniqueness

> # of zero in xold

>

> unique:=xunique or

> ((xold[1] < xnew[1][1]) and (xnew[1][2] < xold[2]));

> xnew[2]:=[];

> fi:

>

> # Bisection, if no reduction of interval is achieved

>

> if (xnew[1] = xold) or (xnew[2] = xold) then

> c:=c[2];

>

> if (xold[1] < c) and (c < xold[2]) then

> xnew[1]:=[xold[1],c];

> xnew[2]:=[c,xold[2]];

> else

> N:=N+1;

> zeros[N]:=xold;

> infos[N]:=unique;

> RETURN();

> fi;

>

> fi:

>

> # Output of computed enclosures or further iteration steps

>

> for i from 1 to 2 do

> if (xnew[i] <> []) then

> if (rel_diam(xnew[i]) < eps) then

> if (is_in(0.,f(xnew[i]))) then

> if unique=false then

> insert:=true;

>

> # Combine enclosure intervals with equal end points

>

> for j from 1 to N do

> if ((xnew[i] &intersect zeros[j]) <> []) and

> (infos[j] = false) then

> zeros[j]:=zeros[j] &union xnew[i];

> insert:=false;

> break;

> fi;

> od;

>

> if insert then

> N:=N + 1;

> zeros[N]:=xnew[i];

> infos[N]:=unique;

> fi;

> else

> N:=N + 1;

> zeros[N]:=xnew[i];

> infos[N]:=unique;

> fi;

> fi:

> else

> newton(f,df,xnew[i],eps,unique);

> fi:

> od:

> RETURN();

> end:

-----------------------------------------------------------------------------
PROCEDURE COMPUTE_ALL_ZEROS

>

> compute_all_zeros:=proc(g::{function,procedure},xstart::interval,eps)

> global zeros,infos,N,iter_counter;

> local f,df,xold,unique,i,x,x1,x2,n;

>

> # Convert g and g' to interval functions

>

> x:='x';

>

> if type(simplify(g(x)),numeric) or (simplify(diff(g(x),x)) = 0) then

> lprint(`Function is constant = `, simplify(g(x)));

> RETURN();

> fi;

>

> f:=inapply(g(x),x);

> df:=inapply(diff(g(x),x),x);

>

> if not(type(f(xstart),interval)) then

> ERROR(`First argument must be real function of one real variable.`);

> fi;

>

> if not(type(xstart[1],numeric) and type(xstart[2],numeric)) then

> ERROR(`Start_interval_must_be_real.`);

> fi:

>

> x1:=xstart[1];

> x2:=xstart[2];

>

> # Determine precision of computation

>

> if (nargs > 3) and type(args[4],nonnegint) then

> Digits:=max(args[4],10,Digits);

> fi;

>

> while not(is_in(Evalf(x1),x1) and is_in(Evalf(x2),x2))

> and (Digits <= 100) do

> Digits:=eval(Digits) + 5;

> od;

>

> Digits:= max(ceil(abs(Evalf(log10(eps)))) + 5,

> Digits);

> lprint(`Digits = `,Digits);

>

> lprint(` `);

>

> if width(xstart)=0 then

> xold:=[rd(x1),ru(x2)];

> Digits:=Digits+2;

> else xold:=xstart;

> fi;

>

> # Initialization

>

> zeros:=table();

> infos:=table();

> N:=0;

> iter_counter:=0;

> unique:=false;

>

> # Perform Iteration

>

> newton(f,df,xold,Evalf(eps),unique);

>

> # Output

>

> if N = 0 then

> lprint(`No zero contained in start interval.`);

> RETURN();

> fi;

>

> for i from 1 to N do

> print(eval(zeros[i]));

> if infos[i]=true then

> print(`contains exactly one zero.`);

> else

> print(`potential enclosure of zero`);

> fi;

> od;

> lprint(` `);

> lprint(`Number of enclosures of zeros: `,N);

> lprint(`Number of Iterations steps: `,iter_counter);

> end:

>

-----------------------------------------------------------------------------

PROCEDURE NEWTON_PLOT

Graphical representation of one interval newton step.

> newton_plot:=proc(g,f,xold::interval,ynew::list,xnew::list)

> local

> xrange,yrange,p,hy,h1,h2,i,intxold,intxnew,n,punkt,punkt1,punkt2,punkt3,cx,x,c;

>

> lprint(`xold=`,xold);

> lprint(`xnew1=`,xnew[1]);

>

> if (xnew[2] <> []) then

> lprint(`xnew2=`,xnew[2]);

> fi;

>

> xrange:=(xold[1]-(width(xold)/10.))..(xold[2]+(width(xold)/10.));

>

> yrange:=(f(xold)[1]-width(f(xold))/10.)..(f(xold)[2]+width(f(xold))/10.);

>

> x:='x';

> p[1]:=plot(g(x),x=xrange,color=black,numpoints=200):

> hy:=width(f(xold));

> h1:=hy/25.;

> h2:=hy/50.;

>

> # Representation of old interval - yellow

> intxold:=plottools[rectangle]([xold[1],0.],[xold[2],-h2]);

> p[2]:=PLOT(intxold,COLOR(RGB,1,1,0));

>

> # Boundaries of old interval

>

> p[3]:=plot([[xold[1],-h1],[xold[1],h1]],color=black,thickness=4):

>

> p[4]:=plot([[xold[2],-h1],[xold[2],h1]],color=black,thickness=4):

>

> # Labelling for old interval

>

> p[5]:=textplot([midpoint(xold)[1],-h1,'xold'],align=BELOW);

>

> # Representation of computed new intervals

>

> n:=6;

>

> for i from 1 to 2 do

> if xnew[i] <> [] then

>

> intxnew:=plottools[rectangle]([xnew[i][1],h2],[xnew[i][2],0.]);

> p[n]:=PLOT(intxnew,COLOR(RGB,0,1,0));

> n:=n+1:

>

> p[n]:=plot([[xnew[i][1],-h1],[xnew[i][1],h1]],thickness=4):

> n:=n+1:

>

> p[n]:=plot([[xnew[i][2],-h1],[xnew[i][2],h1]],thickness=4):

> n:=n+1:

> if i=1 then

>

> p[n]:=textplot([midpoint(xnew[i])[1],h1,'xnew1'],align=ABOVE);

> elif i=2 then

>

> p[n]:=textplot([midpoint(xnew[i])[1],h1,'xnew2'],align=ABOVE);

> fi:

> n:=n+1;

> fi:

> od:

>

> # representation of tangent - red.

>

> c:=midpoint(xold)[2];

> punkt[1]:=[c,g(c)];

>

> # Determination of intersection points

>

> if nops(ynew)=2 then

> if is_in(-infinity,ynew[1]) then

> punkt[2]:=[ynew[1][2],0];

> punkt[3]:=[ynew[2][1],0];

> else

> punkt[2]:=[ynew[1][1],0];

> punkt[3]:=[ynew[2][2],0];

> fi;

> else

> if ynew[1] = [-infinity,infinity] then

>

> p:=display([seq(p[i],i=1..n-1)],axes=framed,xtickmarks=5,

>

> view=[xrange,yrange],font=[TIMES,ITALIC,10],axesfont=[TIMES,ROMAN,10],

> scaling=unconstrained);

>

> RETURN(p);

>

> elif is_in(-infinity,ynew[1]) then

> punkt[2]:=[ynew[1][2],0];

> punkt[3]:=[0.,g(c)];

> elif is_in(infinity,ynew[1]) then

> punkt[2]:=[ynew[1][1],0];

> punkt[3]:=[0.,g(c)];

> else

> punkt[2]:=[ynew[1][1],0];

> punkt[3]:=[ynew[1][2],0];

> fi;

>

> point(punkt1,punkt[1]):

> point(punkt2,punkt[2]):

> point(punkt3,punkt[3]);

> line(g1,[punkt1,punkt2]);

> line(g2,[punkt1,punkt3]);

>

> p[n]:=draw([g1],view=[xrange,yrange],color=red,linestyle=3);

>

> n:=n+1;

>

> p[n]:=draw([g2],view=[xrange,yrange],color=red,linestyle=3);

>

> # Representation of all graphical elements

>

> p:=display([seq(p[i],i=1..n)],axes=framed,xtickmarks=5,view=[xrange,yrange],

> font=[TIMES,ITALIC,10],axesfont=[TIMES,ROMAN,10],scaling=unconstrained);

>

> RETURN(p);

> end:

>

-----------------------------------------------------------------------------

PROCEDURE NEWTON_WITH_PLOT

Newton-Iteration with graphical representation

> newton_with_plot:=proc(g,f,df,xold::interval,eps::float,xunique::boolean)

> local c,ynew,xnew,i,j,insert,unique;

> global zeros,infos,N,iter_counter,max_iter_counter,p,dx;

>

> if not(type(xold[1],numeric) and type(xold[2],numeric)) then

> RETURN(`Start interval_must_be_real.`);

> fi:

>

> if not(is_in(0.,f(xold))) then

> RETURN();

> fi:

>

> # Increase number of iteration steps by 1

>

> iter_counter:=iter_counter + 1;

>

> # Extended interval Newton Step

>

> c:=midpoint(xold);

>

> ynew:=[ext_int_div(f(c),df(xold))];

> if (nops(ynew) >= 2) then

> ynew[1]:=c &- ynew[1];

> xnew[1]:=ynew[1] &intersect xold;

>

> ynew[2]:=c &- ynew[2];

> xnew[2]:=ynew[2] &intersect xold;

>

> unique:=false;

> else

> ynew[1]:=c &- ynew[1];

> xnew[1]:=ynew[1] &intersect xold;

>

> if xnew[1] = [] then RETURN() fi;

>

> # Strict inclusion implies existence and uniqueness

> # of zero in xold

>
unique:=xunique or

> ((xold[1] < xnew[1][1]) and (xnew[1][1] < xold[2]));

>

> xnew[2]:=[];

> fi:

>

> # Bisection, if no reduction of interval is achieved

> if (xnew[1] = xold) or (xnew[2] = xold) then

> c:=c[2];

> if (xold[1] < c) and (c < xold[2]) then

> xnew[1]:=[xold[1],c];

> xnew[2]:=[c,xold[2]];

> else

> N:=N+1;

> zeros[N]:=xold;

> infos[N]:=unique;

> RETURN();

> fi;

>

> fi:

>

> # Representation of Newton step

>

> xnew:=[xnew[1],xnew[2]];

>

> lprint(`Iteration step `, iter_counter);

> p[iter_counter]:=newton_plot(g,f,xold,ynew,xnew);

> print(p[iter_counter]);

>

> # Output of computed enclosures or further iteration steps

>

> for i from 1 to 2 do

> if (xnew[i] <> []) then

> if (rel_diam(xnew[i]) < eps) or (iter_counter=max_iter_counter) then

> if (is_in(0.,f(xnew[i]))) then

> if unique=false then

> insert:=true;

>

> # Combine enclosure intervals with equal end points

>

> for j from 1 to N do

> if ((xnew[i] &intersect zeros[j]) <> []) and

> (infos[j] = false) then

> zeros[j]:=zeros[j] &union xnew[i];

> insert:=false;

> break;

> fi;

> od;

>

> if insert then

> N:=N + 1;

> zeros[N]:=xnew[i];

> infos[N]:=unique;

> fi;

> else

> N:=N + 1;

> zeros[N]:=xnew[i];

> infos[N]:=unique;

> fi;

> fi:

> else

> newton_with_plot(g,f,df,xnew[i],eps,unique);

> fi:

> od:

> RETURN();

> end:

>

-----------------------------------------------------------------------------

PROCEDURE COMPUTE_ALL_ZEROS_WITH_PLOT

Enclosure of the zeros of a continuously differentiable function g with graphical representation

> compute_all_zeros_with_plot:=proc(g::{function,procedure},xstart::interval,eps)

> global zeros,infos,N,iter_counter,max_iter_counter;

> local f,df,xold,unique,i,x,x1,x2;

>

> x:='x';

>

> if type(simplify(g(x)),numeric) or (simplify(diff(g(x),x)) = 0) then

> lprint(`Function is constant = `, simplify(g(x)));

> RETURN();

> fi;

>

> # Conversion of g and g' into interval functions

>

> if (nargs > 4) and type(args[5],nonnegint) then

> max_iter_counter:=args[5];

> else

> max_iter_counter:=readstat("Enter maximal number of iteration steps: ");

>

> fi;

>

> f:=inapply(g(x),x);

> df:=inapply(diff(g(x),x),x);

>

> if not(type(f(xstart),interval)) then

> ERROR(`First argument must be real function in one real variable.`);

> fi;

>

> x1:=xstart[1];

> x2:=xstart[2];

>

> # Determine precsion of computation

>

> if (nargs > 3) and type(args[4],nonnegint) then

> Digits:=max(args[4],10,Digits);

> fi;

>

> while not(is_in(Evalf(x1),x1) and is_in(Evalf(x2),x2))

> and (Digits <= 100) do

> Digits:=eval(Digits) + 5;

> od;

>

> Digits:= max(ceil(abs(Evalf(log10(eps)))) + 5,Digits);

>

> lprint(`Digits = `,Digits);

> lprint(` `);

>

> if width(xstart)=0 then

> xold:=[rd(x1),ru(x2)];

> Digits:=Digits+2;

> else xold:=xstart;

> fi;

>

> # Initialization

>

> zeros:=table();

> infos:=table();

> N:=0;

> iter_counter:=0;

> unique:=false;

>

> # Perform iteration

>

> newton_with_plot(g,f,df,xold,Evalf(eps),unique);

>

> # Output

>

> if N = 0 then

> lprint(`No zero contained in start interval.`);

> RETURN();

> fi;

>

> for i from 1 to N do

> print(eval(zeros[i]));

> if infos[i]=true then

> print(` contains exactly one zero.`);

> else

> print(` potential enclosure of zero`);

> fi;

> od;

> lprint(` `);

> lprint(`Number of enclosures of zeros: `,N);

> lprint(`Number of iteration steps: `,iter_counter);

> end:

>

###########################################################################

COMPLEX DISC ARITHMETIC

###########################################################################

Type definitions and aouxiliary procedures

Type definitions --> see init procedure

-----------------------------------------------------------------------------

Definition of `type/complex_disc`:

--> moved to init procedure!

-----------------------------------------------------------------------------

Definition of `type/complex_interval`:

--> moved to init procedure!

-----------------------------------------------------------------------------

PROCEDURE `&CABS`

This procedure computes a real interval enclosing the abolute value of a complex number
or a complex rectangular interval.

> `&cabs`:=proc(z)

> local w;

> if type(z,complex_interval) then

> w:=&sqrt( (&sqr(z[1])) &+ (&sqr(z[2])) );

> elif type(z,complex) then

> w:=&sqrt( (&sqr(Re(z))) &+ (&sqr(Im(z))) );

> else

> '&cabs(z)';

> fi;

> end:

>

-----------------------------------------------------------------------------

PROCEDURE COMPLEX_DISC_PLOT

Graphical representation of a complex disc

> complex_disc_plot:=proc(z::complex_disc)

> local opts, x, y, r;

> opts:=[args[2..nargs]];

> x:=z[1];

> y:=z[2];

> r:=z[3];

> plot([x + r*cos(t), y +r*sin(t), t=0..2*Pi],op(opts));

> end:

> complex_polynom_plot:=proc(p::polynom(complex,z),u::complex_disc)

> local opts, x, y, r;

>

> opts:=[args[3..nargs]];

> x:=u[1];

> y:=u[2];

> r:=u[3];

> complexplot(subs(z=x + I*y + r * (cos(t) + I*sin(t)),p),t=0..2*Pi,op(opts));

> end:

>

###########################################################################

BASIC ARITHMETIC OPERATIONS

-----------------------------------------------------------------------------

PROCEDURE `&CADD`

Centered addition of two Complex discs

> `&cadd`:=proc(u::{complex,complex_disc},v::{complex,complex_disc})

> local cx, cy, c, r;

>

> if type(u,complex_disc) and type(v,complex_disc) then

> cx:= u[1] &+ v[1];

> cy:= u[2] &+ v[2];

> c:=[mid(cx),mid(cy)];

> if u[3]=0 then

> r:=v[3];

> elif v[3]=0 then

> r:=u[3];

> else

> r:= u[3] &+ v[3];

> r:=r[2];

> fi;

>

> r:= ru(Evalf(r) + width(cx));

> r:= ru(Evalf(r) + width(cy));

> RETURN([op(c),r]);

> elif type(u,complex) and type(v,complex_disc) then

> [Re(u),Im(u),0] &cadd v;

> elif type(u,complex_disc) and type(v,complex) then

> u &cadd [Re(v),Im(v),0];

> else

> [Re(u),Im(u),0] &cadd [Re(v),Im(v),0];

> fi;

>

> end:

-----------------------------------------------------------------------------

PROCEDURE `&CSUB`

Centered subtraction of two Complex discs

> `&csub`:=proc(u::{complex,complex_disc},v::{complex,complex_disc})

> local cx, cy, c, r;

> if type(u,complex_disc) and type(v,complex_disc) then

> cx:=u[1] &- v[1];

> cy:=u[2] &- v[2];

> c:=[mid(cx),mid(cy)];

> if u[3] = 0 then

> r:=v[3];

> elif v[3]=0 then

> r:=u[3];

> else

> r:= u[3] &+ v[3];

> r:=r[2];

> fi;

>

> r:= ru(Evalf(r) + width(cx));

> r:= ru(Evalf(r) + width(cy));

> RETURN([op(c),r]);

> elif type(u,complex) and type(v,complex_disc) then

> [Re(u),Im(u),0] &csub v;

> elif type(u,complex_disc) and type(v,complex) then

> u &csub [Re(v),Im(v),0];

> else

> [Re(u),Im(u),0] &csub [Re(v),Im(v),0];

> fi;

> end:

-----------------------------------------------------------------------------

PROCEDURE `&CMULT`

Centered multiplication of two Complex discs

> `&cmult`:=proc(u::{complex,complex_disc},v::{complex,complex_disc})

> local cx, cy, c, r, r1, r2, r3;

>

> if type(u,complex_disc) and type(v,complex_disc) then

> cx:=(u[1] &* v[1]) &- (u[2] &* v[2]);

> cy:=(u[1] &* v[2]) &+ (u[2] &* v[1]);

> c:=[mid(cx),mid(cy)];

> if u[3]=0 then

> r:=(&cabs(u[1]+I*u[2])) &* v[3];

> elif v[3]=0 then

> r:=(&cabs(v[1]+I*v[2])) &* u[3];

> else

> r1:= (&cabs(u[1]+I*u[2])) &* v[3];

> r2:= (&cabs(v[1]+I*v[2])) &* u[3];

> r3:= u[3] &* v[3];

> r:= r1 &+ r2 &+ r3;

> fi;

> r:=r[2];

> r:= ru(Evalf(r) + width(cx));

> r:= ru(Evalf(r) + width(cy));

> RETURN([op(c),r]);

> elif type(u,complex) and type(v,complex_disc) then

> [Re(u),Im(u),0] &cmult v;

> elif type(u,complex_disc) and type(v,complex) then

> u &cmult [Re(v),Im(v),0];

> else

> [Re(u),Im(u),0] &cmult [Re(v),Im(v),0];

> fi;

> end:

-----------------------------------------------------------------------------

PROCEDURE `&CDIV`

Centered division of two Complex discs

> `&cdiv`:=proc(u::{complex,complex_disc},v::{complex,complex_disc})

> local c, cx, cy, r, nenner, a;

>

> if type(u,complex_disc) and type(v,complex_disc) then

> a:=&cabs(v[1]+I*v[2]);

> if (a[1] <= v[3]) then

> ERROR(`Denominator contains zero! `):

> elif (u[1] = 1) and (u[2]=0) and (u[3]=0) then

> nenner:= (&sqr(&cabs((v[1]+I*v[2])))) &- (&sqr(v[3]));

> if is_in(0,nenner) then

> ERROR(`Formula denominator contains zero! `):

> else

> cx:=v[1] &/ nenner;

> cy:=((-1)*v[2]) &/ nenner;

> r:=v[3] &/ nenner;

> c:=[mid(cx),mid(cy)];

> fi;

> else

> RETURN( u &cmult ( [1,0,0] &cdiv v));

>

> fi;

> r:=r[2];

> r:= ru(Evalf(r) + width(cx));

> r:= ru(Evalf(r) + width(cy));

> RETURN([op(c),r]);

> elif type(u,complex) and type(v,complex_disc) then

> [Re(u),Im(u),0] &cdiv v;

> elif type(u,complex_disc) and type(v,complex) then

> u &cdiv [Re(v),Im(v),0];

> elif type(u,complex) and type(v,complex) then

> [Re(u),Im(u),0] &cdiv [Re(v),Im(v),0];

> fi;

> end:

-----------------------------------------------------------------------------

PROCEDURE X0_START

x0_start:=proc(s1::interval,s2::interval)

> local f, F, xu, xo, s3, s, t;

> s:='s';

> t:='t';

> f:=(s,t)->(sqrt((1+s^2)*(1+t^2))*s*t - s^2*t^2) / (1+s^2+t^2);

> F:=inapply(f(s,t),s,t);

> xu:=F(s1,s2);

> xu:=xu[1];

> s3:=s1 &* s2;

> xo:=min(s1[2],s2[2],s3[2]);

> if xu <= xo then

> RETURN([xu,xo]);

> else

> ERROR(`Kann Startintervall nicht bestimmen`);

> fi;

>

> end:

-----------------------------------------------------------------------------

PROCEDURE `&CMULT_OPT`

Surface optimal multiplication of two complex discs.

> `&cmult_opt`:=proc(u::{complex,complex_disc},v::{complex,complex_disc})

> local c,cx, cy, r, x0, s1, s2, r1, r2, r3, z1, z2, z1z2, r1r2,

> r1z2, r2z1, P, dP,t, i, xold, m;

>

> if type(u,complex_disc) and type(v,complex_disc) then

> if (u[1]=0 and u[2]=0 and u[3]=0) or

> (v[1]=0 and v[2]=0 and v[3]=0) then

> RETURN([0,0,0]);

> # if one of the discs is zero, zero is returned

>
# u and/or w are circles with radius 0 ("a point")

>

> elif u[3]=0 then

> # Determine midpoint

> cx:=(u[1] &* v[1]) &- (u[2] &* v[2]);

> cy:=(u[1] &* v[2]) &+ (u[2] &* v[1]);

> c:=[mid(cx),mid(cy)];

>

> # Determine radius

>

> if v[3]=0 then

> r:=[0,0];

> else

> r:= (&cabs(u[1]+I*u[2])) &* v[3];

> fi;

>

> r:=r[2];

> r:=ru(Evalf(r) + width(cx));

> r:=ru(Evalf(r) + width(cy));

> RETURN([op(c),r]);

> elif v[3]=0 then

> # Determine midpoint

> cx:=(u[1] &* v[1]) &- (u[2] &* v[2]);

> cy:=(u[1] &* v[2]) &+ (u[2] &* v[1]);

> c:=[mid(cx),mid(cy)];

>

> # Determine radius

> r:=(&cabs(v[1]+I*v[2])) &* u[3];

>

> r:=r[2];

> r:=ru(Evalf(r) + width(cx));

> r:=ru(Evalf(r) + width(cy));

> RETURN([op(c),r]);

>

> # One factor is circle with midpoint 0.

>

> elif (u[1]=0) and (u[2]=0) then

> c:=[0,0];

> if (v[1]=0) and (v[2]=0) then

> r:= u[3] &* v[3];

> else

> r:= u[3] &* ((&cabs(v[1]+I*v[2])) &+ v[3]);

> fi;

>

> r:=r[2];

> RETURN([0,0,r]);

> elif (v[1]=0) and (v[2]=0) then

> c:=[0,0];

> r:=v[3] &* ((&cabs(u[1]+I*u[2])) &+ u[3]);

> r:=r[2];

> RETURN([0,0,r]);

> else

>

> s1:= u[3] &/ (&cabs(u[1]+I*u[2]));

> s2:= v[3] &/ (&cabs(v[1]+I*v[2]));

> z1:=u[1] + I*u[2];

> z2:=v[1] + I*v[2];

>

> x0:=x0_start(s1,s2);

> #lprint(x0);

> z1z2:=&intpower(&cabs(z1),2) &* &intpower(&cabs(z2),2);

> r1r2:=&intpower(u[3],2) &* &intpower(v[3],2);

>

> r1z2:=&intpower(u[3],2) &* &intpower(&cabs(z2),2);

> r2z1:=&intpower(v[3],2) &* &intpower(&cabs(z1),2);

>

> # Determine x0 using Newton method,

> # start interval = x0start

>

> t:='t';

> P:=(2 &* z1z2 &* (&intpower(t,3))) &+

> ((z1z2 &+ r1z2 &+ r2z1) &* (&intpower(t,2))) &- r1r2;

> #lprint(eval(subs(t=x0,P)));

> dP:=(6 &* z1z2 &* (&intpower(t,2))) &+ (2 &* (z1z2 &+ r1z2 &+ r2z1) &* t);

> #lprint(`dP=`,eval(subs(t=x0,dP)));

>

> for i from 1 to 10 do

> #lprint(x0);

> xold:=x0;

> m:=midpoint(xold);

> x0:=(m &-

> (eval(subs(t=m,P)) &/ eval(subs(t=xold,dP))) );

> x0:=x0 &intersect xold;

> if x0=xold then break; fi;

> od;

> #lprint(`x0=`,x0);

>

> # Determine midpoint

> cx:=((u[1] &* v[1]) &- (u[2] &* v[2])) &* (1 &+ x0);

> cy:=((u[1] &* v[2]) &+ (u[2] &* v[1])) &* (1 &+ x0);

> c:=[mid(cx),mid(cy)];

> #lprint(`c=`,c);

>

> # Determine radius

>

> r1:=3 &* z1z2 &* &intpower(x0,2);

> r2:=2 &* (z1z2 &+ r1z2 &+ r2z1) &* x0;

> r3:=r1z2 &+ r2z1 &+ r1r2;

> r:=&sqrt(r1 &+ r2 &+ r3);

> r:=r[2];

> r:=ru(Evalf(r) + width(cx));

> r:=ru(Evalf(r) + width(cy));

> RETURN([op(c),r]);

> fi;

> elif type(u,complex) and type(v,complex_disc) then

> [Re(u),Im(u),0] &cmult_opt v;

> elif type(u,complex_disc) and type(v,complex) then

> u &cmult_opt [Re(v),Im(v),0];

> else

> [Re(u),Im(u),0] &cmult_opt [Re(v),Im(v),0];

> fi;

> end:

-----------------------------------------------------------------------------

PROCEDURE `&CDIV_OPT`

Surface optimal division of two complex discs.

> `&cdiv_opt`:=proc(u::{complex,complex_disc},v::{complex,complex_disc})

> if type(u,complex_disc) and type(v,complex_disc) then

> RETURN( u &cmult_opt ([1,0,0] &cdiv v));

> elif type(u,complex) and type(v,complex_disc) then

> [Re(u),Im(u),0] &cdiv_opt v;

> elif type(u,complex_disc) and type(v,complex) then

> u &cdiv_opt [Re(v),Im(v),0];

> else

> [Re(u),Im(u),0] &cdiv_opt [Re(v),Im(v),0]

> fi;

> end:

>

###########################################################################

RANGE ENCLOSURES FOR COMPLEX POLYNOMIALS

-----------------------------------------------------------------------------

PROCEDURE HORNER_EVAL_CENT

Horner-Scheme with centered multiplication

> horner_eval_cent:=proc(p::polynom(complex,z),u::{complex,complex_disc})

> local n, Z, s, i, koeff, sx, sy;

> ### WARNING: degree(0,x) now returns -infinity

> n:=degree(p,z);

>

> # p ist constant polynomial

>

> if n=0 or n=-infinity then

> RETURN(p);

> fi;

>

> if type(u,complex) then

> Z:=[Re(u),Im(u),0];

> else

> Z:=u;

> fi;

> koeff:=coeff(p,z^n);

> s:=[Re(koeff),Im(koeff),0];

>

> for i from (n-1) to 0 by -1 do

> if i>0 then

> koeff:=coeff(p,z^i);

> if koeff=0 then

> s:=s &cmult Z;

> else

> s:=(s &cmult Z) &cadd [Re(koeff),Im(koeff),0];

> fi;

> elif (i=0) then

> ### WARNING: ldegree(0,x) now returns infinity

> if ldegree(p,z) = 0 then

> s:= (s &cmult Z) &cadd [Re(tcoeff(p)),Im(tcoeff(p)),0];

> else

> s:= s &cmult Z;

> fi;

>

> od;

>

> RETURN(s);

> end:

-----------------------------------------------------------------------------

PROCEDURE HORNER_EVAL_OPT

Horner-Scheme with surface optimal multiplication

> horner_eval_opt:=proc(p::polynom(complex,z),u::{complex,complex_disc})

> local n, Z, s, i, koeff;

> ### WARNING: degree(0,x) now returns -infinity

> n:=degree(p,z);

>

> if n=0 or n=-infinity then

> RETURN(p);

> fi;

>

> if type(u,complex) then

> Z:=[Re(u),Im(u),0];

> else

> Z:=u;

> fi;

>

> koeff:=coeff(p,z^n);

> s:=[Re(koeff),Im(koeff),0];

> for i from (n-1) to 0 by -1 do

> if i>0 then

> koeff:=coeff(p,z^i);

> if koeff=0 then

> s:=s &cmult_opt Z;

> else

> s:=(s &cmult_opt Z) &cadd [Re(koeff),Im(koeff),0];

> fi;

> elif (i=0) then

> ### WARNING: ldegree(0,x) now returns infinity

> if ldegree(p,z) = 0 then

> s:= (s &cmult_opt Z) &cadd [Re(tcoeff(p)),Im(tcoeff(p)),0];

> else

> s:= s &cmult_opt Z;

> fi;

>

> od;

> RETURN(s);

> end:

-----------------------------------------------------------------------------

PROCEDURE CENTRED_FORM_EVAL

Centered form using centered multiplication

> centred_form_eval:=proc(p::polynom(complex,z),

> u::{complex_disc,complex})

> local c, cx, cy, r, r1, r2, cu, pc, pcx, pcy, i, poly, n, Z;

>

> if type(u,complex) then

> Z:=[Re(u),Im(u),0];

> cu:=u;

> else

> Z:=u;

> cu:=u[1]+I*u[2];

> fi;

> ### WARNING: degree(0,x) now returns -infinity

> n:=degree(p);

>

> if n=0 or n=-infinity then

> RETURN(p);

> fi;

>

> # Determine midpont of p(Z)

>

> pc:=horner_eval_cent(p,cu);

> cx:=[op(1,pc[1] &- pc[3]),op(2,pc[1] &+ pc[3])];

> cy:=[op(1,pc[2] &- pc[3]),op(2,pc[2] &+ pc[3])];

> c:=[mid(cx),mid(cy)];

>

> # Determine radius of p(Z)

>

> r:=0;

> for i from 1 to n do

> poly:=diff(p,z$i);

> pc:=horner_eval_cent(poly,cu);

>

> if type(pc,complex_disc) then

> pcx:=[op(1,pc[1] &- pc[3]),op(2,pc[1] &+ pc[3])];

> pcy:=[op(1,pc[2] &- pc[3]),op(2,pc[2] &+ pc[3])];

> r1:= &cabs( [pcx,pcy] );

> else

> r1:=&cabs( pc);

> fi;

>

> if i>1 then

> r2:= (&intpower(Z[3],i));

> else

> r2:= Z[3];

> fi;

>

> r:= r &+ ((r1 &* r2) &/ (i!));

> od;

> r:=r[2];

> r:=ru(Evalf(r) + width(cx));

> r:= ru(Evalf(r) + width(cy));

> RETURN([op(c),r]);

> end:

>

###########################################################################

EXPONENTIAL FUNCTION FOR COMPLEX DISCS

-----------------------------------------------------------------------------
PROCEDURE CEXP

> cexp:=proc(u::{complex,complex_disc})

> local cx, cy, c, r, r1, r2;

>

> if type(u,complex_disc) then

>

> # Bestimmung des Mittelpunktes

>

> cx:=(&exp(u[1])) &* (&cos(u[2]));

> cy:=(&exp(u[1])) &* (&sin(u[2]));

> c:=[mid(cx),mid(cy)];

>

> # Bestimmung des Radius

>

> r1:=&cabs( [cx,cy] );

> r2:=(&exp(u[3])) &- 1;

> r:=r1 &* r2;

>

> r:=r[2];

> r:=ru(Evalf(r) + width(cx));

> r:=ru(Evalf(r) + width(cy));

>

> RETURN([op(c),r]);

> else

> cexp([Re(u),Im(u),0]);

> fi;

> end :

###########################################################################

> end module:

###########################################################################

end package intpakX

###########################################################################

intpakX v1.0 (c) Scientific Computing/Software Engineering Research Group,
Department of Mathematics, University of Wuppertal, Germany

Contact:
kraemer@math.uni-wuppertal.de
markus.grimmer@math.uni-wuppertal.de
http://www.math.uni-wuppertal.de/wrswt

###########################################################################

Disclaimer:

While every effort has been made to validate the solutions in this package,

Waterloo Maple Inc. and the contributors are not responsible for any errors

contained and are not liable for any damages resulting from the use of this

material.

###########################################################################

Warning, the name changecoords has been redefined

###########################################################################

You can create a Maple library with the following commands:

i) Tell Maple where to put the library:

> libname:="C:/mylib/interval", libname;
You can use the directory name of your choice in the above line. An empty directory is recommended.

ii) Create the new library:

> march('create',libname[1],100);

iii) Save the library

> savelib('intpakX');

###########################################################################
Now you can load the library as follows:

> restart;

> libname:="C:/mylib/interval", libname;

> with(intpakX);

Warning, the name changecoords has been redefined

Warning, the assigned name midpoint now has a global binding

Warning, the protected names ilog10, max and min have been redefined and unprotected

>