2607.19283v1 / ENOTV/Localization.lean
all files
import ENOTV.Energy
import Mathlib.LinearAlgebra.Lagrange
/-!
# ENO localization and the paper's local-amplitude lemma
This file separates the earlier Fjordholm--Mishra--Tadmor theorem used by
the paper from the elementary localization argument proved in the paper.
The former is exposed as a proposition that can be passed as a hypothesis;
the latter is proved below from the stencil recursion.
-/
noncomputable section
open scoped BigOperators
open Polynomial
namespace ENOTV
/-! ## A literal reconstruction, used to state the external FMT theorem -/
/-- Cumulative interface value at the boundary indexed by `n`. -/
def cumulative (u : CSeq) (n : ℤ) : ℝ :=
u.sum fun j x => if j < n then x else 0
/-- The degree-`k` Lagrange interpolant of cumulative interface values on
the `k+1` boundaries of the selected `k`-cell stencil. -/
def cumulativeInterpolant (k : ℕ) (u : CSeq) (i : ℤ) : ℝ[X] :=
let r := enoLeft (u : Seq) i (k - 1)
Lagrange.interpolate Finset.univ
(fun q : Fin (k + 1) => ((r + (q : ℕ) : ℤ) : ℝ))
(fun q : Fin (k + 1) => cumulative u (r + (q : ℕ)))
/-- Literal reconstructed trace: the derivative of the cumulative
interpolant evaluated at boundary `i+1`. -/
def literalTrace (k : ℕ) (u : CSeq) (i : ℤ) : ℝ :=
(derivative (cumulativeInterpolant k u i)).eval (i + 1)
/-- Literal right-minus-left interface jump of the two neighboring ENO
reconstruction polynomials. -/
def literalReconstructedJump (k : ℕ) (u : CSeq) (i : ℤ) : ℝ :=
literalTrace k u (i + 1) - literalTrace k u i
/-- The ENO source written naively from the two Lagrange-reconstructed
interface traces, before applying any localization theorem. -/
def literalSource (k : ℕ) (u : CSeq) : ℝ :=
(cdiff u).sum fun i ai => ai * literalReconstructedJump k u i
/-- The exact external theorem about classical ENO used by the paper:
the literal cell-average reconstruction equals the localized FMT formula,
and every term in that formula has the sign of the underlying jump.
Results which use this proposition take a value of it as an explicitly
declared hypothesis.
-/
def FMTReconstructionTheorem : Prop :=
(∀ (k : ℕ), 2 ≤ k → ∀ (u : CSeq) (i : ℤ),
literalReconstructedJump k u i =
reconstructedJump k (u : Seq) i) ∧
(∀ (k : ℕ), 2 ≤ k → ∀ (u : CSeq) (i j : ℤ),
selectedInterval k (u : Seq) i j →
0 ≤ cdiff u i *
(gamma k (Int.toNat (i - j)) *
cdiffIter (k - 1) (cdiff u) j))
/-- The naive reconstruction source and the localized source coincide.
This is the explicit semantic bridge from the definition above to the
quantity used throughout the estimates. -/
theorem literalSource_eq_cSource
(hFMT : FMTReconstructionTheorem) {k : ℕ} (hk : 2 ≤ k) (u : CSeq) :
literalSource k u = cSource k u := by
unfold literalSource cSource
apply Finsupp.sum_congr
intro i ai
rw [hFMT.1 k hk u i]
/-! ## Difference estimates on a finite block -/
theorem diffIter_eq_fwdDiff (n : ℕ) (a : Seq) :
diffIter n a = (fwdDiff (1 : ℤ))^[n] a := by
induction n with
| zero => rfl
| succ n ih =>
rw [diffIter, ih, Function.iterate_succ_apply']
rfl
/-- A forward difference is bounded by `2^n` times a bound for its
`n+1` input values. -/
theorem abs_diffIter_le_of_block (n : ℕ) (a : Seq) (j : ℤ) (M : ℝ)
(hM : 0 ≤ M)
(h : ∀ s : ℕ, s ≤ n → |a (j + s)| ≤ M) :
|diffIter n a j| ≤ (2 : ℝ) ^ n * M := by
induction n generalizing j with
| zero =>
simpa [diffIter] using h 0 (by omega)
| succ n ih =>
rw [diffIter, diff]
calc
|diffIter n a (j + 1) - diffIter n a j|
≤ |diffIter n a (j + 1)| + |diffIter n a j| :=
abs_sub _ _
_ ≤ 2 ^ n * M + 2 ^ n * M := by
gcongr
· apply ih
intro s hs
convert h (s + 1) (by omega) using 1 <;> push_cast <;> ring
· exact ih j (fun s hs => h s (by omega))
_ = 2 ^ (n + 1) * M := by ring
/-- Isolate the leftmost value in a forward difference. -/
theorem abs_head_le_diff_add (n : ℕ) (a : Seq) (j : ℤ) (M : ℝ)
(hM : 0 ≤ M)
(h : ∀ s : ℕ, 1 ≤ s → s ≤ n → |a (j + s)| ≤ M) :
|a j| ≤ |diffIter n a j| + ((2 : ℝ) ^ n - 1) * M := by
induction n generalizing j with
| zero => simp [diffIter]
| succ n ih =>
have htail :
|diffIter n a (j + 1)| ≤ (2 : ℝ) ^ n * M := by
apply abs_diffIter_le_of_block n a (j + 1) M hM
intro s hs
convert h (s + 1) (by omega) (by omega) using 1 <;> push_cast <;> ring
have hold :
|a j| ≤ |diffIter n a j| + ((2 : ℝ) ^ n - 1) * M := by
apply ih j
intro s hs0 hsn
exact h s hs0 (by omega)
have hrec :
|diffIter n a j| ≤
|diffIter (n + 1) a j| + |diffIter n a (j + 1)| := by
rw [diffIter, diff]
calc
|diffIter n a j| =
|(diffIter n a j - diffIter n a (j + 1)) +
diffIter n a (j + 1)| := by ring_nf
_ ≤ |diffIter n a j - diffIter n a (j + 1)| +
|diffIter n a (j + 1)| := abs_add_le _ _
_ = |diffIter n a (j + 1) - diffIter n a j| +
|diffIter n a (j + 1)| := by rw [abs_sub_comm]
calc
|a j| ≤ |diffIter n a j| + (2 ^ n - 1) * M := hold
_ ≤ (|diffIter (n + 1) a j| + |diffIter n a (j + 1)|) +
(2 ^ n - 1) * M := by gcongr
_ ≤ (|diffIter (n + 1) a j| + 2 ^ n * M) +
(2 ^ n - 1) * M := by gcongr
_ = |diffIter (n + 1) a j| + (2 ^ (n + 1) - 1) * M := by ring
/-- Isolate the rightmost value in a forward difference. -/
theorem abs_tail_le_diff_add (n : ℕ) (a : Seq) (j : ℤ) (M : ℝ)
(hM : 0 ≤ M)
(h : ∀ s : ℕ, s < n → |a (j + s)| ≤ M) :
|a (j + n)| ≤ |diffIter n a j| + ((2 : ℝ) ^ n - 1) * M := by
induction n generalizing j with
| zero => simp [diffIter]
| succ n ih =>
have hleft :
|diffIter n a j| ≤ (2 : ℝ) ^ n * M := by
apply abs_diffIter_le_of_block n a j M hM
intro s hs
exact h s (by omega)
have hold :
|a ((j + 1) + n)| ≤ |diffIter n a (j + 1)| +
((2 : ℝ) ^ n - 1) * M := by
apply ih (j + 1)
intro s hs
convert h (s + 1) (by omega) using 1 <;> push_cast <;> ring
have hrec :
|diffIter n a (j + 1)| ≤
|diffIter (n + 1) a j| + |diffIter n a j| := by
rw [diffIter, diff]
calc
|diffIter n a (j + 1)| =
|(diffIter n a (j + 1) - diffIter n a j) +
diffIter n a j| := by ring_nf
_ ≤ |diffIter n a (j + 1) - diffIter n a j| +
|diffIter n a j| := abs_add_le _ _
rw [show j + (n + 1 : ℕ) = (j + 1) + n by push_cast; ring]
calc
|a ((j + 1) + n)| ≤
|diffIter n a (j + 1)| + (2 ^ n - 1) * M := hold
_ ≤ (|diffIter (n + 1) a j| + |diffIter n a j|) +
(2 ^ n - 1) * M := by gcongr
_ ≤ (|diffIter (n + 1) a j| + 2 ^ n * M) +
(2 ^ n - 1) * M := by gcongr
_ = |diffIter (n + 1) a j| + (2 ^ (n + 1) - 1) * M := by ring
/-! ## Stencil intervals and local amplitudes -/
/-- The paper's recursively generated local-amplitude constant. -/
def localConstant : ℕ → ℝ
| 0 => 1
| n + 1 =>
if n = 0 then 1 else
((2 : ℝ) ^ n - 1) * localConstant n
theorem localConstant_succ {n : ℕ} (hn : 1 ≤ n) :
localConstant (n + 1) =
((2 : ℝ) ^ n - 1) * localConstant n := by
rw [localConstant, if_neg (by omega)]
theorem localConstant_nonneg (n : ℕ) : 0 ≤ localConstant n := by
induction n with
| zero => simp [localConstant]
| succ n ih =>
by_cases hn : n = 0
· simp [localConstant, hn]
· rw [localConstant, if_neg hn]
exact mul_nonneg (sub_nonneg.mpr (one_le_pow₀ (by norm_num))) ih
theorem localConstant_one : localConstant 1 = 1 := by
simp [localConstant]
theorem localConstant_le_succ {n : ℕ} (hn : 1 ≤ n) :
localConstant n ≤ localConstant (n + 1) := by
rw [localConstant_succ hn]
have hfactor : (1 : ℝ) ≤ 2 ^ n - 1 := by
have : (2 : ℝ) ^ 1 ≤ 2 ^ n :=
pow_le_pow_right₀ (by norm_num) hn
norm_num at this
linarith
exact (le_mul_of_one_le_left (localConstant_nonneg n) hfactor)
theorem enoLeft_monotone (u : Seq) (ell : ℕ) :
Monotone fun i : ℤ => enoLeft u i ell := by
intro i q hiq
apply Int.le_induction (motive := fun n _ =>
enoLeft u i ell ≤ enoLeft u n ell)
· exact le_rfl
· intro n hin ih
exact ih.trans (enoLeft_mono_step u n ell)
· exact hiq
theorem selectedInterval_nonempty_prev {u : Seq} {i : ℤ} {ell : ℕ}
(hnew : enoLeft u i (ell + 1) <
enoLeft u (i + 1) (ell + 1)) :
enoLeft u i ell < enoLeft u (i + 1) ell := by
by_contra h
have heq : enoLeft u i ell = enoLeft u (i + 1) ell := by
exact le_antisymm (enoLeft_mono_step u i ell) (not_lt.mp h)
simp [enoLeft, heq] at hnew
theorem enoLeft_succ_cases (u : Seq) (i : ℤ) (n : ℕ) :
enoLeft u i (n + 1) = enoLeft u i n - 1 ∨
enoLeft u i (n + 1) = enoLeft u i n := by
simp only [enoLeft]
split_ifs <;> simp
theorem enoLeft_succ_moved_iff (u : Seq) (i : ℤ) (n : ℕ) :
enoLeft u i (n + 1) = enoLeft u i n - 1 ↔
|diffIter n (jumps u) (enoLeft u i n - 1)| <
|diffIter n (jumps u) (enoLeft u i n)| := by
simp only [enoLeft]
split_ifs with h
· exact ⟨fun _ => h, fun _ => rfl⟩
· constructor
· omega
· exact fun h' => (h h').elim
theorem enoLeft_succ_stayed_iff (u : Seq) (i : ℤ) (n : ℕ) :
enoLeft u i (n + 1) = enoLeft u i n ↔
¬ |diffIter n (jumps u) (enoLeft u i n - 1)| <
|diffIter n (jumps u) (enoLeft u i n)| := by
simp only [enoLeft]
split_ifs with h
· constructor
· omega
· exact fun h' => (h' h).elim
· exact ⟨fun _ => h, fun _ => rfl⟩
/-- Pointwise form of Lemma 2.3: every entry of the selected length-`ell`
block is controlled by the interface jump selecting that block. -/
theorem selected_block_bound
(ell : ℕ) (hell : 1 ≤ ell) (u : Seq) (i j : ℤ)
(hj : enoLeft u i (ell - 1) ≤ j ∧
j < enoLeft u (i + 1) (ell - 1)) :
∀ s : ℕ, s < ell →
|jumps u (j + s)| ≤ localConstant ell * |jumps u i| := by
induction ell, hell using Nat.le_induction generalizing i j with
| base =>
have hji : j = i := by
simp [enoLeft] at hj
omega
subst j
intro s hs
have : s = 0 := by omega
subst s
simp [localConstant_one]
| succ ell hell ih =>
intro s hs
let a := jumps u
let L := enoLeft u i (ell - 1)
let R := enoLeft u (i + 1) (ell - 1)
change |a (j + s)| ≤ localConstant (ell + 1) * |a i|
have hell0 : 0 < ell := by omega
have hlevel : ell + 1 - 1 = ell := by omega
rw [hlevel] at hj
have hLR : L < R := by
have hstrict :
enoLeft u i ell < enoLeft u (i + 1) ell :=
lt_of_le_of_lt hj.1 hj.2
have hp := selectedInterval_nonempty_prev
(ell := ell - 1) (u := u) (i := i)
(by simpa [Nat.sub_add_cancel hell] using hstrict)
exact hp
have hM : 0 ≤ localConstant ell * |a i| :=
mul_nonneg (localConstant_nonneg ell) (abs_nonneg _)
by_cases hjL : j < L
· have hjEq : j = L - 1 := by
have hc := enoLeft_succ_cases u i (ell - 1)
rw [Nat.sub_add_cancel hell] at hc
change enoLeft u i ell = L - 1 ∨ enoLeft u i ell = L at hc
omega
have hmove :
|diffIter (ell - 1) a (L - 1)| <
|diffIter (ell - 1) a L| := by
have hc := enoLeft_succ_cases u i (ell - 1)
rw [Nat.sub_add_cancel hell] at hc
change enoLeft u i ell = L - 1 ∨ enoLeft u i ell = L at hc
have heq : enoLeft u i ell = L - 1 := by omega
have hm := (enoLeft_succ_moved_iff u i (ell - 1)).mp
(by simpa [Nat.sub_add_cancel hell, L] using heq)
simpa [a, L] using hm
have hold : ∀ q : ℕ, q < ell →
|a (L + q)| ≤ localConstant ell * |a i| :=
ih i L ⟨le_rfl, hLR⟩
have hdiff :
|diffIter (ell - 1) a L| ≤
2 ^ (ell - 1) * (localConstant ell * |a i|) := by
apply abs_diffIter_le_of_block _ _ _ _ hM
intro q hq
exact hold q (by omega)
rw [hjEq]
by_cases hs0 : s = 0
· subst s
have hhead := abs_head_le_diff_add (ell - 1) a (L - 1)
(localConstant ell * |a i|) hM
(fun q hq0 hq =>
by
have hv := hold (q - 1) (by omega)
convert hv using 1
congr 2
rw [Nat.cast_sub hq0]
push_cast
ring)
calc
|a (L - 1 + (0 : ℕ))| = |a (L - 1)| := by norm_num
_ ≤ |diffIter (ell - 1) a (L - 1)| +
(2 ^ (ell - 1) - 1) *
(localConstant ell * |a i|) := hhead
_ ≤ 2 ^ (ell - 1) * (localConstant ell * |a i|) +
(2 ^ (ell - 1) - 1) *
(localConstant ell * |a i|) := by gcongr; exact hmove.le.trans hdiff
_ = localConstant (ell + 1) * |a i| := by
rw [localConstant_succ hell]
have hpow : (2 : ℝ) ^ ell =
2 * (2 : ℝ) ^ (ell - 1) := by
conv_lhs => rw [show ell = (ell - 1) + 1 by omega, pow_succ]
ring
rw [hpow]
ring
· have hsq : s - 1 < ell := by omega
have hsge : 1 ≤ s := by omega
have hv := (hold (s - 1) hsq).trans
(mul_le_mul_of_nonneg_right (localConstant_le_succ hell)
(abs_nonneg _))
convert hv using 1
congr 2
rw [Nat.cast_sub hsge]
push_cast
ring
· by_cases hjR : j < R - 1
· have hjold : L ≤ j ∧ j < R := ⟨le_of_not_gt hjL, by omega⟩
have hj1old : L ≤ j + 1 ∧ j + 1 < R :=
⟨by omega, by omega⟩
by_cases hsell : s < ell
· exact (ih i j hjold s hsell).trans
(mul_le_mul_of_nonneg_right (localConstant_le_succ hell)
(abs_nonneg _))
· have hsEq : s = ell := by omega
subst s
have hv := (ih i (j + 1) hj1old (ell - 1) (by omega)).trans
(mul_le_mul_of_nonneg_right (localConstant_le_succ hell)
(abs_nonneg _))
have heq : j + (ell : ℤ) =
(j + 1) + ((ell - 1 : ℕ) : ℤ) := by
rw [Nat.cast_sub hell]
push_cast
ring
simpa [heq] using hv
· have hjEq : j = R - 1 := by
have hright : j < R := by
have hc := enoLeft_succ_cases u (i + 1) (ell - 1)
rw [Nat.sub_add_cancel hell] at hc
change enoLeft u (i + 1) ell = R - 1 ∨
enoLeft u (i + 1) ell = R at hc
omega
omega
have hstay : enoLeft u (i + 1) ell = R := by
have hc := enoLeft_succ_cases u (i + 1) (ell - 1)
rw [Nat.sub_add_cancel hell] at hc
change enoLeft u (i + 1) ell = R - 1 ∨
enoLeft u (i + 1) ell = R at hc
omega
have hcomp :
|diffIter (ell - 1) a R| ≤
|diffIter (ell - 1) a (R - 1)| := by
have hs' := (enoLeft_succ_stayed_iff u (i + 1) (ell - 1)).mp
(by simpa [Nat.sub_add_cancel hell, R] using hstay)
simpa [a, R] using le_of_not_gt hs'
have hold : ∀ q : ℕ, q < ell →
|a (R - 1 + q)| ≤ localConstant ell * |a i| :=
ih i (R - 1) ⟨by omega, by omega⟩
have hdiff :
|diffIter (ell - 1) a (R - 1)| ≤
2 ^ (ell - 1) * (localConstant ell * |a i|) := by
apply abs_diffIter_le_of_block _ _ _ _ hM
intro q hq
exact hold q (by omega)
rw [hjEq]
by_cases hsLast : s = ell
· subst s
have htail := abs_tail_le_diff_add (ell - 1) a R
(localConstant ell * |a i|) hM
(fun q hq =>
by
convert hold (q + 1) (by omega) using 1 <;>
push_cast <;> ring)
calc
|a (R - 1 + (ell : ℕ))| =
|a (R + (ell - 1 : ℕ))| := by
congr 2
rw [Nat.cast_sub hell]
push_cast
ring
_ ≤ |diffIter (ell - 1) a R| +
(2 ^ (ell - 1) - 1) *
(localConstant ell * |a i|) := htail
_ ≤ 2 ^ (ell - 1) * (localConstant ell * |a i|) +
(2 ^ (ell - 1) - 1) *
(localConstant ell * |a i|) := by
gcongr
exact hcomp.trans hdiff
_ = localConstant (ell + 1) * |a i| := by
rw [localConstant_succ hell]
have hpow : (2 : ℝ) ^ ell =
2 * (2 : ℝ) ^ (ell - 1) := by
conv_lhs => rw [show ell = (ell - 1) + 1 by omega, pow_succ]
ring
rw [hpow]
ring
· exact (hold s (by omega)).trans (by
gcongr
exact localConstant_le_succ hell)
theorem local_amplitude_bound {k : ℕ} (hk : 1 ≤ k) (u : Seq) (i j : ℤ)
(hj : selectedInterval k u i j) :
blockAmplitude k (jumps u) j ≤ localConstant k * |jumps u i| := by
unfold blockAmplitude
simp only [dif_pos (show 0 < k by omega)]
apply Finset.sup'_le
intro s hs
simp only [Finset.mem_range] at hs
exact selected_block_bound k hk u i j hj s hs
/-! ## The unique owner of a selected difference -/
theorem owner_offset_exists {k : ℕ} (hk : 1 ≤ k) (u : Seq) (j : ℤ) :
∃ n : ℕ, j <
enoLeft u (j + (n : ℤ) + 1) (k - 1) := by
refine ⟨k - 1, ?_⟩
have hb := (enoLeft_bounds u
(j + ((k - 1 : ℕ) : ℤ) + 1) (k - 1)).1
push_cast at hb
omega
/-- The unique interface whose final endpoint interval contains `j`. -/
def selectedOwner (k : ℕ) (u : Seq) (j : ℤ) : ℤ :=
if hk : 1 ≤ k then
j + (Nat.find (owner_offset_exists hk u j) : ℤ)
else j
theorem selectedOwner_mem {k : ℕ} (hk : 1 ≤ k) (u : Seq) (j : ℤ) :
selectedInterval k u (selectedOwner k u j) j := by
let n := Nat.find (owner_offset_exists hk u j)
have hn := Nat.find_spec (owner_offset_exists hk u j)
have hright :
j < enoLeft u (j + (n : ℤ) + 1) (k - 1) := hn
have hleft :
enoLeft u (j + (n : ℤ)) (k - 1) ≤ j := by
by_cases hn0 : n = 0
· rw [hn0]
simpa using (enoLeft_bounds u j (k - 1)).2
· have hmin := Nat.find_min
(owner_offset_exists hk u j) (show n - 1 < n by omega)
simp only [not_lt] at hmin
have heq :
j + (((n - 1 : ℕ) : ℤ)) + 1 = j + (n : ℤ) := by
rw [Nat.cast_sub (show 1 ≤ n by omega)]
push_cast
ring
simpa [heq] using hmin
simp only [selectedOwner, dif_pos hk, selectedInterval]
constructor
· exact hleft
· convert hright using 1 <;> ring
theorem selectedInterval_owner_unique {k : ℕ} (hk : 1 ≤ k) (u : Seq)
(i j : ℤ) (hj : selectedInterval k u i j) :
selectedOwner k u j = i := by
have ho := selectedOwner_mem hk u j
by_contra hne
by_cases hoi : selectedOwner k u j < i
· have hm := enoLeft_monotone u (k - 1)
(show selectedOwner k u j + 1 ≤ i by omega)
exact (not_lt_of_ge hm) (lt_of_le_of_lt hj.1 ho.2)
· have hio : i < selectedOwner k u j := by omega
have hm := enoLeft_monotone u (k - 1)
(show i + 1 ≤ selectedOwner k u j by omega)
exact (not_lt_of_ge hm) (lt_of_le_of_lt ho.1 hj.2)
end ENOTV