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