Services & Resources / Wolfram Forums / MathGroup Archive
-----

MathGroup Archive 2010

[Date Index] [Thread Index] [Author Index]

Search the Archive

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


  • 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 ?