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
- roots