roots
maarten.v...@icos.be
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
Andrzej Kozlowski
roots will all be between -20<x<20 (since Abs[Cos[x]]<1). So start with
In[5]:=
Plot[0.05 x +Cos[x],{x,-20,20}]
You will see all 13-teen roots. Now click on the graphic to select it.
Holding the Command key down on the Mac (I guess the Control key under
Windows?) click on the individual points where the curve intersects the x
axis and Copy. Paste the result into the FindRoot command (removing the
y-value) as below. You will get pretty good approximations to all the roots.
Here are just the first four values obtained in this way:
In[7]:=
FindRoot[0.05*x + Cos[x] == 0, {x, -18.8966}]
Out[7]=
{x -> -19.1433}
In[8]:=
FindRoot[0.05*x + Cos[x] == 0, {x, {-18.394}}]
Out[8]=
{x -> -18.4538}
In[9]:=
FindRoot[0.05*x + Cos[x] == 0, {x, -13.4688}]
Out[9]=
{x -> -13.4028}
In[10]:=
FindRoot[0.05*x + Cos[x] == 0, {x, -11.6596}]
Out[10]=
{x -> -11.6152}
--
Andrzej Kozlowski
Toyama International University
JAPAN
http://platon.c.u-tokyo.ac.jp/andrzej/
http://sigma.tuins.ac.jp/~andrzej/
on 01.6.16 3:47 PM, maarten.v...@icos.be at
BobH...@aol.com
In a message dated 2001/6/16 2:54:32 AM, maarten.v...@icos.be writes:
>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.
>
Here is an approach
findAllRoots[eqn_, {x_Symbol, sp_?NumericQ},
Interval[{xmin_?NumericQ, xmax_?NumericQ}], n_Integer?Positive,
opts___?OptionQ] :=
Module[
{soln = {FindRoot[eqn, {x, sp}, opts]}, cr},
While[Length[soln] < n,
cr = FindRoot[eqn, {x, Random[Real, {xmin, xmax}]}, opts];
If[xmin <= (x/.cr) <= xmax,
soln = Union[Append[soln, cr],
SameTest -> (Abs[(x/.#1)-(x/.#2)]<10^-5&)]]];
Sort[soln]];
expr = 0.05 x +Cos[x];
nbrRoots = 13;
ans = findAllRoots[expr==0, {x, 0}, Interval[{-20, 18}], nbrRoots];
The maximum absolute error is
Max[Abs[expr /. ans]]
4.5600301290527057*^-7
Plot[expr, {x, -26, 23}, PlotStyle -> RGBColor[0, 0, 1],
Epilog -> {AbsolutePointSize[4], RGBColor[1, 0, 0],
Point[{x,expr}] /.ans}];
Bob Hanlon
Chantilly, VA USA
David Park
In a problem like this, FindRoot needs some help.
f[x_] := 0.05x + Cos[x]
Look for the turning points in the function.
f'[x]
0.05 - Sin[x]
The turning points will occur at ArcSin[0.05] +/- 2n Pi and
(Pi-ArcSin[0.05])+/- 2n Pi. This generates the turning points over a
somewhat wider range than we need.
tpts = Sort[Flatten[Table[{ArcSin[0.05] + 2*Pi*n,
Pi - ArcSin[0.05] + 2*Pi*n}, {n, -4, 3}]]];
Between successive turning points we have monotonic functions. If we use the
turning points as limits and start at the midpoint, FindRoot just might get
the root. The following generates the second argument for FindRoot using
this approach.
Partition[tpts, 2, 1]
rootspecs = {x, Plus @@ #/2, First[#], Last[#]} & /@ %
This is the second item on a list of 15 generated rootspecs.
{x, -20.4204, -22.0412, -18.7995}
I then Map FindRoot onto this list of rootspecs. I used the Check command to
weed out cases where FindRoot went outside the specified domain. You obtain
error messages for those cases but no returned value.
rootsols = Check[FindRoot[f[x], #], Unevaluated[Sequence[]]] & /@ rootspecs
{{x -> -19.14330451998876}, {x -> -18.453754054424195},
{x -> -13.402772345741981},
{x -> -11.615238480230595},
{x -> -7.471140913774326}, {x -> -4.963166669457849},
{x -> -1.4959299134300144},
{x -> 1.6535692756905693}, {x -> 4.486156653004949},
{x -> 8.280873579694441}, {x -> 10.446026873732087},
{x -> 14.98402202417216}, {x -> 16.323959951197935}}
We could also have weeded our the domains with no roots by checking if f
evaluated at the turning points changed signs.
It would be nice if FindRoot had an option or mode such that it wouldn't
quit when it goes outside the specified domain, but used a bisection step
and resumed from the midpoint of the resulting domain. I think this would
generally be of value when we were dealing with monotonic functions in the
specified range.
David Park
dj...@earthlink.net
http://home.earthlink.net/~djmp/
> From: maarten.v...@icos.be [mailto:maarten.v...@icos.be]
>
> 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.
>
Hermann Meier
> Maarten van der Burgt wrote:
> What is the easiest solution to finding the 13 (real) roots of 0.05 x +
> Cos[x] == 0. (....)
<< NumericalMath`IntervalRoots`
IntervalNewton[0.05 x + Cos[x], x, Interval[{-30., 30.}], 10^-10.,
MaxRecursion -> 100]//InputForm
gives
Interval[{-19.14330451808022, -19.14330451805268},
{-18.45375405471094, -18.453754054710902},
{-13.402771471945435, -13.402771471945416},
{-11.615238614490243, -11.615238614490224},
{-7.471140913733229, -7.471140913733216},
{-4.963167688817952, -4.9631676888179435},
{-1.4959299132767379, -1.4959299132677975},
{1.6535692761064251, 1.6535692761313203},
{4.4861562872107585, 4.4861562872107665},
{8.280873580136642, 8.28087358013666},
{10.446026863672607, 10.44602686367263},
{14.984022024364082, 14.984022024364215},
{16.323959949381134, 16.323959949381166}]
Hermann Meier
hme...@webshuttle.ch
Adam Strzebonski
Cos[x] with initial terms of its Taylor series (we don't
get a proof that there are only 13 solutions, but we get
13 good approximations of solutions).
In[1]:= poly=x/20+Normal[Series[Cos[x], {x, 0, 100}]];
We need to find approximations of real roots of poly.
I don't want to use NSolve here, because it tries to find
all complex roots. For high degree polynomials, if the number
of real roots is much smaller than the degree, it is faster
to isolate the real roots using exact methods than to find
all the complex roots numerically.
We count the real roots first, so that we don't create any
nonreal Root objects: the exact isolation of nonreal roots
is slow.
In[2]:= <<Algebra`RootIsolation`
In[3]:= CountRoots[poly, {x, -Infinity, Infinity}]//Timing
Out[3]= {4.58 Second, 14}
The default root isolation method is based on computing numerical
approximations of all complex roots, so we need to change it.
In[4]:= SetOptions[Root, ExactRootIsolation->True];
Now we compute the real roots of poly,
In[5]:= rts=Table[Root[poly, i], {i, 14}];//Timing
Out[5]= {2.55 Second, Null}
find their numerical approximations,
In[6]:= nrs=N[rts, 25];//Timing
Out[6]= {4.09 Second, Null}
and check which of them satisfy the equation x/20+Cos[x]==0
In[7]:= #/20+Cos[#]&/@nrs
-24 -24 -24
-24
Out[7]= {-1.65832863853270619235960, 0. 10 , 0. 10 , 0. 10 , 0.
10 ,
-24 -24 -25 -25 -24 -24
-24
> 0. 10 , 0. 10 , 0. 10 , 0. 10 , 0. 10 , 0. 10 , 0. 10 ,
-24 -24
> 0. 10 , 0. 10 }
Here are the 13 roots of x/20+Cos[x]==0
In[8]:= Drop[nrs, 1]//InputForm
Out[8]//InputForm=
{-19.14330451806538667880986570564361032264`25,
-18.45375405471092374532000068940118616129`25,
-13.40277147194542547333729073937332385855`25,
-11.61523861449023313015711996574049279572`25,
-7.47114091373322226214343658073674762609`25,
-4.96316768881794772985279294265779822657`25,
-1.49592991327175815500746007291204533048`25,
1.65356927612222660064136483786947282365`25,
4.48615628721076275861724952418303472826`25,
8.28087358013665167538895805572450457893`25,
10.44602686367261821801301652620945695607`25,
14.98402202436415195182719178207711523666`25,
16.32395994938114931284258124546947389267`25}
Best Regards,
Adam Strzebonski
Wolfram Research
Andrzej Kozlowski
In[2]:=
<<Algebra`RootIsolation`
In[3]:=
CountRoots[poly, {x, -20,20}]//Timing
Out[3]=
{4.98333 Second,13}
The roots can also be found by using InequalitySolve, though this seems to
take longer than your method:
In[5]:=
Timing[N[InequalitySolve[poly == 0 && -20 < x < 20, x]]]
Out[5]=
{68.8833 Second, x == -19.1433 || x == -18.4538 ||
x == -13.4028 || x == -11.6152 || x == -7.47114 ||
x == -4.96317 || x == -1.49593 || x == 1.65357 ||
x == 4.48616 || x == 8.28087 || x == 10.446 ||
x == 14.984 || x == 16.324}
Of course one would still need to verify that they do approximate the roots
of the original equation.
--
Andrzej Kozlowski
Toyama International University
JAPAN
http://platon.c.u-tokyo.ac.jp/andrzej/
http://sigma.tuins.ac.jp/~andrzej/
Daniel Lichtblau
>
> 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
Alan Mason
ZerosAll[cos[x] + 0.05 x, x, -30, 30, 0.08], where mesh size = 0.08 here is
arbitrary. Make changes to taste. The Zeros routine uses the so-called
"Illinois algorithm" and works quite well.
Code (uses procedural program style) follows:
Zeros[expr_, var_,range_List, prec_: 16, eps_:10^-12] :=
Module[{x0, x1, y0, y1,xold, yold,x, y, h, f,func, ran = range}, i=0;
f[arg_] := N[expr //. var->arg, prec];
x0 = ran[[1]]; x1 = ran[[2]];
y0 = f[x0]; y1 = f[x1];
If[y0 y1 > 0, Print["y0 is ", y0];
Print["y1 is ", y1];
Print["Try another initial range"]; Return[]];
While[Abs[x0-x1] > eps,i++;
If[i > 69,Print["MaxIter = 70 Exceeded"]; Return[]];
If[y0 y1 <= 0, h = -y1 (x1-x0)/(y1-y0);
xold = x0;x0=x1; yold =
y0;
y0=y1;x1+= h;y1=f[x1],
While[y0 y1 > 0, i++;
If[i > 69,Print["MaxIter = 70 Exceeded"]; Return[]];
yold=yold/2;
h=-y1 (x1-xold)/(y1-yold);
x0 = x1;x1+=h; y=f[x1]; y0=y1; y1=y]]
];
{N[x0, prec], N[y0, prec]}]
ZerosAll[expr_, var_, lower_, upper_, mesh_, prec_:16, eps_:10^-12] :=
Module[{f, sgn=1, i,L={}, n = Floor[(upper-lower)/mesh]},
f[arg_] := N[expr //. var->arg, prec];
For[i=0, i< n, i++,sgn =Sign[f[lower + i* mesh] f[lower + (i+1)*mesh]];
If[sgn < 0,
L = Append[L,
Zeros[expr, var, {lower+i*mesh, lower+(i+1)*mesh}, prec, eps]]];
];
Return[L];
]
roots = First[#]& /@ L;
errors = Last[#]& /@ L;
<maarten.v...@icos.be> wrote in message
news:9gevki$7k8$1...@smc.vnet.net...