2607.19283v1 / Computations.lean
all files
import Mathlib
set_option maxHeartbeats 0
set_option maxRecDepth 100000
/-!
# Finite checks from the ENO--TV paper
This file independently checks the finite calculations in Sections 6.3--6.5
and Appendix A of Li--Wu, arXiv:2607.19283v1. All exhaustive certificates use
kernel reduction and theorem-producing tactics, not native-code evaluation.
-/
open scoped BigOperators
namespace ENOTVComputations
/-! ## The Gaussian-binomial coefficient in (6.20) -/
/-- Partitions of `10` fitting in a `3 × 7` rectangle, represented by three rows. -/
def partitions3x7 : List (ℕ × ℕ × ℕ) :=
(List.range 8).flatMap fun a =>
(List.range 8).flatMap fun b =>
(List.range 8).map (fun c => (a, b, c)) |>.filter fun (a, b, c) =>
b ≤ a && c ≤ b && a + b + c = 10
/-- The coefficient `[z^10] {10 \choose 3}_z` is `10`, as asserted in (6.20). -/
theorem cohomology_dimension_2_7 : partitions3x7.length = 10 := by
decide
/-! ## Fourth-order moment coefficients from (6.26) and (6.40) -/
/-- The closed formula (6.26), specialized to natural moment indices. -/
def alpha (mu r : ℕ) : ℚ :=
(-1 : ℚ) ^ (r + 1) * (2 * (Nat.factorial mu : ℚ) ^ 2) /
(r * Nat.factorial (mu - r) * Nat.factorial (mu + r))
/-- At moment index two, `α₁ = 4/3` and `α₂ = -1/6`, as in (6.66). -/
theorem fourth_order_moment_coefficients :
alpha 2 1 = 4 / 3 ∧ alpha 2 2 = -1 / 6 := by
norm_num [alpha]
/-- Formula (6.40) then gives the sole coefficient `β₁ = -1/6`. -/
theorem fourth_order_beta : alpha 2 2 * (2 - 1) = -1 / 6 := by
norm_num [alpha]
/-! ## Direct coefficient evaluation for the degree-seven obstruction -/
/-- Exponents of `u`, `c0`, `c1`, and `c2`, respectively. -/
structure Exp4 where
eu : ℕ
e0 : ℕ
e1 : ℕ
e2 : ℕ
deriving DecidableEq
def Exp4.total (e : Exp4) : ℕ := e.eu + e.e0 + e.e1 + e.e2
/-- A rational linear form in `u,c0,c1,c2`. -/
structure Lin4 where
lu : ℚ
l0 : ℚ
l1 : ℚ
l2 : ℚ
/-- Multinomial coefficient formula for one coefficient of a power of a linear form. -/
def linearPowerCoeff (l : Lin4) (n : ℕ) (e : Exp4) : ℚ :=
if e.total = n then
(Nat.factorial n : ℚ) /
((Nat.factorial e.eu : ℚ) * Nat.factorial e.e0 * Nat.factorial e.e1 * Nat.factorial e.e2) *
l.lu ^ e.eu * l.l0 ^ e.e0 * l.l1 ^ e.e1 * l.l2 ^ e.e2
else 0
/-- Cauchy product of two coefficient functions, evaluated in the finite box below `e`. -/
def coeffMul (p q : Exp4 → ℚ) (e : Exp4) : ℚ :=
∑ au ∈ Finset.range (e.eu + 1),
∑ a0 ∈ Finset.range (e.e0 + 1),
∑ a1 ∈ Finset.range (e.e1 + 1),
∑ a2 ∈ Finset.range (e.e2 + 1),
p ⟨au, a0, a1, a2⟩ * q ⟨e.eu - au, e.e0 - a0, e.e1 - a1, e.e2 - a2⟩
def linU : Lin4 := ⟨1, 0, 0, 0⟩
def linTu : Lin4 := ⟨1, 1, 0, 0⟩
def linT2u : Lin4 := ⟨1, 2, 1, 0⟩
def linTinvU : Lin4 := ⟨1, -1, 1, -1⟩
/-- One coefficient of the cubic two-point entropy-conservative flux (6.65). -/
def fStarCoeff (x y : Lin4) (e : Exp4) : ℚ :=
(1 / 12) *
(linearPowerCoeff x 3 e
+ coeffMul (linearPowerCoeff x 2) (linearPowerCoeff y 1) e
+ coeffMul (linearPowerCoeff x 1) (linearPowerCoeff y 2) e
+ linearPowerCoeff y 3 e)
/-- One coefficient of the moment-index-two fourth-order FMT flux (6.67). -/
def fMomentTwoCoeff (e : Exp4) : ℚ :=
(4 / 3) * fStarCoeff linU linTu e
- (1 / 6) * (fStarCoeff linU linT2u e + fStarCoeff linTinvU linTu e)
/-- One coefficient of the raw total-degree-seven mismatch (6.68)/(A.7). -/
def gRawCoeff (e : Exp4) : ℚ :=
coeffMul
(fun q => linearPowerCoeff linTu 4 q - linearPowerCoeff linU 4 q)
fMomentTwoCoeff e
- (4 / 21) * (linearPowerCoeff linTu 7 e - linearPowerCoeff linU 7 e)
def coeff4 (a b c d : ℕ) : ℚ := gRawCoeff ⟨a, b, c, d⟩
/-- All seven nonzero raw coefficients displayed in (A.11) are reproduced. -/
theorem raw_coefficients :
coeff4 5 1 1 0 = -2 / 3 ∧
coeff4 5 1 0 1 = 1 / 3 ∧
coeff4 4 2 1 0 = -5 / 3 ∧
coeff4 4 2 0 1 = 5 / 18 ∧
coeff4 4 1 2 0 = -4 / 9 ∧
coeff4 3 3 1 0 = -22 / 9 ∧
coeff4 2 5 0 0 = -17 / 6 := by
decide +kernel
/-- The other five coefficients in the support of `L` vanish, as claimed. -/
theorem remaining_functional_coefficients_vanish :
coeff4 6 1 0 0 = 0 ∧
coeff4 6 0 0 1 = 0 ∧
coeff4 5 0 1 1 = 0 ∧
coeff4 6 0 1 0 = 0 ∧
coeff4 5 0 2 0 = 0 := by
decide +kernel
/-- The rational functional listed in Appendix A, evaluated on `G_raw`. -/
def obstructionPairing : ℚ :=
(-180) * coeff4 6 1 0 0
+ 1620 * coeff4 6 0 0 1
- 90 * coeff4 5 1 0 1
- 180 * coeff4 5 0 1 1
+ 36 * coeff4 4 2 0 1
- 45 * coeff4 3 3 1 0
+ 720 * coeff4 6 0 1 0
- 90 * coeff4 5 1 1 0
+ 240 * coeff4 5 0 2 0
- 12 * coeff4 4 2 1 0
+ 72 * coeff4 4 1 2 0
+ 60 * coeff4 2 5 0 0
/-- The decisive pairing in (A.12) is exactly `-32`. -/
theorem obstruction_pairing : obstructionPairing = -32 := by
rcases raw_coefficients with ⟨h₁, h₂, h₃, h₄, h₅, h₆, h₇⟩
rcases remaining_functional_coefficients_vanish with ⟨h₈, h₉, h₁₀, h₁₁, h₁₂⟩
norm_num [obstructionPairing, h₁, h₂, h₃, h₄, h₅, h₆, h₇,
h₈, h₉, h₁₀, h₁₁, h₁₂]
/-! ## Exhaustive verification that the functional kills the nilpotent range -/
/-- The table `λ_{a,b,c,d}` from Appendix A, extended by zero. -/
def lambda (a b c d : ℕ) : ℚ :=
if (a, b, c, d) = (6, 1, 0, 0) then -180 else
if (a, b, c, d) = (6, 0, 0, 1) then 1620 else
if (a, b, c, d) = (5, 1, 0, 1) then -90 else
if (a, b, c, d) = (5, 0, 1, 1) then -180 else
if (a, b, c, d) = (4, 2, 0, 1) then 36 else
if (a, b, c, d) = (3, 3, 1, 0) then -45 else
if (a, b, c, d) = (6, 0, 1, 0) then 720 else
if (a, b, c, d) = (5, 1, 1, 0) then -90 else
if (a, b, c, d) = (5, 0, 2, 0) then 240 else
if (a, b, c, d) = (4, 2, 1, 0) then -12 else
if (a, b, c, d) = (4, 1, 2, 0) then 72 else
if (a, b, c, d) = (2, 5, 0, 0) then 60 else 0
/-- Formula (A.9) for `L (N₇ (u^a c0^b c1^c c2^d))`. -/
def nilpotentPairing (a b c d : ℕ) : ℚ :=
a * (lambda (a - 1) (b + 1) c d
- (1 / 2) * lambda (a - 1) b (c + 1) d
+ (1 / 3) * lambda (a - 1) b c (d + 1))
+ b * (lambda a (b - 1) (c + 1) d
- (1 / 2) * lambda a (b - 1) c (d + 1))
+ c * lambda a b (c - 1) (d + 1)
/-- The complete list of exponent quadruples of monomials of total degree seven. -/
def degreeSevenBasis : List Exp4 :=
(List.range 8).flatMap fun a =>
(List.range 8).flatMap fun b =>
(List.range 8).flatMap fun c =>
(List.range 8).map (fun d => ⟨a, b, c, d⟩) |>.filter fun e => e.total = 7
/-- There are `binom (7+4-1) (4-1) = 120` degree-seven basis monomials. -/
theorem degreeSevenBasis_card : degreeSevenBasis.length = 120 := by
decide
/-- Boolean certificate that (A.9) vanishes on every degree-seven monomial. -/
def nilpotentRangeCheck : Bool :=
degreeSevenBasis.all fun e => nilpotentPairing e.eu e.e0 e.e1 e.e2 == 0
/--
Exhaustive check of all 120 degree-seven basis monomials: `L ∘ N_7 = 0`,
the finite cokernel calculation behind (A.10).
-/
theorem functional_annihilates_nilpotent_range : nilpotentRangeCheck = true := by
decide +kernel
/-! ## A concrete fourth-order Euler-block instance -/
def listDiff : List ℤ → List ℤ
| x :: y :: xs => (y - x) :: listDiff (y :: xs)
| _ => []
def iterListDiff : ℕ → List ℤ → List ℤ
| 0, xs => xs
| n + 1, xs => iterListDiff n (listDiff xs)
def eulerP2 (L t : ℕ) : ℤ := (t : ℤ) ^ 2 - 2 * (L : ℤ) * t
def eulerBlocks4 (L N : ℕ) : List ℤ :=
(List.range N).flatMap fun n =>
(List.range (2 * L)).map fun t => if Even n then eulerP2 L t else -eulerP2 L t
def listMaxNatAbs (xs : List ℤ) : ℕ := xs.foldl (fun m x => max m x.natAbs) 0
def jumpMoment5 (xs : List ℤ) : ℕ := xs.foldl (fun acc x => acc + x.natAbs ^ 5) 0
def cumulativeAmplitude (xs : List ℤ) : ℕ :=
listMaxNatAbs (xs.scanl (· + ·) 0)
def localizedE4 (xs : List ℤ) : ℕ :=
let padded := [0, 0, 0] ++ xs ++ [0, 0, 0]
let d3 := iterListDiff 3 padded
(List.range d3.length).foldl (fun acc j =>
acc + listMaxNatAbs ((padded.drop j).take 4) * (d3.getD j 0).natAbs) 0
def eulerBlockSummary (L N : ℕ) : ℕ × ℕ × ℕ :=
let xs := eulerBlocks4 L N
(jumpMoment5 xs, cumulativeAmplitude xs, localizedE4 xs)
/-- Two exact finite instances of the fourth-order Euler-block construction. -/
theorem eulerBlock_instances :
eulerBlockSummary 2 4 = (6040, 10, 106) ∧
eulerBlockSummary 4 16 = (49577664, 84, 1034) ∧
(eulerBlocks4 2 4).sum = 0 ∧ (eulerBlocks4 4 16).sum = 0 := by
decide +kernel
/-- The normalized moment/localized-source quotient already increases in these two instances. -/
theorem eulerBlock_normalized_ratio_increases :
6040 * 84 ^ 3 * 1034 < 49577664 * 10 ^ 3 * 106 := by
norm_num
/-! ## The strictly convex entropy's numerical bound -/
/-- The second derivative from (6.63). -/
def etaStarSecond (x : ℝ) : ℝ := 1 + 4 * x ^ 3 + 7 * x ^ 6
/-- The paper's global strict-convexity bound `η⋆'' ≥ 3/7`. -/
theorem etaStarSecond_lower_bound (x : ℝ) : (3 : ℝ) / 7 ≤ etaStarSecond x := by
dsimp [etaStarSecond]
nlinarith [sq_nonneg (x ^ 3 + (2 : ℝ) / 7)]
#print axioms cohomology_dimension_2_7
#print axioms fourth_order_moment_coefficients
#print axioms raw_coefficients
#print axioms obstruction_pairing
#print axioms functional_annihilates_nilpotent_range
#print axioms eulerBlock_instances
#print axioms etaStarSecond_lower_bound
end ENOTVComputations