Re: Root Finding Methods Gaurenteed to Find All Root

• To: mathgroup at smc.vnet.net
• Subject: [mg112788] Re: Root Finding Methods Gaurenteed to Find All Root
• From: Daniel Lichtblau <danl at wolfram.com>
• Date: Thu, 30 Sep 2010 04:53:11 -0400 (EDT)
• References: <001801cb5fd0\$32430380\$96c90a80\$@md.metrocast.net>

Ted Ersek wrote:
> It was mentioned that some of the methods to find all roots in a certain
> region require interval analysis.
> However, Interval arithmetic in Mathematica is limited. First of all
> Mathematica only knows how to evaluate
> elementary function and simple functions such as (Abs, Min) when given an
> Interval as an argument.
> So for example Mathematica does nothing with   Gamma[ Interval[ {2,3} ] ]
> even though Gamma[x] is continuous and monotonic over this interval with
> Gamma[2] == 1; Gamma[3] ==2
>
> Besides that Mathematica performs calculations as if each instance of
> Interval[{xmin_, xmax_}]
> is a value between (xmin, xmax) that is unrelated to any other instance of
> Interval[{xmin_, xmax_}].
> So I present the following as a worst case example of what this can give.
>
> In[1]:= poly=Function[x, Evaluate[ChebyshevT[15,x]]]
> Out[1]= Function[x,-15 x+560 x^3-6048 x^5+28800 x^7-70400 x^9+92160
> x^11-61440 x^13+16384 x^15]
>
> In[2]:= poly[Interval[{-1,1}]]
> Out[2]= Interval[{-275807,275807}]
>
> The above interval isn't very useful when we can see from  Plot[ poly[x],
> {x,-1,1} ]
> that  poly[x] oscillates between (-1,1) over this interval.
> For a general function we can approximate the interval we really want using
>
> NMinimize[ { f[x], xmin<x<xmax}, x ],  NMaximize[ { f[x], xmin<x<xmax}, x ].
>
> However, that's about as difficult as finding a root of f[x].
> Actually that may be more difficult than finding a root of f[x].
>
> Ted Ersek

For analytic univariate functions one can also use quadrature, subject
to numerical issues. Below is code for this purpose. It is a less buggy
version of code from

Caveat 1: It probably still has bugs.

Caveat 2: It uses a numeric epsilon. Intervals said to contain a unique
root might be off by epsilon (that is, the root might be outside). With
even more hackery than I already do, one could correct for this deficiency.

rectangleIntegrate[func_, a_, b_, c_, d_, z_, f2_:1] := Module[
{t1=0, t2=0, t3=0, t4=0, integrand=f2*D[func,z]/func, res},
t1 = NIntegrate[Evaluate[integrand], {z,a+I*c,b+I*c}];
t2 = NIntegrate[Evaluate[integrand], {z,b+I*c,b+I*d}];
t3 = NIntegrate[Evaluate[integrand], {z,b+I*d,a+I*d}];
t4 = NIntegrate[Evaluate[integrand], {z,a+I*d,a+I*c}];
res = (t1+t2+t3+t4)/(2*Pi*I);
res]

isolateRootsInRectangle[func_, a_, b_, c_, d_, z_, eps_] := Module[
{numroots, avgx=(a+b)/2, avgy=(c+d)/2,
r1, r2, r3, r4, rinter, len=0, allrts},
numroots =
Round[rectangleIntegrate[func,a-eps,b+eps,c-eps,d+eps,z]];
If[numroots==0, Return[{}]];
If[numroots==1, Return[{{a,b,c,d}}]];
If[c==d,
r1 = isolateRootsInRectangle[func,a,avgx,c,d,z,eps];
r2 = isolateRootsInRectangle[func,avgx,b,c,d,z,eps];
If [Length[r1]>=1 && Length[r2]>=1,
rinter = isolateRootsInRectangle[func,avgx,avgx,c,d,z,eps];
If [Length[rinter]==1,
r1 = isolateRootsInRectangle[func,a,avgx-2*eps,c,d,z,eps];
r2 = isolateRootsInRectangle[func,avgx+2*eps,b,c,d,z,eps]];
r2 = Join[r2,rinter];
];
,
r1 = isolateRootsInRectangle[func,a,avgx,c,avgy,z,eps];
r2 = isolateRootsInRectangle[func,a,avgx,avgy,d,z,eps];
If [Length[r1]>=1 && Length[r2]>=1,
rinter = isolateRootsInRectangle[func,a,avgx,avgy,avgy,z,eps];
If [Length[rinter]==1,
r1 = isolateRootsInRectangle[a,avgx,c,avgy-2*eps,z,eps];
r2 = isolateRootsInRectangle[func,a,avgx,avgy+2*eps,d,z,eps];
r2 = Join[r2,rinter];
];
];
];
If[a==b,
r3 = r4 = {};
,
If[c==d,
r3 = r4 = {};
,
r3 = isolateRootsInRectangle[func,avgx,b,c,avgy,z,eps];
r4 = isolateRootsInRectangle[func,avgx,b,avgy,d,z,eps];
If [Length[r3]>=1 && Length[r4]>=1,
rinter = isolateRootsInRectangle[func,avgx,b,avgy,avgy,z,eps];
If [Length[rinter]==1,
r3 = isolateRootsInRectangle[func,avgx,b,c,avgy-2*eps,z,eps];
r4 = isolateRootsInRectangle[func,avgx,b,avgy+2*eps,d,z,eps];
r4 = Join[r4,rinter];
];
];
];
];
If [a!=b && c!=d && Length[r1]>=1 && Length[r3]>=1,
rinter = isolateRootsInRectangle[func,avgx,avgx,avgy,avgy,z,eps];
If [Length[rinter]==1,
r1 = isolateRootsInRectangle[func,a,avgx-2*eps,c,avgy-2*eps,z,eps];
r3 = isolateRootsInRectangle[func,avgx+2*eps,b,avgy+2*eps,d,z,eps];
r3 = Join[r3, rinter];
];
];
Flatten[Join[{r1,r2,r3,r4}/.{}:>Sequence[]],1]
]

Here is the Chebyshev example.

poly = Function[x, Evaluate[ChebyshevT[15,x]]];

In[8]:= InputForm[rts = isolateRootsInRectangle[poly[z], -1, 1, 0, 0, z,
1/10^3]]

NIntegrate::slwcon:
Numerical integration converging too slowly; suspect one of the
following:
singularity, value of the integration is 0, highly oscillatory
integrand,
or WorkingPrecision too small.

NIntegrate::ncvb:
NIntegrate failed to converge to prescribed accuracy after 9
recursive bisections in z near {z} = {-0.864176 - 0.001 I}
-12
. NIntegrate obtained -5.37348 10    + 46.7339 I and 0.345719
for the integral and error estimates.

NIntegrate::slwcon:
Numerical integration converging too slowly; suspect one of the
following:
singularity, value of the integration is 0, highly oscillatory
integrand,
or WorkingPrecision too small.

NIntegrate::ncvb:
NIntegrate failed to converge to prescribed accuracy after 9
recursive bisections in z near {z} = {0.864176 + 0.001 I}
-12
. NIntegrate obtained -5.37348 10    + 46.7339 I and 0.345719
for the integral and error estimates.

NIntegrate::slwcon:
Numerical integration converging too slowly; suspect one of the
following:
singularity, value of the integration is 0, highly oscillatory
integrand,
or WorkingPrecision too small.

General::stop: Further output of NIntegrate::slwcon
will be suppressed during this calculation.

NIntegrate::ncvb:
NIntegrate failed to converge to prescribed accuracy after 9
recursive bisections in z near {z} = {-0.406078 - 0.001 I}
. NIntegrate obtained -4.07502 + 24.1523 I and 0.00208387
for the integral and error estimates.

General::stop: Further output of NIntegrate::ncvb
will be suppressed during this calculation.

Out[8]//InputForm=
{{-1, -15501/16000, 0, 0}, {-15501/16000, -7501/8000, 0, 0},
{-3501/4000, -1501/2000, 0, 0}, {-1501/2000, -2503/4000, 0, 0},
{-2503/4000, -501/1000, 0, 0}, {-501/1000, -503/2000, 0, 0},
{-503/2000, -1/500, 0, 0}, {1/500, 503/2000, 0, 0},
{503/2000, 501/1000, 0, 0}, {501/1000, 2503/4000, 0, 0},
{2503/4000, 1501/2000, 0, 0}, {1501/2000, 3501/4000, 0, 0},
{7501/8000, 15501/16000, 0, 0}, {15501/16000, 1, 0, 0}, {0, 0, 0, 0}}

Obviously the messages could be supressed. Also one might tweak
NIntegrate options e.g. via

SetOptions[NIntegrate,PrecisionGoal->6, AccuracyGoal->6,
WorkingPrecision->20];

in order to improve reliability.

Daniel Lichtblau
Wolfram Research

• Prev by Date: Re: Question on Solve
• Next by Date: Re: Interpolate in polar coordinates or cartesian
• Previous by thread: How to unflatten an array ?