2026-06-08 04:14:45
Kepler solved his eponymous equation
M = E − e sin(E)
by finding a fixed point of
E = M + e sin(E).
So guess a value of E and stick it into the right hand side. Then plug that value into the right hand side again. Kepler said a couple iterations should be enough. And a couple iterations are enough if the eccentricity e is small and you don’t need much accuracy.
The rate of convergence is determined by e. Kepler implicitly had in mind small values of e because he wasn’t aware of anything orbiting the sun in a highly elliptical orbit. Here’s an example with eccentricity 0.05, about the eccentricity of the orbits of Jupiter and Saturn.
from math import sin
M, E, e = 1, 1, 0.05
for _ in range(5):
E = M + e*sin(E)
residual = M - (E - e*sin(E))
print(residual)
The residual after just two iterations is 2.77 × 10−5. If you change e to 0.2, the eccentricity of Mercury’s orbit, it takes three iterations to get comparable accuracy. Mercury has the most eccentric orbit of any object Kepler would have known about.
Now suppose you’d like to solve for E when M = 1 for Halley’s comet, and you’d like an error of less than 10−8. Now you need 16 iterations.
C. F. W. Peters discovered a faster algorithm in 1891.
E1 = M + e sin(E0)
E2 = M + e sin(E1)
E3 = (E2E0 − E1²)/(E2 − 2E1 + E0)
Let’s look at the results of doing three iterations of Peters’ method for Halley’s comet.
M, E0, e = 1, 1, 0.967
for _ in range(3):
E1 = M + e*sin(E0)
E2 = M + e*sin(E1)
E3 = (E2*E0 - E1**2)/(E2 - 2*E1 + E0)
residual = M - (E - e*sin(E3))
print(residual)
E0 = E3 # for next iteration
This gives a residual of −7.23 × 10−10. Each iteration of Peters’ method requires a little more than twice as much work as an iteration of Kepler’s method, but 3 iterations of Peters’ method accomplished more than 16 iterations of Kepler’s method.
Peters’ algorithm from 1891 was a special case of Alexander Aitken’s series acceleration method published in 1926.
2026-06-08 03:06:36
An earlier post discussed how to solve Kepler’s equation
M = E − e sin(E)
using a sine series. You could also solve Kepler’s equation using a power series, which Lagrange did in 1771. Both approaches express E as a function of e and M, but from different perspectives. Bessel thought of his solution as a sum of sines in M, with coefficients that depend on e. Lagrange thought of his solution as a power series in e whose coefficients involve sines in M. You can rearrange the terms of either solution into the other.
The most interesting thing about the power series solution, in my opinion, is that it only converges for e less than roughly 2/3 while the sine series solution is valid for all e < 1. In astronomical terms, this means the power series solution works for the orbit of some planets but not others!
In our solar system, the planets all have eccentricity well below 2/3, but not all minor planets do. For example, the orbit of Eris has eccentricity 0.4407 but the orbit of Sedna has eccentricity 0.8549. And in other solar systems there are planets with eccentricity much greater than 2/3.
The radius of convergence for Lagrange’s power series solution is called the Laplace limit. Its value is eL = 0.6627…. There’s no obvious reason why there’s anything special about this value. There’s no astronomical reason for this value. It’s an artifact of the power series form of the solution.
If the series works for e = 0.66, you would reasonably think it works for e = 0.67, but that’s not the case. And if you’re observant, you might notice that although the series works for e = 0.66, it takes longer to converge than for smaller values of e; the rate of convergence is slowing down, warning you of danger ahead.
The exact value of eL is the unique real solution to the equation
There’s no obvious reason for this either. It has to do with finding the largest circle that can fit in a lens-shaped region of convergence. More on that here.
We can calculate eL with the following Python code.
from math import exp
from scipy.optimize import root_scalar
def f(x):
t = (1 + x*x)**0.5
return x*math.exp(t) - 1 - t
sol = root_scalar(f, bracket=[0, 1], method='brentq')
print(sol.root)
This prints 0.6627434193491817.
We can use the Lagrange inversion formula to find the series, just as Lagrange did two and a half centuries ago.
The powers of sine can be expanded into the sum of sines of various frequencies and differentiated, leading to the equation
The post The Laplace limit first appeared on John D. Cook.
2026-06-08 00:34:26
I ran across a cranky formula for π based on physical constants here
![\pi = \left( \frac{E}{\frac{1}{2}mc^{2}} \right)^{1/2} \left[ J_{\lambda} \cdot \lambda^{5} \left( e^{\frac{hc}{\lambda kT}} - 1 \right) \right]](https://www.johndcook.com/crank1.svg)
and decided to play around with it.
The source describes λ as “wavelength (chosen in the microwave region)” and I thought perhaps you could chose a value of λ to make the equation work. But as a comment pointed out, the bracketed expression is simply 2hc², independent of λ, due to Planck’s blackbody law. That means we can simplify the expression above to
Now the values of h and c are known. In fact, they’re now exactly known by definition: other SI units are defined in terms of h and c. The mass of an electron is known to 11 significant figures.
But E in the equation above is “Total energy of the universe.” I don’t even know what that means. Does it refer to the observable universe? Does it include dark energy? Does it include the energy equivalent of mass?
I asked a couple LLMs that the total energy of the universe might mean and what its value might be, and they said something like “Depends. It might be zero. It might be infinite. But if I had to say, I’d say around 1070 Joules.”
If we solve the equation above for E we get 2.8480347886530404 × 1019 Joules. I have no idea how to justify that.
The expression for π is not dimensionless. I suppose you could choose some nonstandard units that would make the equation work.
The source I linked to above cites Mathematical Cranks by Underwood Dudley, but I couldn’t find it in the book.
2026-06-07 02:46:02
The previous post very briefly said that the integral representation for Bessel functions was motived by solving Kepler’s equation. This post will go into more detail.
There are multiple ways to describe the position of a planet in an elliptical orbit around a star. For historical reasons, these descriptions have arcane names such as mean anomaly, true anomaly, and eccentric anomaly. This post explains how these three are related.
For this post, it is enough to say that often you know mean anomaly M and want to know eccentric anomaly E. These are related via Kepler’s equation
where e is the eccentricity of the orbit. You’d like to solve for E as a function of M and e, but there’s no elementary way to do that.
One way to solve Kepler’s equation is to take a guess at E and plug it into the right hand side of
to get a new E, and keep iterating until the two sides are closer together. I write more about this here.
Another approach to solving Kepler’s equation is to use Newton’s method. I write more about that here.
Still another approach is to expand E in a sine series and find the series coefficients. An advantage to this approach is that once you have the coefficients, you have an expression for E as a function of M, and you can plug in more values of M without having to solve Kepler’s equation for each value of M separately.
Kepler’s equation is easy to solve at E = 0 and at E = π. In both cases, E = M. So the function E − M is zero at both ends of [0, π], which suggests we try to expand E − M in a sine series
We then calculate the Fourier coefficients an as usual.
The second line uses integration by parts. The third line uses Kepler’s equation. The last line uses the definition of the Bessel functions Jn given in the previous post.
The post From Kepler to Bessel first appeared on John D. Cook.2026-06-05 20:04:28
Yesterday I wrote a post showing that the trapezoid rule evaluates the integral
very efficiently. But how do we know what the exact integral is for comparison? If you ask Mathematica, it will tell you the integral equals −2π J1(1) where J1 is a Bessel function. This may seem like rabbit out of a hat, but it’s actually a simple calculation given the integral definition of Bessel functions:
Since cosine is even, we can write our integral over [−π, π] as twice the integral over [0, π]. Then a change of variables turns this into the definition of Jn(z) with n = 1 and z = 1.
A deeper question is what have we accomplished by just giving a new name to essentially the same problem we started with. Another question is why in the world are Bessel functions defined as above.
As for what we’ve accomplished, we’ve related out integration problem to a very well-studied function. Bessel functions have been studied for two centuries and it’s easy to find software to evaluate them. Even the usually minimalist command line calculator bc has a function j(x, n) for evaluating Jn(x) for integer values of n. We could calculate our integral to 50 decimal places as follows.
~$ bc -l >>> scale = 50 >>> -8*a(1)*j(1,1) -2.76491937476833705153256665538788207487495025542883
Note that bc doesn’t have a value of π built in, but a(x) evaluates the arctangent function, and π = 4 arctan(1).
There are multiple ways of defining Bessel functions. The three main ways would be in terms of their power series, in terms of the differential equations they solve, and in terms of their integral representation. Friedrich Bessel defined what we now call Bessel functions of the first kind, the Jn functions, in terms of their integral representations.
Why did Bessel care about these integrals? They came out of his calculations in celestial mechanics. One example of this is solving Kepler’s equation with Fourier series; the Fourier coefficients are given by Bessel functions. This is worked out in detail in the next post.
Bessel functions had occurred in applications before Mr. Bessel drew attention to them. The functions are named after him because he was the first to systematically study them.
Mathematics is developed inductively but taught deductively. It’s common for things to be kicked around for years before someone decides they deserve a name and systematic study. See this post on the central limit theorem for another example. The CLT is older than the Gaussian distribution, even older than Gauss.
2026-06-05 04:12:11
One reason people study Latin is that it is the ancestor of many modern languages. English derives from West Germanic languages, not from Latin, but much of English vocabulary, perhaps as much as 60%, derives from Latin, either directly or indirectly through French.
Knowing a bit of Latin makes sense of many things that would otherwise seem completely arbitrary, such as why the symbols for gold, silver, and lead are Au, Ag, and Pb respectively.
Similarly, ed(1) is the Latin of Linux [1]. Many conventions in command line utilities follow conventions that go back to the ed(1) line editor. They may go back even further. Just as Latin didn’t come out of nowhere, neither did ed(1), but you can’t go back indefinitely. It’s convenient to start history somewhere, and this post will start with ed(1) just as much discussion of Western linguistics starts with Latin.
The following are features of ed(1) that live on in sed, awk, grep, vi, perl, bash, etc.
[1] Because the name “ed” is so short, and looks so much like the name Ed, it’s convenient to use its full Unix name ed(1). The parenthesized number is used to disambiguate different things that have the same name, such as the user command kill(1) and the system call kill(2). There is no ed(2) or any other higher-numbered ed. The number is there to make the name stand out, not to disambiguate anything.
The post The Latin of Linux first appeared on John D. Cook.
