Formally proving a calculation with Claude and Lean

Anthropic's Claude AI generated a Lean formal proof for a Fourier coefficient calculation involving Bessel functions, requiring eight iterations to fix errors and produce a working proof with four 'sorry' placeholders for standard analysis lemmas.

I ran an experiment today to see whether Claude 1 could generate Lean code to prove a calculation at the bottom of this post https://www.johndcook.com/blog/2026/06/06/from-kepler-to-bessel/ , 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 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. --------------------------------------------------------------- -/ 1 I used the default setting: Sonnet 4.6, low effort, “thinking” turned off.