Formally proving a calculation with Claude and Lean

(832) 422-8646

Contact

I ran an experiment today to see whether Claude [1] 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 sorrys?

Real.besselJ definition or a Kepler namespace. The three remaining sorrys 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.

/-<br>KeplerBessel.lean<br>Formal proof that the n-th Fourier sine-series coefficient of (E − M),<br>where E and M are related by Kepler's equation M = E − e · sin E,<br>equals aₙ = (2/n) · Jₙ(n·e),<br>with Jₙ the Bessel function of the first kind of integer order n.

Mathematical content<br>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<br>= (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<br>= (2/n) · Jₙ(n·e).

The last step uses the Bessel integral representation<br>Jₙ(x) = (1/π) ∫₀^π cos(n·θ − x·sin θ) dθ.<br>-/

import Mathlib

open Real MeasureTheory intervalIntegral Filter Set

noncomputable section

/-! ---------------------------------------------------------------<br>§1 Variables<br>--------------------------------------------------------------- -/

variable (e : ℝ) (he : 0 ≤ e) (he1 : e x - e * sin x, so HasDerivAt follows directly<br>-- from sub-rule and const_mul applied to hasDerivAt_sin.<br>(hasDerivAt_id E).sub ((hasDerivAt_sin E).const_mul e)

/-- The derivative of `keplerMap e` is positive when e keplerMap_hasDerivAt e' E)<br>(fun E => keplerMap_deriv_pos he' he1' E)

/-!<br>We axiomatise the inverse eccAnom : ℝ → ℝ → ℝ and its key<br>properties, all of which follow from the Inverse Function Theorem<br>applied to the smooth, strictly monotone map keplerMap e.<br>-/

/-- The eccentric anomaly: the smooth right-inverse of `keplerMap e`. -/<br>axiom eccAnom (e : ℝ) : ℝ → ℝ

/-- `eccAnom e M` satisfies Kepler's equation. -/<br>axiom eccAnom_kepler (e M : ℝ) :<br>keplerMap e (eccAnom e M) = M

/-- `eccAnom e` is differentiable, derivative = 1/(1 − e·cos(eccAnom e M)). -/<br>axiom eccAnom_hasDerivAt (e M : ℝ) :<br>HasDerivAt (eccAnom e) (1 / (1 - e * cos (eccAnom e M))) M

/-- Boundary value at 0. -/<br>axiom eccAnom_zero (e : ℝ) : eccAnom e 0 = 0

/-- Boundary value at π. -/<br>axiom eccAnom_pi (e : ℝ) : eccAnom e π = π

/-! ---------------------------------------------------------------<br>§3 Bessel function of the first kind (integer order)

Defined by the classical integral representation.<br>--------------------------------------------------------------- -/

/-- Bessel function J_n(x) via its integral representation. -/<br>def besselJ (n : ℕ) (x : ℝ) : ℝ :=<br>(1 / π) * ∫ θ in (0 : ℝ)..π, cos (↑n * θ - x * sin θ)

/-! ---------------------------------------------------------------<br>§4 Fourier coefficient

Named keplerFourierCoeff to avoid clashing with Mathlib's own<br>`fourierCoeff` which is defined on AddCircle.<br>--------------------------------------------------------------- -/

/-- The n-th Fourier sine coefficient of eccAnom e M − M on [0,π]. -/<br>def keplerFourierCoeff (e : ℝ) (n : ℕ) : ℝ :=<br>(2 / π) * ∫ M in (0 : ℝ)..π,<br>(eccAnom e M - M) * sin (↑n * M)

/-! ---------------------------------------------------------------<br>§5 Main theorem<br>--------------------------------------------------------------- -/

/--<br>**Main theorem.** For n ≥ 1, the Fourier sine coefficient of the<br>eccentric-anomaly displacement satisfies aₙ = (2/n) · Jₙ(n·e).<br>-/<br>theorem keplerFourierCoeff_eq_besselJ (n : ℕ) (hn : 1 ≤ n) :<br>keplerFourierCoeff e n = (2 / (n : ℝ)) * besselJ n (↑n * e) := by<br>simp only [keplerFourierCoeff, besselJ]<br>-- Goal:<br>-- (2/π) · ∫₀^π (E(M)−M)·sin(nM) dM<br>-- = (2/n) · (1/π) · ∫₀^π cos(nθ − ne·sinθ) dθ

-- ── Step 1: Integration by parts ─────────────────────────────────────<br>-- u = E(M)−M, dv = sin(nM)dM → v = −cos(nM)/n<br>-- Boundary: [uv]₀^π = 0 by eccAnom_zero, eccAnom_pi.<br>-- Result:...