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
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 page — index — recent changes — DEMZ
This page was written in the "embarrassingly readable" markup language RHTF, and was last updated on 2019 Feb 01.
s.27