(* :Name: NumericalMath`IntervalRoots` *)

(* :Title: Interval Methods for Rootfinding *)

(* :Author: Jerry B. Keiper *)

(* :Summary:
This package implements bisection, secant, and Newton interval
methods for rootfinding
*)

(* :Context: NumericalMath`IntervalRoots` *) 

(* :Package Version: Mathematica 1.0 *)

(* :Copyright: Copyright 1990-1999,  Wolfram Research, Inc.  *)

(* :History:
	Written by Jerry B. Keiper, January 1994.
*)

(* :Keywords: interval arithmetic, interval analysis, rootfinding,
	bisection, secant, Newton
*)

(* :Source:
	R. E. Moore, Interval Arithmetic, Prentice-Hall, Englewood Cliffs,
	NJ, 1966.

	H. Ratschek and J. Rokne, New Computer Methods for Global
	Optimation, Ellis Horwood, Ltd., Chichester, England, 1988.

	E. Hansen, Global Optimization using Interval Analysis, Marcel
	Dekker, Inc., New York, 1992.

	H. Ratschek and J. Rokne, Computer Methods for the Range of
	Functions, Ellis Horwood, Ltd., Chichester, England, 1988.

	Y. Akyildiz and A. Hayes, Interval Analysis and Root
	Finding, Presented at the 1993 Mathematica Developer's
	Conference, Champaign, Illinois.
*)

(* :Mathematica Version: 3.0 *)

(* :Limitation:
*)

(* :Discussion:
	The three most important methods of rootfinding: bisection,
	secant, and Newton, also have there interval equivalents.
	This package implements those methods.

	The ordinary secant method is a way to find roots without
	finding derivatives.  The interval analog of the secant method
	requires the evaluation of (f[X]-f[x])/(X-x), which must be
	evaluated using Taylor series methods or the result will be
	the entire real line, because x is an element of the interval
	X.  (An alternative in the case of polynomial functions is to
	cancel the (X-x) from the denominator.)
*)

BeginPackage["NumericalMath`IntervalRoots`"];

IntervalBisection::usage =
"IntervalBisection[f, x, int, eps] uses the interval bisection method
to find subintervals of the given interval int of length less than eps
(possibly) containing roots of the expression f."

IntervalSecant::usage =
"IntervalSecant[f, x, int, eps] uses an interval secant method to find
subintervals of the given interval int of length less than eps (possibly)
containing roots of the expression f."

IntervalNewton::usage =
"IntervalNewton[f, x, int, eps] uses an interval Newton method to find
subintervals of the given interval int of length less than eps (possibly)
containing roots of the expression f."

Begin["NumericalMath`IntervalRoots`Private`"]; 

Options[IntervalBisection] =
Options[IntervalSecant] =
Options[IntervalNewton] =
	 {MaxRecursion -> 7, WorkingPrecision -> $MachinePrecision};

IntervalBisection::rec =
IntervalSecant::rec =
IntervalNewton::rec = "MaxRecursion exceeded."

IntervalBisection[f_, x_, intab_, eps_, opts___] :=
    Block[{int, n, p},
	{n, p} = {MaxRecursion, WorkingPrecision} /. {opts} /.
		 Options[IntervalBisection];
	int = Interval /@ (List @@ intab);
	IntervalUnion @@ Flatten[IntBis[f, x, #, eps, 1, n, p]& /@ int]
	];

IntBis[f_, x_, int_Interval, eps_, infd_, n_, p_] :=
    Block[{a, b, c},
	If[!IntervalMemberQ[f /. x -> int, 0], Return[{}]];
	a = int[[1,1]];
	b = int[[1,2]];
	If[b-a < eps, Return[int]];
	If[n == 0, Message[IntervalBisection::rec]; Return[int]];
	c = N[(a+b)/2, p];
	If[!NumberQ[c],
	    If[c === Indeterminate,
		(* both a and b are infinite *)
		c = 0,
		(* else one or the other of a and b is infinite *)
		c = If[c > 0, a + infd, b - infd]
		]
	    ];
	IntBis[f, x, #, eps, 2infd, n-1, p]& /@
		 {Interval[{a,c}], Interval[{c,b}]}
	];

IntervalSecant[f_, x_, int_Interval, eps_, opts___] :=
    Block[{jac, n, p, xm},
	{n, p} = {MaxRecursion, WorkingPrecision} /. {opts} /.
		 Options[IntervalBisection];
	jac = D[f /. x -> xm, xm] + (x - xm) D[f, {x, 2}]/2;
	IntervalUnion @@
	    Select[Flatten[(IntSec[f,jac,x,xm,#,eps,n,p])& /@ (List @@ int)],
		IntervalMemberQ[f /. x -> #, 0]&]
	];

IntSec[f_, jac_, x_, xm_, {a_, b_}, eps_, n_, p_] :=
    Block[{xmid, int = Interval[{a, b}]},
	If[b-a < eps, Return[int]];
	If[n == 0, Message[IntervalSecant::rec]; Return[int]];
	xmid = Interval[(a+b)/2];
	int = IntervalIntersection[int,
		xmid - N[f /. x -> xmid, p]/
			N[jac /. {x -> int, xm -> xmid}, p]];
	(IntSec[f, jac, x, xm, #, eps, n-1, p])& /@ (List @@ int)
	];

IntervalNewton[f_, x_, int_Interval, eps_, opts___] :=
    Block[{jac, n, p},
	{n, p} = {MaxRecursion, WorkingPrecision} /. {opts} /.
		 Options[IntervalBisection];
	jac = D[f, x];
	IntervalUnion @@
	    Select[Flatten[(IntNewt[f,jac,x,#,eps,n,p])& /@ (List @@ int)],
		IntervalMemberQ[f /. x -> #, 0]&]
	];

IntNewt[f_, jac_, x_, {a_, b_}, eps_, n_, p_] :=
    Block[{xmid, int = Interval[{a, b}]},
	If[b-a < eps, Return[int]];
	If[n == 0, Message[IntervalNewton::rec]; Return[int]];
	xmid = Interval[(a+b)/2];
	int = IntervalIntersection[int,
		xmid - N[f /. x -> xmid, p]/N[jac /. x -> int, p]];
	(IntNewt[f, jac, x, #, eps, n-1, p])& /@ (List @@ int)
	];

End[ ] (* "NumericalMath`IntervalRoots`Private`" *)

Protect[IntervalBisection, IntervalSecant, IntervalNewton];

EndPackage[ ] (* "NumericalMath`IntervalRoots`" *)

