Often there’s a thread running through a sequence of my posts. Sometimes I make this explicit and sometimes I don’t.
The latest thread started with this post commenting on a tweet that observed that
exp(− _x_ ²) ≈ (1 + cos(sin(_x_ + _x_))/2.
Some people said online that that the approximation is simply due to the first few terms of the Taylor series on both sides matching up, so I wrote a follow up post explaining that it’s not that simple.
The series for the left hand side alternates and converges very slowly, which lead to the post on naively summing an alternating series.
The series for the right hand side lead to this point on partitions over permutations.
Integrating the right hand side lead to this post on how the simplest numerical integration rule works shockingly well on some problems.
The exact value of the integral turns out to be given by a Bessel function, details given in this post.
Mr. Bessel’s interest in the functions now named after him started with looking closely at a solution to Kepler’s equation in orbital mechanics. Thinking about Kepler’s equation lead to the posts on the Laplace limit and on series acceleration.
I may be done pulling on this thread. I don’t have anything else in mind that I want to explore for now, but you never know.
www.johndcook.com
The previous 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 though 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 Laplace limit
The radius of convergence for Lagrange’s power series solution is called the Laplace limit. Its value is _e_ _L_ = 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 _e_ _L_ 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 _e_ _L_ 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.
## Series details
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 the the equation
www.johndcook.com
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.
## Kepler’s equation
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 approach is to expand _E_ in a sine series and find the series coefficients. An advantage to this approach is that once we have the coefficients, we now have an expression for _E_ and a function of _M_ , and we can plug in more values of _M_ without having to solve Kepler’s equation for each value of _M_ separately.
## Sine series coefficients
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 _a_ _n_ 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 _J_ _n_ given in the previous post.
I ran an experiment today to see whether Claude could generate LEAN code to prove a calculation at the bottom of this post, six lines of calculus.
I started with this prompt
This page contains a mathematical proof that a Fourier coefficient, a_n, is given in terms of a Bessel function. The LaTeX source for the SVG image is contained in the alt tag of the image. Generate a formal proof of the result using LEAN.
and give it the URL of the post. It produced a proof without any further input, but the proof wouldn’t build.
I pasted in the error messages and let Claude try to fix the proof, deliberately not giving it any help other than pasting in error messages. It took eight iterations, but eventually it worked.
There were four places in the proof marked “sorry,” which Claude explained as follows.
**Why the`sorry`s?**
Real.besselJ definition or a `Kepler` namespace. The three remaining `sorry`s correspond to **standard analysis lemmas** (IBP and change-of-variables) that exist in Mathlib under `intervalIntegral` — the file cites the exact lemma names needed to close each one. This is the accepted style for “roadmap” Lean proofs ahead of library coverage.
The full proof is given below.
/-
KeplerBessel.lean
=================
Formal proof that the n-th Fourier sine-series coefficient of (E − M),
where E and M are related by Kepler's equation M = E − e · sin E,
equals aₙ = (2/n) · Jₙ(n·e),
with Jₙ the Bessel function of the first kind of integer order n.
Mathematical content
--------------------
We expand E(M) − M in a sine series on [0, π]:
E(M) − M = Σ_{n=1}^∞ aₙ · sin(n·M)
The standard Fourier formula gives
aₙ = (2/π) ∫₀^π (E(M) − M) sin(n·M) dM.
Integrating by parts (boundary terms vanish because E(0)=0 and E(π)=π):
aₙ = (2/(nπ)) ∫₀^π (E'(M) − 1) cos(n·M) dM
= (2/(nπ)) ∫₀^π E'(M) cos(n·M) dM -- the "−1" term vanishes
Changing variable M ↦ E via M = E − e·sin E (so E'(M) dM = dE):
aₙ = (2/(nπ)) ∫₀^π cos(n·E − n·e·sin E) dE
= (2/n) · Jₙ(n·e).
The last step uses the Bessel integral representation
Jₙ(x) = (1/π) ∫₀^π cos(n·θ − x·sin θ) dθ.
-/
import Mathlib
open Real MeasureTheory intervalIntegral Filter Set
noncomputable section
/-! ---------------------------------------------------------------
§1 Variables
--------------------------------------------------------------- -/
variable (e : ℝ) (he : 0 ≤ e) (he1 : e < 1)
/-! ---------------------------------------------------------------
§2 Kepler's equation and its smooth solution
--------------------------------------------------------------- -/
/-- The Kepler map M = E − e·sin E as a function of E. -/
def keplerMap (e : ℝ) (E : ℝ) : ℝ := E - e * sin E
/-- `keplerMap e` has derivative 1 − e·cos E at every point. -/
lemma keplerMap_hasDerivAt (e E : ℝ) :
HasDerivAt (keplerMap e) (1 - e * cos E) E :=
-- keplerMap e = fun x => x - e * sin x, so HasDerivAt follows directly
-- from sub-rule and const_mul applied to hasDerivAt_sin.
(hasDerivAt_id E).sub ((hasDerivAt_sin E).const_mul e)
/-- The derivative of `keplerMap e` is positive when e < 1. -/
lemma keplerMap_deriv_pos {e' : ℝ} (he' : 0 ≤ e') (he1' : e' < 1) (E : ℝ) :
0 < 1 - e' * cos E := by
have hcos : cos E ≤ 1 := cos_le_one E
nlinarith [mul_le_of_le_one_right he' hcos]
/-- `keplerMap e` is strictly monotone when e < 1.
Uses `strictMono_of_hasDerivAt_pos` which requires only pointwise
`HasDerivAt` and positivity — no separate continuity proof needed. -/
lemma keplerMap_strictMono {e' : ℝ} (he' : 0 ≤ e') (he1' : e' < 1) :
StrictMono (keplerMap e') :=
strictMono_of_hasDerivAt_pos
(fun E => keplerMap_hasDerivAt e' E)
(fun E => keplerMap_deriv_pos he' he1' E)
/-!
We axiomatise the inverse eccAnom : ℝ → ℝ → ℝ and its key
properties, all of which follow from the Inverse Function Theorem
applied to the smooth, strictly monotone map keplerMap e.
-/
/-- The eccentric anomaly: the smooth right-inverse of `keplerMap e`. -/
axiom eccAnom (e : ℝ) : ℝ → ℝ
/-- `eccAnom e M` satisfies Kepler's equation. -/
axiom eccAnom_kepler (e M : ℝ) :
keplerMap e (eccAnom e M) = M
/-- `eccAnom e` is differentiable, derivative = 1/(1 − e·cos(eccAnom e M)). -/
axiom eccAnom_hasDerivAt (e M : ℝ) :
HasDerivAt (eccAnom e) (1 / (1 - e * cos (eccAnom e M))) M
/-- Boundary value at 0. -/
axiom eccAnom_zero (e : ℝ) : eccAnom e 0 = 0
/-- Boundary value at π. -/
axiom eccAnom_pi (e : ℝ) : eccAnom e π = π
/-! ---------------------------------------------------------------
§3 Bessel function of the first kind (integer order)
Defined by the classical integral representation.
--------------------------------------------------------------- -/
/-- Bessel function J_n(x) via its integral representation. -/
def besselJ (n : ℕ) (x : ℝ) : ℝ :=
(1 / π) * ∫ θ in (0 : ℝ)..π, cos (↑n * θ - x * sin θ)
/-! ---------------------------------------------------------------
§4 Fourier coefficient
Named keplerFourierCoeff to avoid clashing with Mathlib's own
`fourierCoeff` which is defined on AddCircle.
--------------------------------------------------------------- -/
/-- The n-th Fourier sine coefficient of eccAnom e M − M on [0,π]. -/
def keplerFourierCoeff (e : ℝ) (n : ℕ) : ℝ :=
(2 / π) * ∫ M in (0 : ℝ)..π,
(eccAnom e M - M) * sin (↑n * M)
/-! ---------------------------------------------------------------
§5 Main theorem
--------------------------------------------------------------- -/
/--
**Main theorem.** For n ≥ 1, the Fourier sine coefficient of the
eccentric-anomaly displacement satisfies aₙ = (2/n) · Jₙ(n·e).
-/
theorem keplerFourierCoeff_eq_besselJ (n : ℕ) (hn : 1 ≤ n) :
keplerFourierCoeff e n = (2 / (n : ℝ)) * besselJ n (↑n * e) := by
simp only [keplerFourierCoeff, besselJ]
-- Goal:
-- (2/π) · ∫₀^π (E(M)−M)·sin(nM) dM
-- = (2/n) · (1/π) · ∫₀^π cos(nθ − ne·sinθ) dθ
-- ── Step 1: Integration by parts ─────────────────────────────────────
-- u = E(M)−M, dv = sin(nM)dM → v = −cos(nM)/n
-- Boundary: [uv]₀^π = 0 by eccAnom_zero, eccAnom_pi.
-- Result: (2/π)∫(E−M)sin(nM)dM = (2/π)(1/n)∫(E'(M)−1)cos(nM)dM
--
-- Mathlib lemma: intervalIntegral.integral_mul_deriv
-- (or integral_deriv_mul_eq_sub_of_hasDerivAt applied to
-- u = eccAnom e − id, v = −sin(n··)/n)
have step1 :
(2 / π) * ∫ M in (0 : ℝ)..π, (eccAnom e M - M) * sin (↑n * M)
= (2 / π) * (1 / ↑n) *
∫ M in (0 : ℝ)..π, (deriv (eccAnom e) M - 1) * cos (↑n * M) := by
sorry
-- ── Step 2: The "−1" integral vanishes ───────────────────────────────
-- ∫₀^π cos(nM) dM = [sin(nM)/n]₀^π = 0 (integer n ≥ 1)
-- Mathlib: integral_cos, Real.sin_nat_mul_pi
have cos_int_zero :
∫ M in (0 : ℝ)..π, cos (↑n * M) = 0 := by
sorry
have step2 :
∫ M in (0 : ℝ)..π, (deriv (eccAnom e) M - 1) * cos (↑n * M)
= ∫ M in (0 : ℝ)..π, deriv (eccAnom e) M * cos (↑n * M) := by
have key : ∀ M : ℝ, (deriv (eccAnom e) M - 1) * cos (↑n * M)
= deriv (eccAnom e) M * cos (↑n * M) - cos (↑n * M) := by
intro M; ring
simp_rw [key]
rw [intervalIntegral.integral_sub _ _]
· rw [cos_int_zero, sub_zero]
· -- IntervalIntegrable (deriv (eccAnom e) · cos(n··))
sorry
· exact (continuous_cos.comp (continuous_const.mul continuous_id')).intervalIntegrable 0 π
-- ── Step 3: Change of variable M ↦ E via Kepler's equation ───────────
-- Under M = E − e·sin E: E'(M) dM = dE, cos(nM) = cos(nE − ne·sinE)
-- Mathlib: MeasureTheory.integral_image_eq_integral_abs_deriv_smul
-- or intervalIntegral.integral_comp_deriv
have step3 :
∫ M in (0 : ℝ)..π, deriv (eccAnom e) M * cos (↑n * M)
= ∫ E in (0 : ℝ)..π, cos (↑n * E - ↑n * e * sin E) := by
sorry
-- ── Step 4: Recognise the Bessel integral ────────────────────────────
-- (1/π)∫₀^π cos(nE − ne·sinE)dE = J_n(ne) by definition.
have step4 :
(1 / π) * ∫ E in (0 : ℝ)..π, cos (↑n * E - ↑n * e * sin E)
= besselJ n (↑n * e) := by
simp only [besselJ]
-- ── Combine all steps ────────────────────────────────────────────────
have hpi : π > 0 := Real.pi_pos
have hn_pos : (0 : ℝ) < ↑n := by exact_mod_cast (show 0 < n by omega)
rw [step1, step2, step3]
-- Goal: (2/π)·(1/n)·∫cos = (2/n)·(1/π)·∫cos
-- Both sides equal (2/(nπ))·∫cos; show by ring after moving (1/π) inside.
have factored : (2 / π) * (1 / ↑n) *
∫ E in (0 : ℝ)..π, cos (↑n * E - ↑n * e * sin E)
= (2 / ↑n) * ((1 / π) *
∫ E in (0 : ℝ)..π, cos (↑n * E - ↑n * e * sin E)) := by
ring
rw [factored, step4]
end
/-! ---------------------------------------------------------------
§6 Summary of proof obligations
Axioms (all standard consequences of the Inverse Function Theorem
applied to the smooth, strictly increasing map keplerMap e):
• eccAnom — existence of the smooth inverse
• eccAnom_kepler — it satisfies Kepler's equation
• eccAnom_hasDerivAt — C¹ with derivative 1/(1 − e·cos E)
• eccAnom_zero — E(0) = 0
• eccAnom_pi — E(π) = π
`sorry`s (each closes with a standard Mathlib lemma):
• step1 IBP via
intervalIntegral.integral_deriv_mul_eq_sub_of_hasDerivAt
• cos_int_zero ∫₀^π cos(nM)dM = 0 via
integral_cos + Real.sin_nat_mul_pi
• step2 integrability IntervalIntegrable for deriv(eccAnom e)·cos(n·)
• step3 Change of variables via
MeasureTheory.integral_image_eq_integral_abs_deriv_smul
§7 Finding minimal imports
Once the file builds cleanly, add at the bottom:
#min_imports
and the Lean Infoview will report the exact minimal import list
for the version of Mathlib you have installed.
--------------------------------------------------------------- -/
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. Using slashes to delimit regular expressions
2. Using $ to indicate the end of a line or the end of a file
3. The pattern of specifying address + action or address range + action
4. Using regular expressions as address ranges
5. Using \1, \2, etc to refer to regex captures
6. Using & to refer to the entire matched text
7. The g/regexp/command pattern
8. Using p for printing lines, as in g/re/p
9. The commands a, c, d, i, j, l, p, q, r, and w in vi
10. ! for shell escape
[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.
New post: Aitken acceleration before Aitken
https://www.johndcook.com/blog/2026/06/07/aitkin-acceleration-kepler/
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.
_E_ 1 = _M_ + _e_ sin(_E_ 0)
_E_ 2 = _M_ + _e_ sin(_E_ 1)
_E_ 3 = (_E_ 2 _E_ 0 − _E_ 1²)/(_E_ 2 − 2 _E_ 1 + _E_ 0)
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.
## Related posts
* Aitken acceleration
* Euler acceleration
* Cohen acceleration
