Re: Newton's method.
- To: mathgroup at smc.vnet.net
- Subject: [mg5237] Re: Newton's method.
- From: danl (Daniel Lichtblau)
- Date: Fri, 15 Nov 1996 03:33:52 -0500
- Organization: Wolfram Research, Inc.
- Sender: owner-wri-mathgroup at wolfram.com
In article <5619te$2c7 at dragonfly.wolfram.com> Gandalf <Gandalf at Empire.Net>
writes:
> Can someone provide a little help on this problem?:
>
> Suppose I wanted to write a sequence to perform the Newton inversion of
> a polynomial which inverts a series 1 + a(1)x + a(2)x^2+... through a
> particular degree. The sequence would look like the following:
>
> (* Newton's inversion of a polynomial *)
> guess = 1;
>
> precision = 1;
>
> While[precision < degree,
> precision *=2;
> If[precision > degree, precision = degree];
> temp = iterate[poly, precision - 1];
> guess = guess + guess*(1 - temp*guess); (*The Newton iteration*)
> guess = iterate[guess, precision -1];
> ];
>
> (* Now guess is 1 / poly through degree of 'degree' *)
>
> ----------------
> iterate would be a function, iterate[f, n] that returns a polynomial f
> taken through degree n. This sequence effectively doubles the precision
> each pass through it. I'm a little stumped on how to implement the
> iterate funtion. Any help greatly appreciated. Thanks.
>
This is along the lines of what you seek. You'll need to make some
adjustments for notational differences. I also did not isolate the
'iterate' as a separate routine, and furthermore in effect assumed the
degree we want is a power of two.
seriesInvert[poly_, degree_Integer, x_] := Module[
vars},
vars][[1]];
example:
In[120]:= poly = 1 + 2*x;
In[121]:= seriesInvert[poly, 16, x]
2 3 4 5 6 7 8
Out[121]= 1 - 2 x + 4 x - 8 x + 16 x - 32 x + 64 x - 128 x + 256 x
-
9 10 11 12 13 14
> 512 x + 1024 x - 2048 x + 4096 x - 8192 x + 16384 x -
15
> 32768 x
I am certain that one can write code that is cleaner and more efficient.
But basically I'd expect the time to roughly a-bit-worse-than-double for
each doubling of degree. This is because the expansion, coefficient list
extraction, and solving each deal with roughly twice as many terms, and
the complexity of these gets a bit worse in each iteration. Below are
timings for version 3.0, on a Pentium Pro 200.
In[173]:= Table[Timing[seriesInvert[poly,2^j, x];], {j,1,8}]
Out[173]= {{0.01 Second, Null}, {0.03 Second, Null}, {0.05 Second, Null},
> {0.12 Second, Null}, {0.25 Second, Null}, {0.53 Second, Null},
> {1.19 Second, Null}, {2.89 Second, Null}}
Daniel Lichtblau
Wolfram Research, Inc.
danl at wolfram.com