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