MathGroup Archive 2006

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

Search the Archive

Re: Strange problem with InverseSeries

  • To: mathgroup at smc.vnet.net
  • Subject: [mg64719] Re: [mg64709] Strange problem with InverseSeries
  • From: "Carl K. Woll" <carlw at wolfram.com>
  • Date: Wed, 1 Mar 2006 04:11:38 -0500 (EST)
  • References: <200602280649.BAA17661@smc.vnet.net>
  • Sender: owner-wri-mathgroup at wolfram.com

Dan Goodman wrote:
> Hi all,
> 
> I'm running Mathematica 4 and I have come across a very strange problem 
> with the InverseSeries function. The following works fine, and does 
> exactly what you'd expect.
> 
> p[X_] := -2 + 6*X - 9*X^2 + 8*X^3 - 6*X^4 + 2*X^5 - X^6;
> X0 = X /. FindRoot[p[X] == 2, {X, 2 I}, MaxIterations -> 100];
> s = Series[p[X], {X, X0, 10}];
> InverseSeries[s];
> 
> Clear[c];
> is = InverseSeries[SeriesData[X, X0, c /@ Range[0, 10]]];
> (c[#1] = SeriesCoefficient[s, #1];) & /@ Range[0, 10];
> is;
> 
> The second part just computes the general case of the inverse series, 
> and then substitutes the particular values. The two methods give the 
> same result.
> 
> However, if I replace this polynomial with the following one
> 
> p[X_] := I*(-2 + 7*X - 12*X^2 + 14*X^3 - 10*X^4 + 7*X^5 - 2*X^6 + X^7);
> 
> then InverseSeries[s]; just runs and runs and (as far as I can tell) 
> never stops. The second method works fine, but is very slow 
> (particularly if you try to compute more terms of the series than 10).
> 
> Any ideas why InverseSeries chokes on this second polynomial but not on 
> the first? My problem is that I want to run this calculation for much 
> higher degree polynomials (the ones I'm looking at all have the same 
> problem if the degree is 7 or higher), and compute many more than the 
> first 10 coefficients of the inverse series.
> 
> If someone can think of a clever and speedy workaround that would make 
> me a very happy and contented individual.
> 
> Many thanks,
> 
> Dan Goodman
> University of Warwick

Dan,

I had no problems using InverseSeries on either of your polynomials. 
Which version of Mathematica were you using?

At any rate, if your roots are simple, than it is possible to use newton 
iteration to compute the series expansions of roots. Here is a function 
which carries out this idea:

simpleroot[f_, r_, z_, steps_: 1] :=
  simpleroot[f, Normal[r - f[r]/f'[r] + O[z]^(2(Exponent[r, z]+1))], z,
   steps - 1]
simpleroot[_, r_, _, 0] := r

The function simpleroot solves for the roots of a bivariate function f, 
starting with the zeroth order solution r. The function f is a pure 
function of one of the variables, with the other variable being the 
parameter z. The output, ans, is a polynomial in z, such that f[ans] is 
equal to 0 to some order in z. Each newton iteration doubles the number 
of terms in the output, so a value of 5 for steps will result in an 
order 32 polynomial. In other words, don't choose too large a value.

For your example, we have

p[x_] := I (-2 + 7 x - 12 x^2 + 14 x^3 - 10 x^4 + 7 x^5 - 2 x^6 +
      x^7);

You were basically interested in solving the equation

p[x]-2==y

for x. First we find the zeroth order roots to p[x]-2:

In[24]:=
rts = NSolve[p[x]-2 == 0,x]

Out[24]=
                   -16
{{x -> -5.82867 10    - 1.80194 I},

                    -17
   {x -> -3.09345 10    - 0.445042 I},

                   -16
   {x -> 7.07571 10    + 1.24698 I}, {x -> 0.0962758 + 1.84628 I},

   {x -> 0.355946 - 1.43588 I}, {x -> 0.67077 + 0.890932 I},

   {x -> 0.877008 - 0.301327 I}}

You were interested in the root closest to 2I, which is the 4th root. 
Use simpleroot:

In[25]:=
ans = simpleroot[Function[{x},p[x]-2-y],x/.rts[[4]],y,3]

Out[25]=
(0.0962758 + 1.84628 I) + (0.0158496 + 0.0173101 I) y -

                           -6     2
   (0.00211246 - 6.48792 10   I) y  +

                                  3
   (0.000209137 - 0.000225001 I) y  +

              -6                    4
   (3.49868 10   + 0.0000538197 I) y  -

              -6             -6     5
   (8.05371 10   + 6.70792 10   I) y  +

              -6           -7     6              -7             -7     7
   (2.15628 10   - 2.504 10   I) y  - (2.82683 10   - 3.75039 10   I) y

Let's check:

In[26]:=
p[ans]-2-y+O[y]^10//Chop

Out[26]=
            -6             -6     8
(3.74191 10   + 2.44352 10   I) y  -

              -7            -8     9       10
   (2.16994 10   - 3.0917 10   I) y  + O[y]

You will need to substitute x-2 for y if you want an answer in the same 
form. If your polynomials are algebraic (have integer or even algebraic 
coefficients), then it is possible to compute the coefficients of the 
answer exactly. If your roots aren't simple, i.e., you have repeated 
roots, then you will need to perform additional, newton polygon, steps.

Carl Woll
Wolfram Research


  • Prev by Date: Re: WolframSSH vs mathssh and private key command line option
  • Next by Date: Re: Re: Limit
  • Previous by thread: Re: WolframSSH vs mathssh and private key command line option
  • Next by thread: Re: Re: Limit