Solving equations and inequalities using validated numeric methods

Adam Strzebonski, Wolfram Research Inc., USA
adams@wolfram.com

Abstract

Validated numeric methods can be used to significantly speed up exact solving of equations and inequalities and to find and represent exact solutions of problems that cannot be solved by known purely symbolic algorithms. Mathematica solvers use such methods to find and represent roots of univariate polynomial and transcendental functions. The cylindrical algebraic decomposition algorithm, used to solve systems of real polynomial equations and inequalities, applies validated numeric methods in the lifting phase. In my talk I will discuss the types of validation of numeric computation results used by the Mathematica solvers. I will present the validation criteria used and describe the algorithms using them.

Problem

Given an expression representing a holomorphic function f:ℂ⊇U->ℂ  find the solutions of f(x)==0 in a closed rectangle R⊂U.

A symbolic-numeric heuristic

SolveHeuristic(f, R, prec)
(1) Use a black-box numeric method with working precision prec to find approximate solutions.
(2) Apply symbolic criteria to the approximate solutions in an attempt to find exact isolating sets and multiplicities of the solutions and to prove that we have found all solutions.
(3) If (2) succeeds return a symbolic representation of the solutions otherwise return failed.

A symbolic-numeric algorithm

If certain additional conditions are satisfied (shown in blue text on the following slides) the following algorithm finds the solutions.

SolveAlgorithm(f, R)
(1) Set prec=machine double precision, answer=failed, 
(2) While (answer==failed) 
	(a) Set answer=SolveHeuristic(f, R, prec).
	(b) Set prec=2*prec.
(3) Return answer.

Example

Find the solutions of  3 (x^3)+(e^x)-sin(cos(2 x-1))+1==0  in the rectangle [-1,1]×[-1,1].

Picture

f=(E^x)-Sin[Cos[2x-1]]+3 (x^3)+1;
ff=f/. x->a+b I;
{ref,imf}=ComplexExpand[{Re[ff],Im[ff]}];
pl=ContourPlot[{ref,imf},{a,-1.5,1.5},{b,-1.5,1.5}];
rts=({Red,PointSize[Large],Point[{Re[#1],Im[#1]}]}&)/@(x/. Solve[f==0&&-1≤Re[x]≤1&&-1≤Im[x]≤1,x]);
Show[{pl,Graphics[{Green,Opacity[0.1],Rectangle[{-1,-1},{1,1}]}],Graphics[rts]}]

Solution

Solve[(E^x)-Sin[Cos[2x-1]]+3 (x^3)+1==0&&-1≤Re[x]≤1&&-1≤Im[x]≤1,x]

The solutions represent exact numbers.

x/.%[[1]]

N[%,100]

Assumptions

Given an expression representing a holomorphic function f:ℂ⊇U->ℂ  find the solutions of f(x)==0 in a closed rectangle R⊂U.

We assume that for any n∈ℕ

• an expression representing (f^((n))) can be computed,

• an interval-arithmetic evaluation is available for (f^((n))).

Interval-arithmetic evaluation

A complex interval is an element of the set 𝔻:={(c,r):c∈ℚ[𝕚], r∈(SubPlus[ℚ])⋃{∞}}.

A complex interval I==(c(I),r(I))∈𝔻 represents the set (SubStar[I])=={x∈ℂ:|x-c(I)|≤ r(I)}.

Let g be an expression representing a holomorphic function. An interval-arithmetic evaluation for g is an algorithm (Subscript[E, g]) which for any I∈𝔻 computes (Subscript[E, g])(I)∈𝔻 with the following properties:

• g((SubStar[I]))⊆(SubStar[(Subscript[E, g])(I)])

• (SubStar[I])⊆(SubStar[J])⟹(SubStar[(Subscript[E, g])(I)])⊆(SubStar[(Subscript[E, g])(J)])

• For any c∈D(g) and ϵ>0 there exists δ>0 such that 
  
   r<δ∧(d,s)==(Subscript[E, g])((c,r)) ⟹ s<ϵ

Approximate solutions

Given an expression representing a holomorphic function f:ℂ⊇U->ℂ  find the solutions of f(x)==0 in a closed rectangle R⊂U.

Let (Subscript[ξ, 1]),…,(Subscript[ξ, k]) be the distinct solutions of f(x)==0  in R.

We assume that there is an algorithm ApproximateSolutions such that

• ApproximateSolutions(f, R, prec) gives a list of numbers in R.

• For prec≥(Subscript[p, 0])  card(ApproximateSolutions(f,R,prec))==k

• For any ϵ>0 there exists (Subscript[p, ϵ])≥(Subscript[p, 0]) such that for any prec≥(Subscript[p, ϵ])
  ApproximateSolutions(f,R,prec)=={(Subscript[r, 1]),…,(Subscript[r, k])}  and (max_(1≤i≤k))|(Subscript[r, i])-(Subscript[ξ, i])|<ϵ.

Root existence criterion

For  I∈𝔻  let us denote (Subscript[I, min]):=(min_(z∈(SubStar[I])))|z| and (Subscript[I, max]):=(max_(z∈(SubStar[I])))|z|.

Criterion: Let f:ℂ⊇U->ℂ  be a holomorphic function, let I∈𝔻 be such that (SubStar[I])⊂U. Suppose that, for some k∈(SubPlus[ℕ]) and J==(c(I),r)∈𝔻,  (SubStar[J])⊂U

((Subscript[((Subscript[E, (f^((k)))])(J)), min])/(k!))(r^k)>(∑_(i=0)^(k-1))((Subscript[((Subscript[E, (f^((i)))])(I)), max])/(i!))(r^i)

then f has exactly k roots in (SubStar[J]) (counted with multiplicities).

Proof

Let (Subscript[z, 0]):=c(I), let φ:=Arg(c((Subscript[E, (f^((k)))])(J))) and let g(z):=exp(-𝕚 φ)f(z).

For z∈(SubStar[J]) we have

|Re((((g^((k)))(z))/(k!)))|≥((Subscript[((Subscript[E, (f^((k)))])(J)), min])/(k!))>(∑_(i=0)^(k-1))((Subscript[((Subscript[E, (f^((i)))])(I)), max])/(i!))(r^(i-k)) ≥(∑_(i=0)^(k-1))|(((f^((i)))((Subscript[z, 0])))/(i!))|(r^(i-k))=(∑_(i=0)^(k-1))|(((g^((i)))((Subscript[z, 0])))/(i!))|(r^(i-k))

The criterion follows from applying the following theorem to g and (SubStar[J]).

Theorem (Neumaier, 1988) Let f:ℂ⊇U->ℂ  be a holomorphic function and let (Subscript[z, 0])∈U and r>0 be such that D:= {z:|z-(Subscript[z, 0])|≤ r}⊂U. If for all z∈D
|Re((((f^((k)))(z))/(k!)))|>(∑_(i=0)^(k-1))|(((f^((i)))((Subscript[z, 0])))/(i!))|(r^(i-k)) 
then f has exactly k roots in D (counted with multiplicities).

Application of the root existence criterion

Criterion: Let f:ℂ⊇U->ℂ  be a holomorphic function, let I∈𝔻 be such that (SubStar[I])⊂U. Suppose that, for some k∈(SubPlus[ℕ]) and J==(c(I),r)∈𝔻,  (SubStar[J])⊂U

((Subscript[((Subscript[E, (f^((k)))])(J)), min])/(k!))(r^k)>(∑_(i=0)^(k-1))((Subscript[((Subscript[E, (f^((i)))])(I)), max])/(i!))(r^i)

then f has exactly k roots in (SubStar[J]) (counted with multiplicities).

We attempt to find an isolated solution of f==0 starting with an approximate solution (Subscript[z, 0])∈ApproximateSolutions(f,R,prec). If successful, we get r>0 and k∈(SubPlus[ℕ]) such that D((Subscript[z, 0]),r) contains exactly k solutions of f==0. If k>1 the method cannot distinguish between one solution of multiplicity k and a cluster of solutions. If f is an elementary function one could use zero-testing to prove multiplicity of the solution (the currently known zero-testing algorithm [Richardson] relies on Schanuel’s conjecture for proof of termination).

(1) Set k:=1.   If (Subscript[z, 0])≠0 set ϵ:=|(Subscript[z, 0])|(10^(-2prec)) else set  ϵ:=(10^(-2 prec)).   Put I:=((Subscript[z, 0]),ϵ) and (Subscript[M, 0]):=((Subscript[((Subscript[E, f])(I)), max])/(i!)).

(2) While k≤(Subscript[k, max]) do
	(a) Compute (Subscript[m, k]):=((Subscript[((Subscript[E, (f^((k)))])(I)), min])/(k!)) and (Subscript[M, k]):=((Subscript[((Subscript[E, (f^((k)))])(I)), max])/(k!)).
	(b) If (Subscript[m, k])==0, increment k and continue the loop.
	(c) Set r:= (((max_(0≤i≤k-1))(2k((Subscript[M, i])/(Subscript[m, k]))))^(1/(k-i))) .
	(d) If  ((Subscript[((Subscript[E, (f^((k)))])(((Subscript[z, 0]), r))), min])/(k!))(r^k)>(∑_(i=0)^(k-1))(Subscript[M, i])(r^i) return (r,k).
	(e) Increment k and continue the loop.

(3) Return failed.

(Subscript[k, max]) is a bound on the total number of solutions, if known, otherwise it is a fixed parameter. ϵ and r are rational numbers -- floating point number approximations of the values given in (1) and (2c). If (Subscript[z, 0]) is close to a root ξ of f of multiplicity μ, k≤μ and ((Subscript[m, k])/(Subscript[M, k]))≈1 then r≲2(μ^2)|(Subscript[z, 0])-ξ|. Hence if for some k, ((Subscript[m, k])/(Subscript[M, k]))≈1 and r is large we can return failed without continuing the loop all the way to (Subscript[k, max]).

Verifying the number of solutions

Given an expression representing a holomorphic function f:ℂ⊇U->ℂ  find the solutions of f(x)==0 in a closed rectangle R⊂U.

To find the number of solutions of f(x)==0 in R (counted with multiplicities) we compute the winding number of f along the boundary of R. We subdivide the boundary into segments such that for each segment either 0∉Re((Subscript[E, f])(I)) or 0∉Im((Subscript[E, f])(I)). This requires that there are no zeros of f on the boundary of R. We use rectangular interval arithmetic here.

Finding approximate solutions

ApproximateSolutions(f, R, prec) uses either computational geometry methods or the Kravanja-Van Barel algorithm to find machine double precision solutions and uses them as starting points to the Newton method. Satisfying the SolveAlgorithm conditions would require an arbitrary-precision version of the part finding the starting points (currently not implemented in Mathematica).

Semialgebraic sets

A real polynomial condition is a formula

f((Subscript[x, 1]),…,(Subscript[x, n]))ρ 0

where f∈ℝ[(Subscript[x, 1]),…,(Subscript[x, n])] and ρ is one of <,≤,>,≥,=, or ≠.

A quantified real polynomial formula is a formula constructed with real polynomial conditions using Boolean operators and quantifiers over real variables.

A subset of  (ℝ^n) is semialgebraic if it is a solution set of a real polynomial formula with n free variables.

Cell decomposition

Every semialgebraic set can be represented as a finite union of disjoint cells defined recursively as follows.

• A cell in ℝ is a point or an open interval.

• A cell in (ℝ^k) has one of the two forms

{((Subscript[a, 1]),…,(Subscript[a, k]),(Subscript[a, k+1])):((Subscript[a, 1]),…,(Subscript[a, k]))∈(Subscript[C, k])∧(Subscript[a, k+1])==r((Subscript[a, 1]),…,(Subscript[a, k]))}

{((Subscript[a, 1]),…,(Subscript[a, k]),(Subscript[a, k+1])):((Subscript[a, 1]),…,(Subscript[a, k]))∈(Subscript[C, k])∧(Subscript[r, 1])((Subscript[a, 1]),…,(Subscript[a, k]))<(Subscript[a, k+1])<(Subscript[r, 2])((Subscript[a, 1]),…,(Subscript[a, k]))}

where (Subscript[C, k]) is a cell in (ℝ^k), r  is a continuous algebraic function, and (Subscript[r, 1]) and (Subscript[r, 2]) are continuous algebraic functions, -∞, or ∞, and (Subscript[r, 1])<(Subscript[r, 2]) on (Subscript[C, k]).

By an algebraic function we mean a function

(Subscript[Root, (Subscript[x, k+1]),p])f:(Subscript[C, k])∋((Subscript[x, 1]),..., (Subscript[x, k]))⟶(Subscript[Root, (Subscript[x, k+1]),p])f((Subscript[x, 1]),..., (Subscript[x, k]))∈ℝ

where

f=(Subscript[c, 0]) (Subscript[x, k+1]^m)+(Subscript[c, 1]) (Subscript[x, k+1]^(m-1))+…+(Subscript[c, m])∈ℝ[(Subscript[x, 1]),…,(Subscript[x, k]),(Subscript[x, k+1])]

is a polynomial and

(Subscript[Root, (Subscript[x, k+1]),p])f((Subscript[x, 1]),..., (Subscript[x, k]))

is the p-th real root (multiplicities counted) of f treated as a univariate polynomial in (Subscript[x, k+1]).

Cylindrical algebraic decomposition

The CAD algorithm, introduced by Collins (1975), allows to compute a cell decomposition of any semialgebraic set presented by a real polynomial formula.

• Cells are arranged cylindrically.

• Polynomials whose roots bound cells extending C are delineable over C.

A set A⊂ℝ[(Subscript[x, 1]),…,(Subscript[x, k]),(Subscript[x, k+1])] is delineable over C⊆(ℝ^k) if each f∈A has a fixed number of real roots on C as a polynomial in (Subscript[x, k+1]), the roots are continuous functions on C, they have constant multiplicities, and two roots of  f,g∈A  are equal either everywhere or nowhere in C.

ineqs=-5 (((y-3))^3)+7 (x^2) (y-3)<1&&(x^2)+2 (((y-4))^2)<21;

pl=RegionPlot[ineqs,{x,-4.5,4.5},{y,-0.5,7.75},AspectRatio->Automatic,PlotPoints->100];
cad=CylindricalDecomposition[ineqs,{x,y}];
bds=Select[cad/.{Or->List,And->List,HoldPattern[Inequality[a_,__,b_]]:>{a,b}},ListQ[#[[1]]]&];
pts=Sort[Union@@(First/@bds),Less];
cbs=Table[{Graphics[{c=Hue[(i-1)/7],Thickness[0.005],Line[{{pts[[i]],0},{pts[[i+1]],0}}]}],Plot[bds[[i,2]],{x,bds[[i,1,1]],bds[[i,1,2]]},ColorFunction->(c&),PlotStyle->Thickness[0.005],PlotRange->{{-4.5,4.5},{-0.5,7.75}}]},{i,7}];
Show[{pl,cbs,Graphics[{Thickness[0.002],Gray,Line[{{#,-0.5},{#,7.75}}]&/@pts}]},Axes->True,Frame->False]

CAFView[CylindricalDecomposition[ineqs,{x,y}],{x,y}]

CAFView code

cview[cad_,vars_,rul_,dvars_,ind_]:=
Module[{k=Length[ind],cc,ee,a,b,ans=False,ha,hb,x,na,rr},
ha=Hue[k/Length[vars],1/3];
hb=Hue[(k+1)/Length[vars],1/3];
If[rul=!={},
a=cad/.rul;
If[a===True||a===False,Return[Row[{a},Background->ha,Frame->True]]]];
x=dvars[[k+1]];
cc=If[Head[cad]===Or,List@@cad,{cad}];
If[k==Length[vars]-1,
Return[Or@@(Row[{Tooltip[#,N[#/.rul]]},Background->ha,Frame->True]&/@cc)]];
ee=cc[[1]];
If[Head[ee]===And,ee=ee[[1]]];
If[!FreeQ[ee,Alternatives@@Drop[vars,k+1]],
Return[Row[{True},Background->ha,Frame->True] &&Row[{Opener[Dynamic[Evaluate[x]]],PaneSelector[{False->Row[{Subscript["b",Sequence@@ind,1]},Frame->True,Background->hb],True->Row[{cview[cad,vars,rul,dvars,Append[ind,1]]},Frame->True]},Dynamic[Evaluate[x]],ImageSize->Automatic]}]]];
Do[ee=cc[[i]];
{a,b}=If[Head[ee]===And, {ee[[1]],Drop[ee,1]},{ee,True}];
na=N[a/.rul];
rr=If[na===True,rul,{}];
If[na=!=False,ans=ans||Row[{Tooltip[a,na]},Background->ha,Frame->True] &&Row[{Opener[Dynamic[Evaluate[x]]],PaneSelector[{False->Row[{Subscript["b",Sequence@@ind,i]},Frame->True,Background->hb],True->Row[{cview[b,vars,rr,dvars,Append[ind,i]]},Frame->True]},Dynamic[Evaluate[x]],ImageSize->Automatic]}]],
{i,Length[cc]}];
ans]

CAFView[cad_,vars_,rul_:{}]:=With[{dvars=Table[Unique["dynamicx"],{Length[vars]}]},cview[cad,vars,rul,dvars,{}]]

The projection phase of the CAD algorithm

Finding a cell decomposition of a semialgebraic set using the CAD algorithm consists of two phases, projection and lifting. In the projection phase we start with the set (Subscript[A, n]) of  factors of the polynomials present in the system, and eliminate variables one by one using a projection operator P such that

(Subscript[P, k+1]):ℝ[(Subscript[x, 1]),…,(Subscript[x, k]),(Subscript[x, k+1])]⊃(Subscript[A, k+1])⟶(Subscript[A, k])⊂ℝ[(Subscript[x, 1]),…,(Subscript[x, k])]

and, generally speaking, if all polynomials of (Subscript[A, k]) have constant signs on a cell C⊂(ℝ^k), then all polynomials of (Subscript[A, k+1]) are delineable over C. This way the roots of polynomials of (Subscript[A, 1]),…,(Subscript[A, n]) are the algebraic functions needed in the construction of the cell decomposition of the semialgebraic set.

The lifting phase of the CAD algorithm

In the lifting phase we find a cell decomposition of the semialgebraic set. We start with cells in (ℝ^1) consisting of all distinct roots of (Subscript[A, 1]) and the open intervals between the roots. We find a sample point in each of the cells and remove the cells whose sample points do not satisfy the system describing the semialgebraic set (the system may contain conditions involving only (Subscript[x, 1])). Now we lift the cells to cells in (ℝ^n), one dimension at a time. Suppose we have lifted the cells to (ℝ^k).  To lift a cell C⊂(ℝ^k) to (ℝ^(k+1)) we find the real roots of (Subscript[A, k+1])(c), where c is a sample point in C and the elements of (Subscript[A, k+1])(c) are the elements of  (Subscript[A, k+1]) with (Subscript[x, 1]),(Subscript[…x, k]) replaced with the coordinates of c.  If r is the p-th root of f(c,(Subscript[x, k+1])) for some f∈(Subscript[A, k+1]), then (Subscript[Root, (Subscript[x, k+1]),p])f((Subscript[x, 1]),..., (Subscript[x, k])) is a continuous algebraic function on C, because the polynomials of (Subscript[A, k+1]) are delineable on C. The lifting of C to (ℝ^(k+1)) will consist of graphs of such algebraic functions, and of the slices of C×ℝ between the consecutive graphs. The sample points in each of the new cells will be obtained by adding the k+1-st coordinate to c, equal to one of the roots or to a number between two consecutive roots. Similarly as in the first step we remove those lifted cells whose sample points do not satisfy the system describing the semialgebraic set.

To lift a cell we need to find the real root structure of (Subscript[A, k+1])(c). That is, we need to find disjoint intervals isolating the distinct real roots of (Subscript[A, k+1])(c), find the multiplicity of each root and identify the common roots of different elements of (Subscript[A, k+1])(c). The coefficients of sample points are in general algebraic numbers. We can avoid algebraic number computations by using a validated numeric method.

Problem (Real root structure)

Given

• a finite set of polynomials A⊂ℤ[(Subscript[x, 1]),…,(Subscript[x, k]),(Subscript[x, k+1])], where k≥0,

• a tuple ((Subscript[I, 1]),…,(Subscript[I, k])) of real intervals such that (Subscript[a, i])∈(Subscript[I, i]) for 1≤i≤k,

find

• a tuple ((Subscript[J, 1]),…,(Subscript[J, l])) of disjoint real intervals,

• a function mult:A×{1, . . . ,l}⟶ℕ

such that

• the polynomials {f(a,(Subscript[x, k+1])):f∈A} have l distinct real roots (Subscript[r, 1]),…,(Subscript[r, l]),

• (Subscript[r, i])∈(Subscript[J, i]) for 1≤i≤l,

• mult(f,i) is the multiplicity of (Subscript[r, i]) as a root of f for all f∈A and 1≤i≤l.

Additional assumptions

For any f,g∈A we can obtain the following information

• the degree of the polynomial f(a,(Subscript[x, k+1])),

• the degree of the g.c.d. of f(a,(Subscript[x, k+1])) and g(a,(Subscript[x, k+1])),

• the degree of the g.c.d. of f(a,(Subscript[x, k+1])) and (∂_(Subscript[x, k+1])) f(a,(Subscript[x, k+1])).

Obtaining the additional information in the CAD algorithm

This information may be deduced from knowing the number of initial coefficients of  f(a,(Subscript[x, k+1])) that are zero and the number of initial principal subresultant coefficients (PSC) of f(a,(Subscript[x, k+1])) and g(a,(Subscript[x, k+1])) or of  f(a,(Subscript[x, k+1])) and (∂_(Subscript[x, k+1])) f(a,(Subscript[x, k+1])) that are zero.

In the CAD algorithm, the projection (Subscript[A, k]) of the polynomial set (Subscript[A, k+1]) contains the nonconstant initial coefficients, the discriminants and the pairwise resultants for all polynomials in (Subscript[A, k+1]). The projection may also contain PSC for pairs of polynomials or for pairs polynomial and its derivative.

The signs of projection polynomials (Subscript[A, k]) on a given cell cell C⊂(ℝ^k) are constant. If a∈C then for a polynomial g∈(Subscript[A, k]), g(a)==0 iff C is the graph of a root of g. Hence if a coefficient or a PSC belongs to the projection, the information on whether it is zero at a comes for free from the construction of  C.  We just need to keep track of where the projection polynomials come from and which projection polynomials are zero on the given cell.

If we need to know whether a PSC that is not a part of the projection is zero at a, we can use zero testing to obtain this information. This may require nontrivial computation, but still it tends to be faster than lifting the cell using exact algebraic number computation.

Timings

Symbolic-numeric real root structure computation

RealRootStructureHeuristic(A, I, prec)
(1) Let F:={g(c((Subscript[I, 1])),…,c((Subscript[I, k])),z):g∈A}⊂ℚ[z] 
(2) Use a black-box numeric method with working precision prec to find approximations of all complex roots of elements of F.
(3) Apply symbolic criteria to the approximate solutions in an attempt to find the real root structure of A at a.

A symbolic-numeric algorithm for cell lifting in CAD

CellLifting(projection data)
(1) Construct cells (Subscript[C, 1]),…,(Subscript[C, m]) in ℝ (exact real root isolation).
(2) Set prec=initial precision, answer=empty.
(3) For 1≤i≤m,
	(a) Let a be a sample point in (Subscript[C, i]). If a∈ℚ, set I:=(a,0), else find I==(c,r) such that c-r<a<c+r and r<|c|(10^(-prec)).
	(b) Append CellLiftingRecursive((Subscript[C, i]),I,prec,projection data) to answer.
(4) Return answer.

CellLiftingRecursive(C,I,prec,projection data)
(1) Set (Subscript[I, C]):=I, newprec=prec.
(2) Compute  root structure=RealRootStructureHeuristic((Subscript[A, k+1]),(Subscript[I, C]),newprec).
(3) While (root structure==failed) 
	(a) Set newprec = 2 newprec.
	(b) Recompute (Subscript[I, C]) working with precision newprec.
	(c) Compute root structure=RealRootStructureHeuristic((Subscript[A, k+1]),(Subscript[I, C]),newprec).
(4) Construct cells (Subscript[C, 1]),…,(Subscript[C, m]) in (ℝ^(k+1)) and intervals (Subscript[I, 1]),…,(Subscript[I, m]) in (𝕀^(k+1))  extending C and (Subscript[I, c]). 
(5) Set answer=empty. For 1≤i≤m, append CellLiftingRecursive((Subscript[C, i]),(Subscript[I, i]),newprec,projection data) to answer.
(6) Return answer.

An alternative version of the algorithm, currently used by the Mathematica CAD implementation, is to revert to exact algebraic number computations for the cell C if the call to RealRootStructureHeuristic in step (2) of CellLiftingRecursive fails. The termination of the version of the algorithm presented here depends on the following property of the numeric root finding method NumericRoots(f, prec).

For any polynomial f∈ℚ[z] and any ϵ>0 there exists (Subscript[p, ϵ]) such that for any prec≥(Subscript[p, ϵ]) 
NumericRoots(f,prec)=={(Subscript[r, 1]),…,(Subscript[r, n])} and (max_(1≤i≤n))|(Subscript[r, i])-(Subscript[ξ, i])|<ϵ
where f(z)==(Subscript[a, n])(z-(Subscript[ξ, 1]))…(z-(Subscript[ξ, n]))

Interval-arithmetic evaluation of polynomials

Let us slightly modify the definition of complex interval.

A complex interval is an element of the set 𝔻:={(c,r):c∈ℚ[𝕚], r∈(SubPlus[ℚ])⋃{0,∞}}.

A real interval is an element of the set 𝕀:={(c,r):c∈ℚ, r∈(SubPlus[ℚ])⋃{0,∞}}. We have 𝕀⊂𝔻.

We allow one-point intervals here -- unlike for analytic functions, we can always compute the exact value of a rational polynomial at a Gaussian rational point.

For I==((Subscript[I, 1]),…,(Subscript[I, k]))∈(𝔻^k) let (SubStar[I]):=(SubStar[((Subscript[I, 1]))])×…×(SubStar[((Subscript[I, k]))])⊆(ℂ^k).

Let g∈ℚ[(Subscript[x, 1]),…,(Subscript[x, k])]. An interval-arithmetic evaluation for g is an algorithm (Subscript[E, g]) which for any I∈(𝔻^k) computes (Subscript[E, g])(I)∈𝔻 with the following properties:

• (SubStar[g])((SubStar[I]))⊆(SubStar[(Subscript[E, g])(I)])

• (SubStar[I])⊆(SubStar[J])⟹(SubStar[(Subscript[E, g])(I)])⊆(SubStar[(Subscript[E, g])(J)])

• For any a∈(ℂ^k) and ϵ>0 there exists δ>0 such that 

  I==(((Subscript[a, 1]),(Subscript[r, 1])),…,((Subscript[a, k]),(Subscript[r, k])))∧(max_(1≤i≤k))(Subscript[r, i])<δ∧(c,r)==(Subscript[E, g])(I) ⟹r<ϵ

• If I∈(𝕀^k) then (Subscript[E, g])(I)∈𝕀.

For f==(Subscript[I, n])(z^n)+…+(Subscript[I, 0])∈𝔻[z] and (Subscript[z, 0])∈ℚ[𝕚] define

f((Subscript[z, 0])):=(Subscript[E, g])((Subscript[I, 0]),…,(Subscript[I, n]),((Subscript[z, 0]),0))∈𝔻

where g==(Subscript[x, n])(z^n)+…+(Subscript[x, 0])∈ℚ[(Subscript[x, 0]),…,(Subscript[x, n]),z].

Root existence criterion

For  I∈𝔻  let us denote (Subscript[I, min]):=(min_(z∈(SubStar[I])))|z| and (Subscript[I, max]):=(max_(z∈(SubStar[I])))|z|.

Criterion: Let f==(Subscript[I, n])(z^n)+…+(Subscript[I, 0])∈𝔻[z] with 0∉(Subscript[I, n]), let (Subscript[z, 0])∈ℚ[𝕚] and let (Subscript[c, i]):=(((f^((i)))((Subscript[z, 0])))/(i!))∈𝔻 for 0≤i≤n. Suppose that

(((max_(0≤i<m))(((n (Subscript[((Subscript[c, i])), max]))/(Subscript[((Subscript[c, m])), min]))))^(1/(m-i)))<r<(((min_(m<i≤n))(((Subscript[((Subscript[c, m])), min])/(Subscript[n((Subscript[c, i])), max]))))^(1/(i-m)))

then for any  (Subscript[a, 0])∈(Subscript[I, 0]),…,(Subscript[a, n])∈(Subscript[I, n])  the polynomial (Subscript[a, n])(z^n)+…+(Subscript[a, 0])∈ℂ[z] has exactly m roots in the disk D((Subscript[z, 0]),r) (counted with multiplicities).

Proof

The criterion follows from:

Proposition  (A.S., JSC 2006, proof based on Rouche’s theorem)
Let f∈ℂ[x] be a polynomial of degree n, let (Subscript[z, 0])∈ℂ and let (Subscript[c, i]):=|(((f^((i)))((Subscript[z, 0])))/(i!))|. Suppose that

(((max_(0≤i<k))(((n (Subscript[c, i]))/(Subscript[c, k]))))^(1/(k-i)))<r<(((min_(k<i≤n))(((Subscript[c, k])/(n (Subscript[c, i])))))^(1/(i-k)))

then f has exactly k roots in the disk D((Subscript[z, 0]),r) (counted with multiplicities).

Application of the root existence criterion

Criterion: Let f==(Subscript[I, n])(z^n)+…+(Subscript[I, 0])∈𝔻[z] with 0∉(Subscript[I, n]), let (Subscript[z, 0])∈ℚ[𝕚] and let (Subscript[c, i]):=(((f^((i)))((Subscript[z, 0])))/(i!))∈𝔻 for 0≤i≤n. Suppose that

(((max_(0≤i<m))(((n (Subscript[((Subscript[c, i])), max]))/(Subscript[((Subscript[c, m])), min]))))^(1/(m-i)))<r<(((min_(m<i≤n))(((Subscript[((Subscript[c, m])), min])/(Subscript[n((Subscript[c, i])), max]))))^(1/(i-m)))

then for any  (Subscript[a, 0])∈(Subscript[I, 0]),…,(Subscript[a, n])∈(Subscript[I, n])  the polynomial (Subscript[a, n])(z^n)+…+(Subscript[a, 0])∈ℂ[z] has exactly m roots in the disk D((Subscript[z, 0]),r) (counted with multiplicities).

Problem

Given

• a finite set of polynomials A⊂ℤ[(Subscript[x, 1]),…,(Subscript[x, k]),(Subscript[x, k+1])], where k≥0,

• a tuple ((Subscript[I, 1]),…,(Subscript[I, k])) of real intervals such that (Subscript[a, i])∈(Subscript[I, i]) for 1≤i≤k,

find

• a tuple ((Subscript[J, 1]),…,(Subscript[J, l])) of disjoint real intervals,

• a function mult:A×{1, . . . ,l}⟶ℕ

such that

• the polynomials {f(a,(Subscript[x, k+1])):f∈A} have l distinct real roots (Subscript[r, 1]),…,(Subscript[r, l]),

• (Subscript[r, i])∈(Subscript[J, i]) for 1≤i≤l,

• mult(f,i) is the multiplicity of (Subscript[r, i]) as a root of f for all f∈A and 1≤i≤l.

Find the root structure of each polynomial

For f∈A do

Find root approximations

(1) Compute (Subscript[f, I])==(Subscript[I, n])(z^n)+…+(Subscript[I, 0])==f((Subscript[I, 1]),…,(Subscript[I, k]),z)∈𝔻[z].
(2) While n≥0 and 0∈(Subscript[I, n])
	(a) If the coefficient at (Subscript[x, k+1]^n) of f(a,(Subscript[x, k+1])) is nonzero, return failed.
	(b) Set (Subscript[f, I])=(Subscript[I, n-1])(z^(n-1))+…+(Subscript[I, 0]) and decrement n.
(3) If (Subscript[f, I])==0 set (Subscript[R, f]^real)=(Subscript[R, f]^complex):=⌀ and continue the loop.
(4) Let (Subscript[f, c]) be f(c((Subscript[I, 1])),…,c((Subscript[I, k])),z) with monomials of degree higher than n removed. 
(5) Compute approximate roots (Subscript[z, 1]),…,(Subscript[z, n])∈ℚ[𝕚] of (Subscript[f, c]).

Use the criterion to find disjoint disks containing roots of f(a,(Subscript[x, k+1]))

(6) For 1≤i≤n find the smallest positive integer (Subscript[m, i]) and a radius (Subscript[r, i]) such that the criterion is satisfied for (Subscript[f, I]), (Subscript[z, i]), (Subscript[m, i]) and (Subscript[r, i]).
(7) Let R:={((Subscript[z, i]),(Subscript[r, i]),(Subscript[m, i])):1≤i≤n}.

Show[Graphics[{Opacity[0.5],Red, Disk[{1, 2},0.2], Disk[{1.05, 2.05},0.2],Disk[{0.9, 1.95},0.25], Disk[{-1.5, 1},0.15], Disk[{2,-0.1},0.2],Disk[{1, -2},0.2], Disk[{1.05, -2.05},0.2],Disk[{0.9, -1.95},0.25], Disk[{-1.5, -1},0.15]}],Axes->True,PlotRange->{{-3,3},{-3,3}}]

(8) For each (z,r,m)∈R such that z∉ℝ∧D(z,r)⋂ℝ≠⌀
	(a) Find the smallest positive integer m' and a radius r' such that the criterion is satisfied for (Subscript[f, I]), z', m' and re(r).
	(b) Replace (z,r,m) with (re(z),r',m') in R. 
(8) While there exist (z,r,m),(z',r',m')∈R, D(z,r)⋂D(z',r')≠⌀ and r≤r', remove (z',r',m') from R.

(9) Set (Subscript[R, f]^real):={(z,r,m)∈R:z∈ℝ} and (Subscript[R, f]^complex):={(z,r,m)∈R:im(z)>0}.

Verify that the disks contain all roots, each element of (Subscript[R, f]^real) contains one root and the root is real

(10) If  (∑_((z,r,m)∈(Subscript[R, f]^real)))m+2(∑_((z,r,m)∈(Subscript[R, f]^complex)))m≠n return failed.
(11) If there exists (z,r,m)∈(Subscript[R, f]^real) such that m>1 find the degree d of the g.c.d. of f(a,(Subscript[x, k+1])) and (∂_(Subscript[x, k+1])) f(a,(Subscript[x, k+1])). If card((Subscript[R, f]^real))+2card((Subscript[R, f]^complex))≠n-d return failed.

Identify the common real roots

For f,g∈A such that there exist (z,r,m)∈(Subscript[R, f]^real) and  (z',r',m')∈(Subscript[R, g]^real)  with D(z,r)⋂D(z',r')≠⌀ do

(1) Set (Subscript[m, tot])=0. 
(2) For each  (z,r,m)∈(Subscript[R, f]^real) and  (z',r',m')∈(Subscript[R, g]^real)  such that D(z,r)⋂D(z',r')≠⌀  set (Subscript[m, tot])=(Subscript[m, tot])+min(m,m').
(3) For each  (z,r,m)∈(Subscript[R, f]^complex) and  (z',r',m')∈(Subscript[R, g]^complex)  such that D(z,r)⋂D(z',r')≠⌀  set (Subscript[m, tot])=(Subscript[m, tot])+2min(m,m').
(4) Find the degree d of the g.c.d. of f(a,(Subscript[x, k+1])) and g(a,(Subscript[x, k+1])).
(5) If (Subscript[m, tot])≠d return failed.

Return isolating intervals for the real roots and the multiplicity function

(1) Find real intervals ((Subscript[J, 1]),…,(Subscript[J, l])) by picking one representative from each set of intersecting intervals in {(z,r): f∈A ∧(z,r,m)∈(Subscript[R, f]^real)}.
(2) For f∈A and  1≤i≤l, if there exists (z,r,m)∈(Subscript[R, f]^real) such that (Subscript[J, i])⋂D(z,r)≠⌀ then mult(f,i)=m else mult(f,i)=0.
(3) Return ((Subscript[J, 1]),…,(Subscript[J, l])) and mult.