Re: Find many Roots
- To: mathgroup at smc.vnet.net
- Subject: [mg32181] Re: [mg32098] Find many Roots
- From: Daniel Lichtblau <danl at wolfram.com>
- Date: Fri, 4 Jan 2002 05:03:56 -0500 (EST)
- References: <200112240832.DAA15698@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
Yoram Pollack wrote:
>
> Hello,
> I am trying to find all roots of the Equation: Tan[x]==1/x, {x,0,12}].
> What is the best way?
> Thanks
> Yoram
Several nice methods have been posted. I thought I would draw attention
to a similar question from last June.
http://library.wolfram.com/mathgroup/archive/2001/Jun/msg00226.html
The biggest difference is that the present function has poles. One can
instead find zeroes of x*Sin[x]-Cos[x] to avoid this issue.
My own reply to the June post (URL below) gave one method for attacking
such problems, using the Cauchy integral formula to count roots in an
interval, bisection to isolate them, and CIF again to approximate
individually isolated roots.
http://library.wolfram.com/mathgroup/archive/2001/Jun/msg00226.html
Since the code was not too long I just include it below. The actual note
goes into more detail regarding some of the advantages, disadvantages,
and potential extensions of this approach.
epsilon = 10^(-1);
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},
PrecisionGoal->3, MaxRecursion->10];
t2 = NIntegrate[Evaluate[integrand], {z,a+I*offset,a-I*offset},
PrecisionGoal->3, MaxRecursion->10];
t3 = NIntegrate[Evaluate[integrand], {z,b-I*offset,b+I*offset},
PrecisionGoal->3, MaxRecursion->10];
(t1+t2+t3)/(2*Pi*I)]
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]]
]
refineRoots[func_,intervals_,eps_,z_]:=
Map[Re[segmentIntegrate[func,#[[1,1]],#[[1,2]],eps,z,z]]&, intervals]
Now we try it out.
f[x_] := x*Sin[x]-Cos[x]
InputForm[Timing[rootintervals =
isolateRoots[f[z],0,12,z,epsilon]]]
This gives some useless messages regarding loss of precision, followed
by
Out[46]//InputForm=
{0.03999999999999915*Second, {Interval[{0, 3}], Interval[{3, 6}],
Interval[{6, 9}], Interval[{9, 12}]}}
In[47]:= Timing[rts=refineRoots[f[z],rootintervals,epsilon,z]]
Out[47]= {0.02 Second, {0.860334, 3.42562, 6.4373, 9.52933}}
We see that the residuals are not fairly small.
In[53]:= InputForm[f[x] /. x->rts]
Out[53]//InputForm=
{-3.984138952084493*^-8, 3.80665711174899*^-8, 3.887810040570372*^-8,
7.042598237916309*^-8}
We could of course do better either by higher precision quadrature or,
more sensibly, by taking these as starting values for a high
precision/accuracy FindRoot.
Other replies to the June query gave code to implement entirely
different methods e.g. polynomial approximation, interval search, and
maybe others as well; these all could be applied to this example as
well.
Daniel Lichtblau
Wolfram Research