Semidiscrete wave error in displacement and velocity

import PDETextbook.Base

open MeasureTheory
open scoped BigOperators ENNReal

noncomputable section

namespace PDETextbook.LarssonThomee.Thm13_1

/-- Closed triangle spanned by three vertices. -/
def triSet (K : Fin 3 → Euclidean 2) : Set (Euclidean 2) := convexHull ℝ (Set.range K)

/-- Ω ⊂ ℝ² is a bounded convex domain whose boundary Γ is a polygon (§13.1 setting). -/
def IsConvexPolygonalDomain (Ω : Set (Euclidean 2)) : Prop :=
  IsOpen Ω ∧ Ω.Nonempty ∧ Convex ℝ Ω ∧ Bornology.IsBounded Ω ∧
    ∃ P : Finset (Euclidean 2), closure Ω = convexHull ℝ (P : Set (Euclidean 2))

/-- Triangulation 𝒯_h of Ω (closed triangles, conforming), with h = max_K diam K. -/
def IsTriangulation (Ω : Set (Euclidean 2)) (T : Finset (Fin 3 → Euclidean 2)) (h : ℝ) :
    Prop :=
  T.Nonempty ∧ (∀ K ∈ T, AffineIndependent ℝ K) ∧
  closure Ω = (⋃ K ∈ T, triSet K) ∧
  (∀ K ∈ T, Metric.diam (triSet K) ≤ h) ∧ (∃ K ∈ T, Metric.diam (triSet K) = h) ∧
  (∀ K ∈ T, ∀ K' ∈ T, triSet K ≠ triSet K' →
    triSet K ∩ triSet K' = ∅ ∨
    (∃ p, p ∈ Set.range K ∧ p ∈ Set.range K' ∧ triSet K ∩ triSet K' = {p}) ∨
    (∃ p q, p ≠ q ∧ p ∈ Set.range K ∧ q ∈ Set.range K ∧ p ∈ Set.range K' ∧
      q ∈ Set.range K' ∧ triSet K ∩ triSet K' = segment ℝ p q)) ∧
  (∀ K ∈ T, ∀ K' ∈ T, ∀ i j k : Fin 3, i ≠ j → K' k ∉ openSegment ℝ (K i) (K j))

/-- Angles of all triangles of all members of {𝒯_h}_{0<h<1} are bounded below uniformly in h. -/
def MinAngleBounded (T : ℝ → Finset (Fin 3 → Euclidean 2)) : Prop :=
  ∃ θ₀ : ℝ, 0 < θ₀ ∧ ∀ h ∈ Set.Ioo (0 : ℝ) 1, ∀ K ∈ T h, ∀ i j k : Fin 3,
    i ≠ j → i ≠ k → j ≠ k → θ₀ ≤ EuclideanGeometry.angle (K j) (K i) (K k)

/-- S_h: continuous piecewise linear functions on 𝒯_h vanishing on Γ. -/
def IsShFun (Ω : Set (Euclidean 2)) (T : Finset (Fin 3 → Euclidean 2))
    (g : Euclidean 2 → ℝ) : Prop :=
  ContinuousOn g (closure Ω) ∧
  (∀ K ∈ T, ∃ A : Euclidean 2 →ᵃ[ℝ] ℝ, Set.EqOn g A (triSet K)) ∧
  ∀ x ∈ frontier Ω, g x = 0

/-- An L² class belongs to S_h. -/
def InSh (Ω : Set (Euclidean 2)) (T : Finset (Fin 3 → Euclidean 2)) (χ : L2 ℝ Ω) : Prop :=
  ∃ g, IsShFun Ω T g ∧ (⇑χ) =ᵐ[domainMeasure Ω] g

/-- c = a(v,w) = (∇v, ∇w), via weak gradients (a.e. unique). -/
def HasEnergyProduct (Ω : Set (Euclidean 2)) (v w : Euclidean 2 → ℝ) (c : ℝ) : Prop :=
  ∃ G H, HasH1Gradient Ω v G ∧ HasH1Gradient Ω w H ∧
    c = ∑ i, ∫ x in Ω, G i x * H i x

/-- Elliptic (Ritz) projection (5.49): R_h v ∈ S_h, a(R_h v, χ) = a(v, χ) ∀ χ ∈ S_h. -/
def IsRitzProjection (Ω : Set (Euclidean 2)) (T : Finset (Fin 3 → Euclidean 2))
    (v r : L2 ℝ Ω) : Prop :=
  InSh Ω T r ∧ ∀ χ, InSh Ω T χ → ∃ c, HasEnergyProduct Ω (⇑r) (⇑χ) c ∧
    HasEnergyProduct Ω (⇑v) (⇑χ) c

/-- Seminorm |v|_1 (⊤ if v ∉ H¹). -/
def h1Seminorm {d : ℕ} (Ω : Set (Euclidean d)) (v : Euclidean d → ℝ) : ℝ≥0∞ :=
  ⨅ (G : Fin d → Euclidean d → ℝ) (_ : HasH1Gradient Ω v G),
    (∑ i, eLpNorm (G i) 2 (domainMeasure Ω) ^ 2) ^ ((1 : ℝ) / 2)

/-- Norm ‖v‖_1 (⊤ if v ∉ H¹). -/
def h1Norm {d : ℕ} (Ω : Set (Euclidean d)) (v : Euclidean d → ℝ) : ℝ≥0∞ :=
  ⨅ (G : Fin d → Euclidean d → ℝ) (_ : HasH1Gradient Ω v G),
    (eLpNorm v 2 (domainMeasure Ω) ^ 2 +
      ∑ i, eLpNorm (G i) 2 (domainMeasure Ω) ^ 2) ^ ((1 : ℝ) / 2)

/-- Norm ‖v‖_2 = (Σ_{|α|≤2} ‖D^α v‖²)^{1/2}; each second-order multi-index counted once (i ≤ j). -/
def h2Norm {d : ℕ} (Ω : Set (Euclidean d)) (v : Euclidean d → ℝ) : ℝ≥0∞ :=
  ⨅ (G : Fin d → Euclidean d → ℝ) (H : Fin d → Fin d → Euclidean d → ℝ)
    (_ : HasH1Gradient Ω v G ∧ ∀ i, HasH1Gradient Ω (G i) (H i)),
    (eLpNorm v 2 (domainMeasure Ω) ^ 2 +
      ∑ i, eLpNorm (G i) 2 (domainMeasure Ω) ^ 2 +
      ∑ p ∈ Finset.univ.filter (fun p : Fin d × Fin d => p.1 ≤ p.2),
        eLpNorm (H p.1 p.2) 2 (domainMeasure Ω) ^ 2) ^ ((1 : ℝ) / 2)

/-- u solves (13.1): u_tt - Δu = f in Ω × ℝ₊ (weak form), u = 0 on Γ (u(t) ∈ H¹₀),
u(0) = v, u_t(0) = w; U', U'' are the L²-valued time derivatives. -/
def IsWaveSolution (Ω : Set (Euclidean 2)) (f : ℝ → L2 ℝ Ω) (v w : L2 ℝ Ω)
    (U U' U'' : ℝ → L2 ℝ Ω) : Prop :=
  (∀ t ≥ 0, HasDerivWithinAt U (U' t) (Set.Ici 0) t ∧
    HasDerivWithinAt U' (U'' t) (Set.Ici 0) t) ∧
  U 0 = v ∧ U' 0 = w ∧
  ∀ t > 0, InH10Closure Ω (⇑(U t)) ∧ ∃ G, HasH1Gradient Ω (⇑(U t)) G ∧
    ∀ φ : Euclidean 2 → ℝ, CompactC1 Ω φ →
      (∫ x in Ω, (U'' t : Euclidean 2 → ℝ) x * φ x) +
        ∑ i, (∫ x in Ω, G i x * partialDeriv φ i x) =
      ∫ x in Ω, (f t : Euclidean 2 → ℝ) x * φ x

/-- u_h solves (13.2): u_h(t) ∈ S_h, (u_{h,tt}, χ) + a(u_h, χ) = (f, χ) ∀ χ ∈ S_h, t > 0,
u_h(0) = v_h, u_{h,t}(0) = w_h. -/
def IsSemidiscreteSolution (Ω : Set (Euclidean 2)) (T : Finset (Fin 3 → Euclidean 2))
    (f : ℝ → L2 ℝ Ω) (vh wh : L2 ℝ Ω) (Uh Uh' Uh'' : ℝ → L2 ℝ Ω) : Prop :=
  (∀ t ≥ 0, InSh Ω T (Uh t) ∧ HasDerivWithinAt Uh (Uh' t) (Set.Ici 0) t ∧
    HasDerivWithinAt Uh' (Uh'' t) (Set.Ici 0) t) ∧
  Uh 0 = vh ∧ Uh' 0 = wh ∧
  ∀ t > 0, ∀ χ, InSh Ω T χ → ∃ c, HasEnergyProduct Ω (⇑(Uh t)) (⇑χ) c ∧
    l2Inner (Uh'' t) χ + c = l2Inner (f t) χ

/-- Theorem 13.1 (Larsson–Thomée). C depends only on Ω and the triangulation family,
not on h, t or the data/solutions. Norms are ENNReal-valued (⊤ outside the space). -/
theorem thm_13_1 (Ω : Set (Euclidean 2)) (hΩ : IsConvexPolygonalDomain Ω)
    (T : ℝ → Finset (Fin 3 → Euclidean 2))
    (hT : ∀ h ∈ Set.Ioo (0 : ℝ) 1, IsTriangulation Ω (T h) h)
    (hang : MinAngleBounded T) :
    ∃ C : NNReal, ∀ h ∈ Set.Ioo (0 : ℝ) 1,
      ∀ (f : ℝ → L2 ℝ Ω) (v w : L2 ℝ Ω) (U U' U'' : ℝ → L2 ℝ Ω),
      IsWaveSolution Ω f v w U U' U'' →
      ∀ (vh wh : L2 ℝ Ω) (Uh Uh' Uh'' : ℝ → L2 ℝ Ω),
      IsSemidiscreteSolution Ω (T h) f vh wh Uh Uh' Uh'' →
      ∀ Rv Rw : L2 ℝ Ω, IsRitzProjection Ω (T h) v Rv → IsRitzProjection Ω (T h) w Rw →
      ∀ t : ℝ, 0 ≤ t →
        ENNReal.ofReal ‖Uh' t - U' t‖ ≤
            (C : ℝ≥0∞) * (h1Seminorm Ω (⇑(vh - Rv)) + ENNReal.ofReal ‖wh - Rw‖) +
            (C : ℝ≥0∞) * ENNReal.ofReal h ^ 2 *
              (h2Norm Ω (⇑(U' t)) + ∫⁻ s in Set.Icc 0 t, h2Norm Ω (⇑(U'' s))) ∧
        ENNReal.ofReal ‖Uh t - U t‖ ≤
            (C : ℝ≥0∞) * (h1Seminorm Ω (⇑(vh - Rv)) + ENNReal.ofReal ‖wh - Rw‖) +
            (C : ℝ≥0∞) * ENNReal.ofReal h ^ 2 *
              (h2Norm Ω (⇑(U t)) + ∫⁻ s in Set.Icc 0 t, h2Norm Ω (⇑(U'' s))) ∧
        h1Seminorm Ω (⇑(Uh t - U t)) ≤
            (C : ℝ≥0∞) * (h1Seminorm Ω (⇑(vh - Rv)) + ENNReal.ofReal ‖wh - Rw‖) +
            (C : ℝ≥0∞) * ENNReal.ofReal h *
              (h2Norm Ω (⇑(U t)) + ∫⁻ s in Set.Icc 0 t, h1Norm Ω (⇑(U'' s))) := by
  sorry

end PDETextbook.LarssonThomee.Thm13_1