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
http://groups.google.com/group/sci.math.symbolic/msg/c88809cbbc7fcb32
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