Services & ResourcesWolfram ForumsMathGroup Archive

[Date Index] [Thread Index] [Author Index]

Re: sorting list of roots af a transcendental function

• To: mathgroup at smc.vnet.net
• Subject: [mg65276] Re: [mg65260] sorting list of roots af a transcendental function
• From: "Carl K. Woll" <carlw at wolfram.com>
• Date: Thu, 23 Mar 2006 06:58:42 -0500 (EST)
• References: <200603221113.GAA10232@smc.vnet.net>
• Sender: owner-wri-mathgroup at wolfram.com

Dule wrote:
> Dear group,
> 
> for calculating a model i need values for x which are given by the 
> transcendental function Cot[x] == x/a - a/(4*x). a is a parameter 0<a<200.
> i obtained the roots with Table and FindRoot:
> Table[FindRoot[Cot[x] == x/a - a/(4*x), {x, i}], {i, 1, 50}]]
> 
> I have two questions:
> 1. Is there a better way to do this?
> 2. How can i construct a list, where the values for x, which appear 
> multiple are dropped?
> 
> Thanks!

One method to find subintervals of an interval guaranteed to contain all 
of the roots of a function f is to use the package IntervalRoots. With 
this package, f must be able to accept an Interval object as an 
argument, which is satisfied with your example. So, load the package:

Needs["NumericalMath`IntervalRoots`"]

We'll use IntervalBisection because your function is discontinuous. With 
IntervalNewton, the roots are found more quickly, but the 
discontinuities may cause some roots to be missed. The default value for 
MaxRecursion is a bit low, so we'll increase it.

In[10]:=
With[{a = 10}, IntervalBisection[Cot[x] - x/a + a/(4*x), x,
      Interval[{0, 50}], 0.1, MaxRecursion -> 30]]

Out[10]=
Interval[{2.24609, 2.34375}, {4.6875, 4.78516}, {7.42188, 7.51953},

   {10.2539, 10.3516}, {13.2812, 13.3789}, {16.2109, 16.3086},

   {19.3359, 19.4336}, {22.3633, 22.4609}, {25.4883, 25.5859},

   {28.6133, 28.7109}, {31.6406, 31.7383}, {34.7656, 34.8633},

   {37.8906, 37.9883}, {41.0156, 41.1133}, {44.1406, 44.2383},

   {47.2656, 47.3633}]

You can refine the search for roots by choosing a smaller eps than 0.1, 
or you can use FindRoot with these smaller intervals. Here we repeat 
with an eps of 10^-6:

In[11]:=
With[{a = 10}, IntervalBisection[Cot[x] - x/a + a/(4*x), x,
      Interval[{0, 50}], 10^-6, MaxRecursion -> 30]]

Out[11]=
Interval[{2.28445, 2.28445}, {4.76129, 4.76129}, {7.46368, 7.46368},

   {10.3266, 10.3266}, {13.2862, 13.2862}, {16.3031, 16.3031},

   {19.3552, 19.3552}, {22.4298, 22.4298}, {25.5197, 25.5197},

   {28.6202, 28.6202}, {31.7285, 31.7285}, {34.8426, 34.8426},

   {37.961, 37.961}, {41.0829, 41.0829}, {44.2075, 44.2075},

   {47.3344, 47.3344}]

Let me repeat the advantage of IntervalBisection. All of the roots to 
your equation between 0 and 50 are guaranteed to lie in one of the above 
intervals. Another nice thing with IntervalBisection is that for 
functions with a finite number of roots, your starting interval can be 
Interval[{-Infinity,Infinity}]. The disadvantage, of course, is that 
most special functions cannot accept Interval objects as input.

Carl Woll
Wolfram Research

• References:
  • sorting list of roots af a transcendental function
    • From: Dule <dule23@gmx.de>