The theory of dynamical systems has been an area of active research amongst physicists and mathematicians for several decades. Perhaps one of the most interesting objects of study in dynamical systems is the phenomenon of chaos.
Chaos is a phenomenon observed in some deterministic systems where the system’s time-dependent trajectory neither diverges, nor converges to a limit point, nor moves periodically. Rather, the system moves on an orbit which looks random.
A simple model system which evidences chaos is the so-called logistic map:
The idea is to choose a value of
Repeat this process for large numbers of iterations and for many values of
A plot of the fixed point
For example, when
Figure 1: The bifurcation diagram for the logistic map. (From Wikipedia.)
This transition from the period-one orbit to chaos is called the “period-doubling” route to chaos. It attracted the attention of many physicists in the 1970s and 80s since it offered the promise of illustrating a mechanism by which we could understand the development of complicated chaotic phenomena such as turbulence in fluids. Indeed, several physical systems were identified which evidenced a period-doubling route to chaos, including Duffing’s oscillator, a dripping faucet, and some simple electronic circuits. Unfortunately, the larger ambition of finally getting a grasp on turbulence by studying the logistic equation did not pan out. Nonetheless, some very interesting mathematics was discovered in the process.
In the late 1970s, Mitchell Feigenbaum, a mathematician at Los Alamos research laboratory, was playing around with the logistic equation using a hand calculator. He noticed some interesting numbers characterize the period doubling transition to chaos in the logistic map. He observed the following:
Examine the length of each interval in
Figure 2: A close-up of the logistic map’s bifurcation diagram showing how
Examine the width of the opening of each parabolic segment when it hits the value superstable” point, and its value depends upon which map is under consideration. The value is 1/2 for the logistic map.) The ratio of successive widths also converges to a value. This is shown in Figure 3. Feigenbaum called this value
Figure 3: A close-up of the logistic map’s bifurcation diagram showing how
Most importantly, Feigenbaum found that the value of these numbers was independent of the exact map used to generate period-doubling and chaos. Specifically, as long as the the map of the unit interval to itself,
is strictly concave and has a quadratic peak, then the map will evidence a period doubling transition to chaos characterized by numbers delta and alpha having exactly the values shown above.
For example, Therefore, the numbers and are universal numerical constants on equal footing as the commonly known constants and
This raises the question: How to calculate high-precision values for
My goal is to compute this function, and then use it to compute a high-precision value for
I am aware of three earlier, high-precision computations of alpha. The first was performed by Australian mathematician Keith Briggs in his PhD thesis . He computed 576 decimal places (of which 344 were later found to be correct). The second was performed by UK mathematician David Broadhurst in 1999 , and is linked from the OEIS. He computed 1018 digits (after the decimal point). Finally, while I was preparing this post Professor Broadhurst pointed me to a paper by Andrea Molteni  who used a Chebyshev expansion instead of a Taylor’s expansion for 10,000 digits.
I became interested in chaos as an undergraduate, my curiosity including learning about the universal numbers
Decades later, when I became aware of Julia and its capabilities, I realized that it has some features which make it a good tool for computing these numbers. Julia’s important features include:
Support for BigFloat. Prior computations used C or other low-level languages and manually linked to arbitrary-precision floating point libraries. This is complicated and error-prone. Julia makes it easy: Just use big() or BigFloat when declaring your high-precision numbers, and the BigFloat type will propagate through the rest of the calculation.
Autodiff support. Like Briggs and Broadhurst, I employ the multidimensional Newton’s method to find the universal function http://www.juliadiff.org/. For my calculation I used the package ForwardDiff.jl.
With these features in mind, here is the method I used to compute alpha. It is the same used by Briggs and by Broadhurst, except that I implemented the computation in Julia.
Define a function
Then use Newton’s method to solve the system. I use gradient.jl (part of the ForwardDiff.jl package) to compute the gradient of
I wrap Newton’s method with an outer loop which calls it repeatedly. I start by asking for only, say, 10 digits. Newton’s method performs the solve then returns the expansion coefficients used to get the 10 digits. I then call Newton’s method again, but ask for more digits, say 30. I use the previously-found coefficients as the starting point for the new computation. This helps keep Newton’s method from wandering away from the desired solution as the number of requested digits increases. It is also useful when validating my result since I have a sequence of converging values which I can compare to each other.
My code to calculate alpha is available on GitHub for you to peruse, Here is the result of a typical run. As you see, I start by requesting a small number of digits, then walk up the requested precision. The variable
N = 10
N = 30
N = 50
If you compare each successive result to the next one, you will see that the sequence converges – the N=10 result contains 11 converged digits (after the decimal point), the N=30 result contains 36 converged digits, and so on.
Here is my best result so far, 1177 digits. This includes new digits at the end, beyond those computed by Broadhurst.
To check my results, I wrote a Julia script which accepts two files holding computed digits. The script then starts at the beginning of each file, and counts the number of matching characters.
stest = readstring("alpha_me.txt") #reading the test file
sref = readstring("alpha_oeis.txt") #reading the reference file
s = ""
count = 0
for i in length(min(stest,sref)) #iterating to count matches
if(stest[i] == sref[i])
count = count + 1
s = s * string(stest[i])
println("alpha_me.txt has $(length(stest)) digits.") #printing results
println("alpha_oeis.txt has $(length(sref)) digits.")
println("Test file agrees with reference file to $count digits.")
println("Digits of agreement: alpha = $s")
Here’s a run comparing one of my results against the digits reported by Broadhurst:
alpha_me.txt has 2441 digits.
alpha_oeis.txt has 1021 digits.
Test file agrees with reference file to 1020 digits.
Digits of agreement: alpha =
This result matches all of Broadhurst’s digits, giving good confidence that my results are correct. Then, as I compute higher and higher number of digits, I can compare successive files against one another to find out how many digits have converged. As mentioned above, my current best result (above) is 1177 digits.
I am not yet the world’s record holder when it comes to computing Andrea Molteni who computed 10000 digits.
However, his calculation was written in C and linked to MPFR and MPI, so undoubtedly took some time and effort to write and debug. By taking advantage of Julia’s high-level features, I was able to write my program inside of a few hours.
Holmes, P., Whitley, D. “On the attracting set for Duffing’s equation”, Physica D, 111–123 (1983)
A. D’Innocenzo and L. Renna, “Modeling leaky faucet dynamics”, Phys. Rev. E, 55, 6676 (1997).
Paul Lindsay, “Period Doubling and Chaotic Behavior in a Driven Anharmonic Oscillator”, Physical Review Letters, 47 1349 (1981).
Mitchell J. Feigenbaum, “Universal Behavior in Nonlinear Systems,” in Los Alamos Science, Summer 1980, pp. 4-27 (available online at http://permalink.lanl.gov/object/tr?what=info:lanl-repo/lareport/LA-UR-80-5007).
Keith Briggs, “Feigenbaum scaling in discrete dynamical systems” (PhD Thesis, University of Melbourne, 1997) (available online via http://keithbriggs.info/thesis.html.
David Broadhurst, private communication. His results are available online at http://www.plouffe.fr/simon/constants/feigenbaum.txt.
Andrea Molteni, “An Efficient Method for the Computation of the Feigenbaum Constants to High Precision”, https://arxiv.org/pdf/1602.02357.pdf. His results are available online at http://converge.to/feigenbaum/.
Philipp H. W. Hoffmann, “A Hitchhiker’s Guide to Automatic Differentiation”, https://arxiv.org/pdf/1411.0583.pdf.
This is a guest blog post, written and edited by Stuart Brorson, Dept of Mathematics, Northeastern University, and formatted by Julia Computing’s Rajshekar Behar
Need help with Julia?
We also provide training and consulting services
and build open source or proprietary packages
for our customers on a consulting basis. Email us:
Julia Computing was founded by all the creators
of the language to provide commercial support
to Julia users. We are based in Boston, New York,
San Francisco, London and Bangalore with
customers across the world.
© 2016 - 2020 Julia Computing, Inc. All Rights Reserved.