Backward Euler–Galerkin heat error

import PDETextbook.Base

open MeasureTheory
open scoped BigOperators ENNReal

noncomputable section

namespace PDETextbook.LarssonThomee.Thm10_5

open Classical

/-- Smooth boundary Γ of Ω, realized by a smooth defining function with nonvanishing gradient on Γ. -/
def HasSmoothBoundary (Ω : Set (Euclidean 2)) : Prop :=
  ∃ ρ : Euclidean 2 → ℝ, ContDiff ℝ (⊤ : ℕ∞) ρ ∧ Ω = {x | ρ x < 0} ∧
    ∀ x ∈ frontier Ω, fderiv ℝ ρ x ≠ 0

/-- Ω_h: union of the closed triangles (convex hulls of vertex triples). -/
def meshDomain (T : Finset (Finset (Euclidean 2))) : Set (Euclidean 2) :=
  ⋃ K ∈ T, convexHull ℝ (K : Set (Euclidean 2))

/-- Triangulation T_h of Ω with mesh size h (Sects. 5.2, 5.3, 10.1): conforming triangles,
angles bounded below by θ₀, triangles in Ω̄, boundary vertices of Ω_h on Γ, and
Ω \ Ω_h of width at most cw·h². -/
structure AdmissibleTriangulation (Ω : Set (Euclidean 2)) (θ₀ cw h : ℝ)
    (T : Finset (Finset (Euclidean 2))) : Prop where
  nonempty : T.Nonempty
  tri : ∀ K ∈ T, K.card = 3 ∧ AffineIndependent ℝ ((↑) : K → Euclidean 2)
  conforming : ∀ K ∈ T, ∀ K' ∈ T, K ≠ K' →
    convexHull ℝ (K : Set (Euclidean 2)) ∩ convexHull ℝ (K' : Set (Euclidean 2)) =
      convexHull ℝ ((K ∩ K' : Finset (Euclidean 2)) : Set (Euclidean 2))
  angle : ∀ K ∈ T, ∀ a ∈ K, ∀ b ∈ K, ∀ c ∈ K, a ≠ b → a ≠ c → b ≠ c →
    θ₀ ≤ EuclideanGeometry.angle b a c
  mesh_le : ∀ K ∈ T, Metric.diam (K : Set (Euclidean 2)) ≤ h
  mesh_eq : ∃ K ∈ T, Metric.diam (K : Set (Euclidean 2)) = h
  inside : ∀ K ∈ T, convexHull ℝ (K : Set (Euclidean 2)) ⊆ closure Ω
  bdry_nodes : ∀ K ∈ T, ∀ p ∈ K, p ∈ frontier (meshDomain T) → p ∈ frontier Ω
  width : ∀ x ∈ Ω, Metric.infDist x (meshDomain T) ≤ cw * h ^ 2

/-- S_h: continuous piecewise linear functions on T_h, vanishing outside Ω_h (hence on ∂Ω_h). -/
def InSh (T : Finset (Finset (Euclidean 2))) (χ : Euclidean 2 → ℝ) : Prop :=
  Continuous χ ∧ (∀ x ∉ meshDomain T, χ x = 0) ∧
    ∀ K ∈ T, ∃ A : Euclidean 2 →ᵃ[ℝ] ℝ, Set.EqOn χ A (convexHull ℝ (K : Set (Euclidean 2)))

/-- (v,w) = ∫_Ω v w dx. -/
def l2Pair (Ω : Set (Euclidean 2)) (w z : Euclidean 2 → ℝ) : ℝ := ∫ x in Ω, w x * z x

/-- a(v,w) = ∫_Ω ∇v·∇w dx. -/
def aForm (Ω : Set (Euclidean 2)) (w z : Euclidean 2 → ℝ) : ℝ :=
  ∫ x in Ω, ∑ i, partialDeriv w i x * partialDeriv z i x

/-- ‖v‖_2 = (Σ_{|α|≤2} ‖D^α v‖²)^{1/2} with weak derivatives; ⊤ if v ∉ H². Mixed α=(1,1) counted once. -/
def h2NormE (Ω : Set (Euclidean 2)) (w : Euclidean 2 → ℝ) : ℝ≥0∞ :=
  ⨅ (G : Fin 2 → Euclidean 2 → ℝ) (H : Fin 2 → Fin 2 → Euclidean 2 → ℝ)
    (_ : HasH1Gradient Ω w G ∧ ∀ i, HasH1Gradient Ω (G i) (H i)),
    (eLpNorm w 2 (domainMeasure Ω) ^ 2 + ∑ i, eLpNorm (G i) 2 (domainMeasure Ω) ^ 2 +
      ∑ i, ∑ j, (if i ≤ j then eLpNorm (H i j) 2 (domainMeasure Ω) ^ 2 else 0)) ^ (1 / 2 : ℝ)

/-- u_t (pointwise time derivative); used only where IsHeatSolution guarantees it exists (t>0, x∈Ω). -/
def timeDeriv (u : ℝ → Euclidean 2 → ℝ) : ℝ → Euclidean 2 → ℝ :=
  fun t x => deriv (fun s => u s x) t

/-- ∫_0^T ‖u_tt‖ ds as an extended value: infimum over genuine witnesses W with
W(s,x) = ∂_t u_t(s,x) for every s∈(0,T), x∈Ω. With no such witness the infimum is empty and
the value is ⊤, so no junk derivative value is ever used. Witnesses agree on (0,T)×Ω, so when
one exists this is exactly ∫_0^T ‖u_tt(s)‖ ds. -/
def uttIntegralE (Ω : Set (Euclidean 2)) (u : ℝ → Euclidean 2 → ℝ) (T : ℝ) : ℝ≥0∞ :=
  ⨅ (W : ℝ → Euclidean 2 → ℝ)
    (_ : ∀ s ∈ Set.Ioo 0 T, ∀ x ∈ Ω, HasDerivAt (fun r => timeDeriv u r x) (W s x) s),
    ∫⁻ s in Set.Ioo 0 T, eLpNorm (W s) 2 (domainMeasure Ω)

/-- Δ = ∂²/∂x₁² + ∂²/∂x₂². -/
def laplacian2 (w : Euclidean 2 → ℝ) (x : Euclidean 2) : ℝ :=
  ∑ i, partialDeriv (fun y => partialDeriv w i y) i x

/-- u solves (10.1) classically: u_t - Δu = f in Ω×R₊, u = 0 on Γ×R₊, u(·,0) = v in Ω. -/
def IsHeatSolution (Ω : Set (Euclidean 2)) (f u : ℝ → Euclidean 2 → ℝ) (v : Euclidean 2 → ℝ) : Prop :=
  ContinuousOn (fun p : ℝ × Euclidean 2 => u p.1 p.2) (Set.Ioi 0 ×ˢ closure Ω) ∧
  ContinuousOn (fun p : ℝ × Euclidean 2 => u p.1 p.2) (Set.Ici 0 ×ˢ Ω) ∧
  (∀ t > (0 : ℝ), ∀ x ∈ Ω, ContDiffAt ℝ 2 (u t) x ∧ DifferentiableAt ℝ (fun s => u s x) t ∧
      timeDeriv u t x - laplacian2 (u t) x = f t x) ∧
  (∀ t > (0 : ℝ), ∀ x ∈ frontier Ω, u t x = 0) ∧
  ∀ x ∈ Ω, u 0 x = v x

/-- Backward Euler–Galerkin scheme (10.20), t_n = n k. -/
def IsBackwardEulerGalerkin (Ω : Set (Euclidean 2)) (T : Finset (Finset (Euclidean 2))) (k : ℝ)
    (f : ℝ → Euclidean 2 → ℝ) (vh : Euclidean 2 → ℝ) (U : ℕ → Euclidean 2 → ℝ) : Prop :=
  (∀ n, InSh T (U n)) ∧ U 0 = vh ∧
  ∀ n : ℕ, 1 ≤ n → InL2 Ω (f ((n : ℝ) * k)) ∧
    ∀ χ, InSh T χ →
      l2Pair Ω (fun x => (U n x - U (n - 1) x) / k) χ + aForm Ω (U n) χ =
        l2Pair Ω (f ((n : ℝ) * k)) χ

/-- Theorem 10.5 (Larsson–Thomée). C depends on Ω, the angle bound θ₀, the width constant cw and
the constant C₀ of (10.10), but not on h, k, n, f, u, v. Norms of u_t, u_tt are extended-valued
(⊤ when the object does not exist), so the bound is trivial in that case as in the source.
No source warnings were listed. -/
theorem thm_10_5 (Ω : Set (Euclidean 2)) (hopen : IsOpen Ω) (hconv : Convex ℝ Ω)
    (hbdd : Bornology.IsBounded Ω) (hne : Ω.Nonempty) (hsmooth : HasSmoothBoundary Ω)
    (θ₀ cw C₀ : ℝ) (hθ : 0 < θ₀) :
    ∃ C : ℝ, ∀ (T : Finset (Finset (Euclidean 2))) (h : ℝ), 0 < h → h < 1 →
      AdmissibleTriangulation Ω θ₀ cw h T →
      ∀ k : ℝ, 0 < k →
      ∀ (f u : ℝ → Euclidean 2 → ℝ) (v vh : Euclidean 2 → ℝ) (U : ℕ → Euclidean 2 → ℝ),
        IsHeatSolution Ω f u v →
        InSh T vh →
        eLpNorm (fun x => vh x - v x) 2 (domainMeasure Ω) ≤
          ENNReal.ofReal (C₀ * h ^ 2) * h2NormE Ω v →
        IsBackwardEulerGalerkin Ω T k f vh U →
        ∀ n : ℕ,
          eLpNorm (fun x => U n x - u ((n : ℝ) * k) x) 2 (domainMeasure Ω) ≤
            ENNReal.ofReal (C * h ^ 2) *
                (h2NormE Ω v + ∫⁻ s in Set.Ioc 0 ((n : ℝ) * k), h2NormE Ω (timeDeriv u s)) +
              ENNReal.ofReal (C * k) * uttIntegralE Ω u ((n : ℝ) * k) := by
  sorry

end PDETextbook.LarssonThomee.Thm10_5