Discrete wave-energy conservation at half steps

import PDETextbook.Base

open MeasureTheory
open scoped BigOperators ENNReal

noncomputable section

namespace PDETextbook.LarssonThomee.Lem13_2

/-- Sect. 13.1 standing assumption: `Ω ⊂ ℝ²` is a bounded convex domain whose boundary `Γ`
is a (closed) polygon with vertices `p 0, …, p (m+2)`. -/
def IsConvexPolygonalDomain (Ω : Set (Euclidean 2)) : Prop :=
  IsOpen Ω ∧ Ω.Nonempty ∧ Bornology.IsBounded Ω ∧ Convex ℝ Ω ∧
  ∃ (m : ℕ) (p : Fin (m + 3) → Euclidean 2),
    frontier Ω = ⋃ i : Fin (m + 3), segment ℝ (p i) (p (i + 1))

/-- A (conforming) triangulation of `Ω̄`: finitely many nondegenerate triangles covering `Ω̄`,
with pairwise disjoint interiors, any two meeting in the hull of their common vertices. -/
structure Triangulation (Ω : Set (Euclidean 2)) where
  tri : Finset (Fin 3 → Euclidean 2)
  nondegenerate : ∀ K ∈ tri, AffineIndependent ℝ K
  cover : (⋃ K ∈ tri, convexHull ℝ (Set.range K)) = closure Ω
  disjoint_interiors : ∀ K ∈ tri, ∀ K' ∈ tri, K ≠ K' →
    Disjoint (interior (convexHull ℝ (Set.range K)))
      (interior (convexHull ℝ (Set.range K')))
  conforming : ∀ K ∈ tri, ∀ K' ∈ tri,
    convexHull ℝ (Set.range K) ∩ convexHull ℝ (Set.range K') =
      convexHull ℝ (Set.range K ∩ Set.range K')

/-- `S_h ⊂ H₀¹`: continuous piecewise linear functions on the triangulation, vanishing on `Γ`
(represented globally by zero extension outside `Ω`). -/
def FESpace {Ω : Set (Euclidean 2)} (𝒯 : Triangulation Ω) : Set (Euclidean 2 → ℝ) :=
  {v | Continuous v ∧ (∀ x, x ∉ Ω → v x = 0) ∧
    ∀ K ∈ 𝒯.tri, ∃ g : Euclidean 2 →ᵃ[ℝ] ℝ, Set.EqOn v g (convexHull ℝ (Set.range K))}

/-- The `L₂(Ω)` inner product `(v, w) = ∫_Ω v w dx`. -/
def innerΩ (Ω : Set (Euclidean 2)) (v w : Euclidean 2 → ℝ) : ℝ :=
  ∫ x in Ω, v x * w x

/-- `‖v‖ = (∫_Ω v² dx)^{1/2}`, the `L₂(Ω)` norm. -/
def normΩ (Ω : Set (Euclidean 2)) (v : Euclidean 2 → ℝ) : ℝ :=
  Real.sqrt (∫ x in Ω, v x ^ 2)

/-- `a(v, w) = (∇v, ∇w) = ∫_Ω ∇v · ∇w dx`; for piecewise linear functions the gradient exists
classically on triangle interiors, i.e. a.e. in `Ω`, and coincides with the weak gradient there. -/
def aForm (Ω : Set (Euclidean 2)) (v w : Euclidean 2 → ℝ) : ℝ :=
  ∫ x in Ω, ∑ i : Fin 2, partialDeriv (𝕜 := ℝ) v i x * partialDeriv (𝕜 := ℝ) w i x

/-- `|v|₁ = (∑_{|α|=1} ‖D^α v‖²)^{1/2} = ‖∇v‖`, the `H¹` seminorm. -/
def seminorm1 (Ω : Set (Euclidean 2)) (v : Euclidean 2 → ℝ) : ℝ :=
  Real.sqrt (∫ x in Ω, ∑ i : Fin 2, (partialDeriv (𝕜 := ℝ) v i x) ^ 2)

/-- Forward difference quotient `∂_t Uⁿ = (Uⁿ⁺¹ - Uⁿ)/k`. -/
def fwdDiff (k : ℝ) (U : ℕ → Euclidean 2 → ℝ) (n : ℕ) : Euclidean 2 → ℝ :=
  fun x => (U (n + 1) x - U n x) / k

/-- The half-step value `U^{n+1/2} = (Uⁿ + Uⁿ⁺¹)/2 = (1/2)(Uⁿ + Uⁿ⁺¹)`, as defined in the
source just before Lemma 13.2 (`halfStep U n` is `U^{n+1/2}`; `halfStep U 0` is `U^{1/2}`). -/
def halfStep (U : ℕ → Euclidean 2 → ℝ) (n : ℕ) : Euclidean 2 → ℝ :=
  fun x => (1 / 2 : ℝ) * (U n x + U (n + 1) x)

/-- `U = (Uⁿ)` with `Uⁿ ∈ S_h` solves scheme (13.5) with time step `k`, `t_n = n k`:
`(∂_t ∂̄_t Uⁿ, χ) + a(¼(Uⁿ⁺¹ + 2Uⁿ + Uⁿ⁻¹), χ) = (f(t_n), χ)` for all `χ ∈ S_h`, `n ≥ 1`,
where `∂_t ∂̄_t Uⁿ = (Uⁿ⁺¹ - 2Uⁿ + Uⁿ⁻¹)/k²`; `U⁰`, `U¹` are the given starting values. -/
def SolvesScheme135 {Ω : Set (Euclidean 2)} (𝒯 : Triangulation Ω) (k : ℝ)
    (f : ℝ → Euclidean 2 → ℝ) (U : ℕ → Euclidean 2 → ℝ) : Prop :=
  (∀ n, U n ∈ FESpace 𝒯) ∧
  ∀ n : ℕ, 1 ≤ n → ∀ χ ∈ FESpace 𝒯,
    innerΩ Ω (fun x => (U (n + 1) x - 2 * U n x + U (n - 1) x) / k ^ 2) χ
      + aForm Ω (fun x => (U (n + 1) x + 2 * U n x + U (n - 1) x) / 4) χ
      = innerΩ Ω (f ((n : ℝ) * k)) χ

/-- Lemma 13.2 (Larsson–Thomée): energy conservation for the solution of (13.5) with `f = 0`:
`‖∂_t Uⁿ‖² + |U^{n+1/2}|₁² = ‖∂_t U⁰‖² + |U^{1/2}|₁²` for `n ≥ 0`. -/
theorem lem_13_2 (Ω : Set (Euclidean 2)) (hΩ : IsConvexPolygonalDomain Ω)
    (𝒯 : Triangulation Ω) (k : ℝ) (hk : 0 < k) (U : ℕ → Euclidean 2 → ℝ)
    (hU : SolvesScheme135 𝒯 k (fun _ _ => 0) U) :
    ∀ n : ℕ, normΩ Ω (fwdDiff k U n) ^ 2 + seminorm1 Ω (halfStep U n) ^ 2
      = normΩ Ω (fwdDiff k U 0) ^ 2 + seminorm1 Ω (halfStep U 0) ^ 2 := by
  sorry

end PDETextbook.LarssonThomee.Lem13_2