Re: roots

*To*: mathgroup at smc.vnet.net*Subject*: [mg29595] Re: [mg29380] roots*From*: Daniel Lichtblau <danl at wolfram.com>*Date*: Wed, 27 Jun 2001 05:12:29 -0400 (EDT)*References*: <200106160647.CAA07743@smc.vnet.net>*Sender*: owner-wri-mathgroup at wolfram.com

maarten.vanderburgt at icos.be wrote: > > hallo, > > What is the easiest solution to finding the 13 (real) roots of 0.05 x + > Cos[x] == 0. > I can be easily seen that there are 13 roots by plotting: > Plot[{0.05 x,-Cos[x]},{x,-30,30}] > but this is not sufficient for finding accurate numerical values. > > FindRoot only gives one value and it is hard to predict which root it will find with a specified start value: > > In[83]:= FindRoot[0.05 x +Cos[x] == 0,{x,0}] > > Out[83]= {x -> -4.96317} > > There are three roots closer to 0 then x == -4.96317 > > An NSolve only works for polynomials. > > Is there no simple way to find all roots of such an equation, eventually within a specified range. > > thanks for your help > > Maarten van der Burgt > Leuven Simplicity is in the eye of the beholder. If you have a plot you trust, perhaps best is to use FindRoot on what appear to the eye to be close approximations, and then verify visually that you got all the roots. Below I show a more general approach that does not rely on Plot to find all the axis crossings. This method makes good use of the Cauchy integral formula. We will look for all roots between -20 and 20. To do so we will use contour integration in a narrow rectangle around this segment, using a vertical offset of 1/10. bound = 20; epsilon = 10^(-1); f[z_] := 1/20*z+Cos[z] Off[NIntegrate::slwcon] segmentIntegrate[func_,a_,b_,offset_,z_, f2_:1] := Module[{t1,t2,t3,integrand=f2*D[func,z]/func}, t1 = NIntegrate[ Evaluate[(integrand/.z->z-I*offset)- (integrand/.z->z+I*offset)], {z,a,b},MaxRecursion->10]; t2= NIntegrate[Evaluate[integrand], {z,a+I*offset,a-I*offset}, MaxRecursion->10]; t3 = NIntegrate[Evaluate[integrand], {z,b-I*offset,b+I*offset}, MaxRecursion->10]; (t1+t2+t3)/(2*Pi*I)] I grouped two terms together below rather than simply doing one path integral. This is to achieve cancellation on terms that might be problematic along long (say, infinite or semi-infinite) segments. For more general paths (see remark (vii) below) one would not do this. The above can be used to count roots on a segment. We use it for interval-bisection in order to isolate roots. isolateRoots[func_,lower_,upper_,z_,eps_] := Module[{numroots,avg=(lower+upper)/2,leftroots}, numroots = Round[segmentIntegrate[func,lower,upper,eps,z]]; If[numroots==0, Return[{}]]; If[numroots==1, Return[{Interval[{lower,upper}]}]]; leftroots = isolateRoots[func,lower,avg,z,eps]; If[Length[leftroots]==numroots,Return[leftroots]]; Join[leftroots,isolateRoots[func,avg,upper,z,eps]] ] We'll get intervals for the example at hand. In[11]:= InputForm[Timing[rootintervals = isolateRoots[f[z],-bound,bound,z,epsilon]]] Out[11]//InputForm= {2.1999999999999997*Second, {Interval[{-20, -75/4}], Interval[{-75/4, -35/2}], Interval[{-15, -25/2}], Interval[{-25/2, -10}], Interval[{-10, -5}], Interval[{-5, -5/2}], Interval[{-5/2, 0}], Interval[{0, 5/2}], Interval[{5/2, 5}], Interval[{5, 10}], Interval[{10, 25/2}], Interval[{25/2, 15}], Interval[{15, 20}]}} Now that we have one root in each interval, we use CIF again but with multiplier function z in order to evaluate the sum of all z-values for z a root. This simply gives the roots themselves. In[12]:= refineRoots[func_,intervals_,eps_,z_]:= Map[Re[segmentIntegrate[func,#[[1,1]],#[[1,2]],eps,z,z]]&, intervals] In[13]:= Timing[rts=refineRoots[f[z],rootintervals,epsilon,z]] Out[13]= {0.65 Second, {-19.1433, -18.4538, -13.4028, -11.6152, -7.47114, -4.96317, -1.49593, 1.65357, 4.48616, 8.28087, 10.446, 14.984, 16.324}} We check that the residuals are suitably small. In[14]:= Max[Abs[f[z] /. MapThread[{z -> #} &, {rts}]]] // InputForm Out[14]//InputForm= 4.854923130181987*^-10 Let me point out some of the subtle issues involved in this method. (i) It depends on numerical integration. I did no error checking but one should check that the root-counting phase gets values that are "close" to integers. (ii) Depending on the integrand, one may want to tune NIntegrate beyond setting the MaxRecursion option to a nondefault value. (iii) The function f[z] must be analytic in a neighborhood of the segment of interest. It must have no zeroes off the segment but in that neighborhood. This affects the choice of epsilon. On the other hand, making it too small may make NIntegrate work quite hard (possibly nondefault WorkingPrecision would be needed). (iv) In particular, branch cuts or poles will cause trouble. Poles, for example, will silently cancel with zeroes and thus mess up the root count. (v) The code I used above assumes that there are no multiple roots. A way to handle thos would be to isolate to some tolerance (say, 10^6*$MachineEpsilon). If the root counting (isolation) step still claims there is more than one root, call it a multiple root, do the refinement computation as above, and divide that value by the multiplicity computed in the isolation step. (vi) The code above assumes there is no root at an endpoint of the segment. It would be good to check endpoints explicitly before each isolation step. Any time a root is found at an endpoint, isolate it (accounting for multiplicity, if any, as per (v) above). Then move inward a small amount (that is, take a neighborhood of that endpoint), check that there are no extra roots in this neighborhood of the endpoint not accounted for by that endpoint, and take for a new endpoint a value midway between the original endpoint and the boundary of the small neighborhood. In this way we overlap with a small region that we know has no roots, but still miss the endpoint that was a root. (vii) These subtleties notwithstanding, the method is quite nice. For one, it generalizes to arbitrary regions provided one has a way to subdivide and to specify the boundaries for the contour integrations. (viii) Moreover it may be extended to higher dimensions (n equations in n variables). I am not conversant with the details, but will note that the theoretical groundwork was done by L. Aizenberg and colleages regarding multivariate residues. Practical issues in implementation (and some theory) were addressed in recent work by P. Kravanja. I believe he considers primarily polydisk regions but the same ideas will work for e.g. real regions in C^n. The main drawback so far as I am aware is that it is a challenge to do sufficiently fast and accurate quadrature as dimension goes up. So a multivariate extension may not be practical at this time beyond low dimensions. All the same, it provides an interesting approach to e.g. the "real solutions" problem wherein one has a polynomial system with many solutions but is only interested in the relatively small number of real ones. In this case complex solvers are generally not of much use, and local root-finders only help when reasonable starting points are known to all desired actual solutions. I would like to make some general remarks of an editorial-page nature. While it may not be a universally held opinion, one sometimes hears in this group that Mathematica is not terribly well suited for computations that are primarily numeric in nature. Indeed, a few months back the question arose as to what might be an ideal language in which to write a "foolproof" (so far as one can obtain) code to find all roots of a univariate function in a given interval; exactly the problem again at hand today. My opinion, which I back with the code above, is that Mathematica is an excellent candidate for this task. Moreover I note that the code is essentially all numeric. Yes, you can instead use C, load QUADPACK for the numeric integrations, and so forth (provided you don't need bignum capabilities). But I think the simplicity of the Mathematica approach, relative ease of development/extension/maintainance (and in particular the easy ability to influence the quadrature steps with options) argue well that Mathematica is a fine choice for the job. Daniel Lichtblau Wolfram Research

**References**:**roots***From:*maarten.vanderburgt@icos.be

**Re: Unappropiate context in a package**

**Problems with NonlinearFit**

**roots**

**Re: roots**