mdbtxt1
mdbtxt2
Proceed to Safety

Laurent Series    

Robert P. Munafo, 2000 Apr 19.



In general, a Laurent series is an infinite series, like a Taylor series, whose terms can be used to approximate some function. One important difference from Taylor series is that the Laurent series method can be used for functions that are less well-behaved, including non-differentiable functions.

In Mu-Ency "Laurent series" refers to a specific Laurent Series described by Ewing and Schober. The function it models is the conformal map psi mapping the unit disk onto the Mandelbrot Set.

Ewing and Schober showed that the area of the Mandelbrot Set could be computed according to the formula:

AM = pi . ( 1 - sigma(n=0 to infinity)( n . bn2 ) )

Where bn are the coefficients of the Laurent series. The calculation of the coefficients is discussed below.

Because the bn values are squared, the sum always increases, and therefore the value of the whole expression always decreases, as successive finite sum approximations are computed. Thus, for a finite n, the formula gives an upper bound to the area.

Attempts to estimate the area using the pixel counting method yield estimates between 1.50 and 1.55. When Ewing and Schober evaluated this series to 240,000 terms, it looked as though the series was going to converge somewhere between 1.60 and 1.65. The series converges so slowly that Ewing and Schober did not try to determine the value of the limit any more precisely than this.

The following table shows the convergence behavior of the series. Scott Huddleston compiled most of the values in this table; the last (for 240,000 terms) is from the Ewing and Schober paper. Each entry in the table represents a doubling of the number of terms used to generate the estimate (except the last two, which are approximately double):

number
of terms
area
upper
bound
delta 1.50659 + e / ln(n)
(see below)
128 2.02781 2.06718 o
0.05029
256 1.97752 1.99710 o
0.05001
512 1.92751 1.94260 o
0.03217
1024 1.89534 1.89900 o
0.04073
2048 1.85461 1.86333 o
0.02009
4096 1.83452 1.83360 u
0.02836
8192 1.80616 1.80844 o
0.01980
16384 1.78636 1.78688 o
0.01953
32768 1.76683 1.76820 o
0.01346
65536 1.75337 1.75184 u
0.01490
115232 1.73847 1.73997 o
0.01097
240000 1.72750 1.72615 u

Notice how the first differences go up and down: the sequence converges erratically. Also notice how, even though we keep doubling the number of terms, each doubling brings the area estimate down almost as much as the previous doubling did. That means if we want to use this method to get even two good digits of the Mandelbrot set's area, we'll have to double many more times.

When this was discussed in 1992, many believed that the limit of this series is greater than the limit of the pixel counting method. Ewing and Schober suggested either that the formula based on the Laurent series might be wrong (due most probably to poor understanding of the boundary of the Mandelbrot set and how the conformal map approaches it), or that perhaps the boundary has positive area. (It is known that the boundary's Hausdorff dimension is 2.0, but that doesn't guarantee that it has positive area.)

However, it seems most likely given the first-difference data, that the series does converge to the same value as given by pixel-counting, but that the rate of convergence is extremely slow. A good formula for the difference between a Ewing/Schober estimate and the pixel-counting area is:

difference = e / ln(number-of-terms)

The value of 1.50659 + e/ln(n) is shown in the 4th column of the table above. Each value is marked with an o if it is an overestimate of column 2, (that is, higher than the actual Laurent series sum), and a u if it is an underestimate. As you can see, it starts out consistently higher, but after the 5th row of the table it continually goes back and forth between being an overestimate and being an underestimate. Thus, it appears to be a good statistical model of the behavior of the Laurent series sum.

Based on that, the number of terms needed to get a higher bound which agrees with the actual Mandelbrot set area to within N digits past the decimal point is

number-of-terms = ee 10N

Which gives 6.39 × 1011 terms to get the first digit, and an inconceivable 1.13 × 10118 terms to get the second!

The Coefficients and Details of Computation

This is the definition of the conformal mapping psi:

psi(w) = w + sigma(n = 0 to infinity)(bn / wn)

This expands to:

psi(w) = w - 1 / 2 + 1 / (8 w) - 1 / (4 w2) + 15 / (128 w3) + 0 / w4 - 47 / (1024 w5) - 1 / (16 w6) + 987 / (32768 w7) + 0 / w8 - 3673 / (262144 w9) + 1 / (32 w10) - 61029 / (4194304 w11) . . .

The coefficients of the series are b0 = -1/2, b1 = 1/8, b2 = -1/4, b3 = 15/128, b4 = 0, b5 = -47/1024, etc. The terms bn are the same in the formula for the Mandelbrot set area which was given above:

AM = pi . ( 1 - sigma(n=0 to infinity)( n . bn2 ) )

The terms bn are defined recursively in terms of three other sets of terms called w(n,m), aj, and u(n,k). The sets w and u are infinite 2-dimensional arrays (kind of like Pascal's triangle) and the other two (b and a) are ordinary infinite series. All of their values are rational numbers, with powers of 2 in the denominator. Here are the definitions Ewing and Schober gave for the four sets of terms:

bn = -1/2, if n = 0
otherwise bn = ( - w(n,n+1) ) / n

w(n,m) = 0, if n = 0
otherwise w(n,m) = am-1 + w(n-1,m) + sigma(j = 0 ... m-2)(aj . w(n-1,m-j-1))

aj = u(0,j+1)

u(n,k) = 1, if 2n-1 = k
otherwise u(n,k) = sigma(j = 0 .. k)(u(n-1,j) . u(n-1,k-j)), if 2n-1 > k
otherwise u(n,k) = 0, if 2n+1-1 > k
otherwise u(n,k) = (u(n+1,k) - sigma(j = 1 ... k-1)(u(n,j) . u(n,k-j))) / 2

(These definitions were adapted from the Maple code given by Gerald Edgar and listed below)

Recursion and the Importance of Storing Values

As you can see, the values of w(n,m) and u(n,k) are usually computed in terms of several other w(n,m) and/or u(n,k) values. When this happens, smaller n, m, n and k are used and the recursive process repeats until they all manage to get down to one of the trivial cases (like w(n,m) when n=0).

This means that if you try to write a computer program that computes bn terms, and defines these four sets of terms as functions that call themselves and each other, the program will take exponentially longer as you go to higher values of n. Here is a table showing how many times one of the four functions will be invoked using this "simple recursive" approach:

to compute
bn for
n from 0 to:
function
calls
ratio series sum
0 1 3.1415
15.00
1 15 3.0925
7.13
2 107 2.6998
6.06
3 648 2.5703
5.57
4 3608 2.5703
5.28
5 19063 2.5372
5.12
6 97570 2.4636
5.01
7 488768 2.4437
4.98
8 2411908 2.4437

Clearly that won't work — each time you add a term it's taking about 4 times as long as all the previous terms put together. The exponential rate is decreasing, but not very quickly. The first entry in the table above is for 128 terms. Even if the average exponential rate for the first 128 terms turned out to be as low as 3.00, the calculation would require something like 1061 steps!

This is a common problem in numerical mathematics, and the solution is to tell the computer to remember each value of w(n,m), aj, and u(n,k) the first time it is computed. The recursive function will not have to be evaluated the second and subsequent times, saving many steps. Here is the same table for the first 8 terms, using this improvement:

to compute
bn for
n from 0 to:
function
calls
ratio series sum
0 1 3.1415
11.00
1 11 3.0925
1.82
2 20 2.6998
1.55
3 31 2.5703
1.42
4 44 2.5703
1.34
5 59 2.5372
1.37
6 81 2.4636
1.25
7 101 2.4437
1.22
8 123 2.4437
16 381 2.2730
32 1295 2.1830
72 5857 2.0928

Now, to get n terms, you only have to perform a little more than n2 steps. That is certainly a significant improvement over the exponential pattern. However, there is a substantial drawback in the fact that you have to store the terms. The w(n,m) and u(n,k) terms create the greatest impact, because they are square arrays: to compute n terms, you end up computing n2 values of w(n,m) and u(n,k). That means that when they computed 240,000 terms of the series, Ewing and Schober's computer had to store about 57,600,000,000 terms of w(n,m) and an equal number of u(n,k).

Given this, and the fact that the series converges so slowly, it is clear that the Ewing-Schober sequence has little practical value as a numerical algorithm for computing the area of the Mandelbrot set. In the example above where we considered computing 2 digits past the decimal point, it took 10118 terms, and that would mean storing 10236 values in memory, give or take a few.


Gerald Edgar has provided Maple code to compute the coefficients, using two different methods. The first method is described as being equivalent to the method of Ewing and Schober:

u := proc(n,k) local j; option remember; if 2^n-1=k then 1; elif 2^n-1>k then sum('u(n-1,j)*u(n-1,k-j)',j=0..k); elif 2^(n+1)-1>k then 0; else (u(n+1,k)-sum('u(n,j)*u(n,k-j)',j=1..k-1))/2; fi; end; a := proc(j) option remember; u(0,j+1) end; w := proc(n,m) local j; option remember; if n=0 then 0; else a(m-1)+w(n-1,m)+sum('a(j)*w(n-1,m-j-1)',j=0..m-2); fi; end; b := proc(m) option remember; if m=0 then -1/2; else -w(m,m+1)/m; fi; end;

Adam Majewski provides this translation into maxima:

u[n,k] := block( [j], if 2^n-1=k then 1 else if 2^n-1>k then sum(u[n-1,j]*u[n-1,k-j],j,0,k) else if 2^(n+1)-1>k then 0 else (u[n+1,k]-sum(u[n,j]*u[n,k-j],j,1,k-1))/2 );    a[j] := u[0,j+1];    w[n,m] := block( [j], if n=0 then 0 else a[m-1]+w[n-1,m]+sum(a[j]*w[n-1,m-j-1],j,0,m-2) );    b(m) := if m=0 then -1/2 else -w[m,m+1]/m ;

Gerald Edgar's second method (also in Maple) is described as being more efficient:

u := proc(n,k) local j; option remember; if 2^n-1=k then 1; elif 2^n-1>k then sum('u(n-1,j)*u(n-1,k-j)',j=0..k); elif 2^(n+1)-1>k then 0; else (u(n+1,k)-sum('u(n,j)*u(n,k-j)',j=1..k-1))/2; fi; end; a := proc(j) option remember; u(0,j+1) end; v := proc(k,j) local l; option remember; if k=1 then -a(j-2)-sum('v(1,l)*a(j-l-1)',l=2..j-1); else v(1,j-k+1) + v(k-1,j-1)+sum('v(1,l)*v(k-1,j-l)',l=2..j-k); fi; end; b := proc(j) local k; option remember; if j=0 then -a(0); else -a(j)-sum('v(k,j)*b(k)',k=1..j-1); fi; end;

Adam Majewski also provides this numerical version in maxima:

betaF[n,m]:=block ( [nnn:2^(n+1)-1], if m=0 then 1.0 else if ((n>0) and (m < nnn)) then 0.0 else (betaF[n+1,m]- sum(betaF[n,k]*betaF[n,m-k], k, nnn, m-nnn) -betaF[0,m-nnn])/2.0 );    b(m):=betaF[0,m+1];


References

Ewing and Schober paper


Acknowledgments

Area estimates table (Laurent Series method):

The 240000 estimate is from the Ewing & Schober article.

The rest are from: Scott Huddleston (scott at crl labs tek com)

Opinions on the discrepancy: Stan Isaacs (isaacs at hpcc01 hp com)

Ewing & Schober reference: Jeff Shallit (shallit at graceland waterloo edu)

Laurent series formula and Maple code: Gerald Edgar (edgar at math ohio-state edu)


Here is Scott Huddleston's complete table of area estimates: # terms: area estimate 72: area < 2.09288 128: area < 2.02781 180: area < 2.01237 256: area < 1.97752 360: area < 1.94961 512: area < 1.92751 720: area < 1.91255 1024: area < 1.89534 1440: area < 1.87172 2048: area < 1.85461 2880: area < 1.84576 4096: area < 1.83452 5760: area < 1.81649 8192: area < 1.80616 11520: area < 1.79642 16384: area < 1.78636 23040: area < 1.7747 32768: area < 1.76683 46080: area < 1.75936 65536: area < 1.75337 92160: area < 1.74321 115232: area < 1.73847


Here is the beginning of the Laurent series, in both reduced and unreduced forms:

n reduced unreduced
0 -1 / 2 -1 / 2
1 1 / 8 1 / 8
2 -1 / 4 -8 / 32
3 15 / 128 15 / 128
4 0 / 1 0 / 512
5 -47 / 1024 -94 / 2048
6 -1 / 16 -512 / 8192
7 987 / 32768 987 / 32768
8 0 / 1 0 / 131072
9 -3673 / 262144 -7346 / 524288
10 1 / 32 65536 / 2097152
11 -61029 / 4194304 -122058 / 8388608
12 0 / 1 0 / 33554432
13 -689455 / 33554432 -2757820 / 134217728
14 -21 / 512 -22020096 / 536870912
15 59250963 / 2147483648
16 0
17 -164712949 / 17179869184
18 39 / 2048
19 -2402805839 / 274877906944
20 -1 / 64
21 -4850812329 / 2199023255552
22 29 / 2048
23 -18151141041 / 70368744177664
24 0
25 3534139462275 / 562949953421312
26 -1039 / 131072
27 -22045971176589 / 9007199254740992
28 -1 / 256
29 -750527255965871 / 72057594037927936
30 -4579 / 524288
31 54146872254247683 / 9223372036854775808
32 0
33 -155379776183158669 / 73786976294838206464
34 2851 / 1048576
35 -6051993294029466699 / 1180591620717411303424
36 -1 / 1024
37 7704579806709870203 / 9444732965739290427392
38 92051 / 16777216
39 -403307733528668035403 / 302231454903657293676544
40 0
41 1650116480759617184697 / 2417851639229258349412352
42 -229813 / 67108864
43 36124726440442241978477 / 38685626227668133590597632
44 -41 / 4096
45 -225851495844149964787753 / 309485009821345068724781056
46 564373 / 67108864
47 -35761228458796476847725737 / 19807040628566084398385987584
48 0
49 362376876750551361794705823 / 158456325028528675187087900672
50 -29407003 / 8589934592
51 -6510398483578238274151194427 / 2535301200456458802993406410752
52 33 / 8192
53 74815618913797220433481657203 / 20282409603651670423947251286016
54 -30057875 / 34359738368
55 -698617278028915809388280344009 / 649037107316853453566312041152512
56 0
57 -8675905413734991085610532783493 / 5192296858534827628530496329220096
58 -27868893 / 68719476736
59 -375687870961637050293461860951517 / 83076749736557242056487941267521536
60 1 / 4096
61 -1418434432207399687114226995905967 / 664613997892457936451903530140172288
62 -11847286243 / 1099511627776
63 1084116104452462070609082665064238307 / 170141183460469231731687303715884105728
64 0




From the Mandelbrot Set Glossary and Encyclopedia, by Robert Munafo, (c) 1987-2024.

Mu-ency main pageindexrecent changesDEMZ


Robert Munafo's home pages on AWS    © 1996-2024 Robert P. Munafo.    about    contact
This work is licensed under a Creative Commons Attribution-NonCommercial 4.0 International License. Details here.

This page was written in the "embarrassingly readable" markup language RHTF, and was last updated on 2019 Feb 01. s.27