ALL roots of non-polynomial equation
DrBob
roots but never changes sign!
Or Sin[1/x], which has infinitely many roots converging on 0, but no
limit for the function at 0.
Or... still worse... 1 + Sin[1/x], which has BOTH problems.
Bobby
-----Original Message-----
From: Andrzej Kozlowski [mailto:and...@platon.c.u-tokyo.ac.jp]
Subject: ALL roots of non-polynomial equation
In your example, yes. Here is one way (adapted from a slightly different
problem in Stan Wagon's "Mathematica in Action")
We make use of Mathematica's ability to plot graphs:
In[1]:=
g = Plot[Sin[x], {x, 0.1, 10.1*Pi}, DisplayFunction ->
Identity];
We make a list of all the coordinates of the points represented on the
graph.
In[2]:=
points = Cases[g, Line[x_] -> x, Infinity][[1]];
We make a list of the signs of the y values:
In[3]:=
signs = Sign /@ Transpose[Cases[g, Line[x_] -> x, Infinity][[
1]]][[2]];
We find the points where the sign changes:
In[4]:=
positions = Position[Rest[signs]*Rest[RotateRight[signs]],
-1]
Out[4]=
{{27}, {51}, {74}, {101}, {126}, {149}, {177}, {200}, {226},
{252}}
We make a list of starting points for FindRoot:
In[5]:=
starts = First[Transpose[Extract[points, positions]]]
Out[5]=
{2.7825096162536145, 6.080185995733974, 8.787418231655966,
12.198138489619575, 15.464841498197309, 18.61672099859868,
21.92859710988888, 24.46767425065356, 27.840417480532142,
31.139545383515845}
We find the roots:
In[6]:=
(FindRoot[Sin[x] == 0, {x, #1}, WorkingPrecision ->
20] & ) /@ starts
Out[6]=
{{x -> 3.141592653589793238462643383255068`20},
{x -> 6.283185307179586476925286766538051`20},
{x -> 9.424777960769379715387930149825109`20},
{x -> 12.566370614359172953850573533079026`20},
{x -> 15.707963267948966192313216916378673`20},
{x -> 18.849555921538759430775860299681079`20},
{x -> 21.991148575128552669238503682979946`20},
{x -> 25.132741228718345907701147066183302`20},
{x -> 28.274333882308139146163790449476032`20},
{x -> 31.415926535897932384626433832775678`20}}
This question has been asked frequently so you can find various
approaches, including this one, in the archives. Of course there is no
guarantee. For very complex functions you may well miss some roots. The
situation can become a lot more complicated if your equation has
multiple roots.
Andrzej
On Thursday, August 8, 2002, at 07:06 PM, Mihajlo Vanevic wrote:
>
> Can Mathematica find (localize) ALL roots of non-polynomial equation
>
> eq[x]==0
>
> on a given segment x \in [a,b], a,b=Real??
>
> (for example Sin[x]==0, for 0.1<x<10.1 Pi )
>
>
>
>
>
>
>
DrBob
changes, to get better initial starts for FindRoot:
g = Plot[Sin[x], {x, 0.1, 10.1*Pi}, DisplayFunction -> Identity];
points = First@Cases[g, Line[x_] -> x, Infinity];
signs = Sign /@ points[[All, 2]];
positions =
Union[#, # + 1] &@Flatten@Position[Rest[signs*RotateRight@signs],
-1]
starts = 1/2Plus @@@ Partition[points[[positions, 1]], 2]
x /. (FindRoot[Sin[x] == 0, {x, #1}] &) /@ starts
{27, 28, 51, 52, 74, 75, 101, 102, 126, 127, 149, 150,
177, 178, 200, 201, 226, 227, 252, 253}
{3.09198, 6.26091, 9.13801, 12.5265, 15.7722, 18.9688, 22.0933, 24.8161,
28.1494, 31.4348}
{3.14159, 6.28319, 9.42478, 12.5664, 15.708, 18.8496, 21.9911, 25.1327,
28.2743, 31.4159}
Bobby Treat
Andrzej Kozlowski
be the one posted some time ago by Adam Strzebonski
In[1]:=
f[x_] = Normal[Sin[10.1*Pi*x] + O[x]^100];
In[2]:=
10.1*Pi*Select[x /. NSolve[f[x] == 0, x],
Im[#1] == 0 && 0.1/(Pi*10.1) <= #1 <= 1 & ]
Out[2]=
{31.41492328081239, 28.274379350873716, 25.132740219910605,
21.991148587034512, 18.84955591555759, 15.707963269122267,
12.56637061424962, 9.424777960775808, 6.283185307178523,
3.141592653590188}
The problem is of course that it is difficult to know how long a Taylor
series to take, whehter one has found all the roots and how accurate the
answers are, although with careful analysis this may be a useful
approach.
Andrzej Kozlowski
Andrzej Kozlowski
Toyama International University
JAPAN
http://platon.c.u-tokyo.ac.jp/andrzej/
DrBob
f[x_] = Normal[Sin[x] + O[x]^100];
Plot[f[x] - Sin[x], {x, 0.1, 10Pi}, PlotRange -> All];
Select[x /. NSolve[f[x] == 0, x], Im[#1] == 0 && 0.1 ? #1 ? 10.1Pi &]
{3.14159, 6.28319, 9.42478, 12.5664, 15.708, 18.8496, 21.9911, 25.1327,
\
28.2744, 31.4156}
The Plot should tell us whether we have enough terms in the Series.
David Park
Ted Ersek has a package on MathSource called Rootsearch that is pretty
decent at finding roots. Here is how it works on your examples.
<< Enhancements`Rootsearch`
RootSearch[Sin[x] + 1 == 0, {x, -30, 30}]
{{x -> -26.7035}, {x -> -20.4204}, {x -> -14.1372}, {x -> -7.85398},
{x -> -1.5708}, {x -> 4.71239}, {x -> 10.9956}, {x -> 17.2788}, {x ->
23.5619}, {x -> 29.8451}}
The last example is of course more difficult because there are an infinite
number of roots clustered around zero. However, if we want to find the roots
in a given range:
roots = RootSearch[1 + Sin[1/x] == 0, {x, 0.005, 0.01}]
Length[roots]
{{x -> 0.00501275}, {x -> 0.00517577}, {x -> 0.00534975}, {x ->
0.00553582}, {x -> 0.00573531}, {x -> 0.00594972}, {x ->
0.00618077}, {x -> 0.0064305}, {x -> 0.00670126}, {x ->
0.00699582}, {x -> 0.00731747}, {x -> 0.00767012}, {x ->
0.00805848}, {x -> 0.00848826}, {x -> 0.00896648}, {x -> 0.00950179}}
16
1 + Sin[1/x] /. roots
{0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}
Plot[1 + Sin[1/x], {x, 0.005, 0.01}];
Not too bad.
David Park
dj...@earthlink.net
http://home.earthlink.net/~djmp/
From: DrBob [mailto:maj...@cox-internet.com]
Good luck with a function like 1 + Sin[x], which has infinitely many
roots but never changes sign!
Or Sin[1/x], which has infinitely many roots converging on 0, but no
limit for the function at 0.
Or... still worse... 1 + Sin[1/x], which has BOTH problems.
Bobby
Selwyn Hollis
NumericalMath`IntervalRoots` has exactly what we need.
<<NumericalMath`IntervalRoots`
In: IntervalNewton[Sin[3*x2] - 0.2*x, x, Interval[{0, 2}],
10^(-12), MaxRecursion -> 20]
Out: Interval[{0, 1.337913459004121`*^-18}, {0.06666864225044185`,
0.06666864225044247`}, {0.9903214888289147`,
0.9903214888289159`}, {1.4814395017188728`,
1.4814395017188753`}, {1.7387361933371692`, 1.7387361933371717`}]
In: First /@ List @@ %
Out: {0, 0.0666686, 0.990321, 1.48144, 1.73874}
---------
Selwyn Hollis
Ersek, Ted R
"Looks like we've been reinventing the wheel, folks. The package
----------
Well that isn't entirely true. The above package only works if the equation
can be evaluated using Mathematica's interval arithmetic which does not
support the evaluation of special functions (BesselJ, Beta, Zeta, ....).
---------
Besides that after Interval[{2, 5}] is substituted for each occurrence of
(x) in an expression we get the largest possible interval that can result if
each instance of Interval[{2, 5}] were an independent value between 2 and
5. This sort of interval arithmetic can be helpful in some situations, but
more often one wants to treat each occurrence of Interval[{2,5}] as the same
value in the given interval. And I expect this will prevent the functions
in the (NumericalMath`IntervalRoots`) package from converging in some
cases. For an example of where Mathematica's interval arithmetic gives
misleading results consider the following:
In[1]:=
Clear[x];
poly[x_]=ChebyshevT[12,x];
poly[Interval[{-1,1}]]
Out[3]=
Interval[{-9799,9801}]
We get a very large interval above even though it's well known that
-1 <= ChebyshevT[n, x] <= 1 ( for -1<x<1; Integer n )
I will point out that the NumericalMath`NMinimize` package in Mathematica
4.2 will allow one to numerically find a global minimum over a specified
interval (even in multiple dimensions).
-------
My RootSearch package on MathSource does a very good job of finding all
roots in a finite real interval of an equation in one variable. I have yet
to find another way to easily solve all the problems in the notebook that
goes with the package.
-------
Regards,
Ted Ersek