Crank–Nicolson finite-element heat error

import PDETextbook.Base

open MeasureTheory
open scoped BigOperators ENNReal

noncomputable section

namespace PDETextbook.LarssonThomee.Thm106

abbrev E2 := Euclidean 2

/-- Closed triangle with vertex triple `K` (Sect. 5.2). -/
def triSet (K : Fin 3 → E2) : Set E2 := convexHull ℝ (Set.range K)

/-- `Ω_h`, the union of the triangles. -/
def meshUnion (T : Finset (Fin 3 → E2)) : Set E2 := ⋃ K ∈ T, triSet K

/-- `h = max_K diam K`. -/
def meshSize (T : Finset (Fin 3 → E2)) : ℝ := ⨆ K ∈ T, Metric.diam (triSet K)

/-- Smooth boundary: locally a regular sublevel set of a C^∞ function. -/
def SmoothBoundary (Ω : Set E2) : Prop :=
  ∀ p ∈ frontier Ω, ∃ V ∈ nhds p, ∃ ρ : E2 → ℝ, ContDiff ℝ (⊤ : ℕ∞) ρ ∧
    (∀ x ∈ V, fderiv ℝ ρ x ≠ 0) ∧ Ω ∩ V = {x | x ∈ V ∧ ρ x < 0}

/-- Member of the triangulation family of Sects. 5.2 and 10: conforming triangles,
boundary vertices of `Ω_h` on `Γ`, `Ω \ Ω_h` of width at most `cw h²`,
all angles at least `θ₀`, and `h < 1`. -/
structure AdmissibleTriangulation (Ω : Set E2) (θ₀ cw : ℝ)
    (T : Finset (Fin 3 → E2)) : Prop where
  nondeg : ∀ K ∈ T, AffineIndependent ℝ K
  conforming : ∀ K ∈ T, ∀ K' ∈ T, triSet K ≠ triSet K' →
    triSet K ∩ triSet K' = convexHull ℝ (Set.range K ∩ Set.range K')
  subset : meshUnion T ⊆ closure Ω
  bdry_nodes : ∀ K ∈ T, ∀ i, K i ∈ frontier (meshUnion T) → K i ∈ frontier Ω
  width : ∀ x ∈ Ω \ meshUnion T, Metric.infDist x (frontier Ω) ≤ cw * meshSize T ^ 2
  min_angle : ∀ K ∈ T, ∀ i j l : Fin 3, i ≠ j → i ≠ l → j ≠ l →
    θ₀ ≤ EuclideanGeometry.angle (K j) (K i) (K l)
  mesh_lt_one : meshSize T < 1

/-- `S_h`: continuous, affine on each triangle, vanishing on `∂Ω_h` and on `Ω \ Ω_h`. -/
def InSh (T : Finset (Fin 3 → E2)) (χ : E2 → ℝ) : Prop :=
  Continuous χ ∧ (∀ K ∈ T, ∃ a : E2 →ᵃ[ℝ] ℝ, Set.EqOn χ a (triSet K)) ∧
    ∀ x, x ∉ interior (meshUnion T) → χ x = 0

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

/-- Classical Laplacian. -/
def lapC (w : E2 → ℝ) (x : E2) : ℝ :=
  ∑ i, partialDeriv (fun y => partialDeriv w i y) i x

/-- `‖w‖₂ = (Σ_{|α|≤2} ‖D^α w‖²)^{1/2}` with weak derivatives (`H i j = ∂_j∂_i`, each
multi-index once); `⊤` if `w ∉ H²`. -/
def h2Norm (Ω : Set E2) (w : E2 → ℝ) : ℝ≥0∞ :=
  ⨅ (G : Fin 2 → E2 → ℝ) (H : Fin 2 → Fin 2 → E2 → ℝ)
    (_ : (∀ i, HasWeakPartial Ω w i (G i)) ∧ ∀ i j, HasWeakPartial Ω (G i) j (H i j)),
    (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 : ℝ)

/-- Distributional Laplacian on `Ω`. -/
def IsWeakLaplacian (Ω : Set E2) (w g : E2 → ℝ) : Prop :=
  LocallyIntegrableOn w Ω volume ∧ LocallyIntegrableOn g Ω volume ∧
  ∀ φ : E2 → ℝ, ContDiff ℝ 2 φ → HasCompactSupport φ → tsupport φ ⊆ Ω →
    (∫ x in Ω, lapC φ x * w x) = ∫ x in Ω, φ x * g x

/-- `‖Δw‖`; `⊤` if `Δw ∉ L₂`. -/
def lapNorm (Ω : Set E2) (w : E2 → ℝ) : ℝ≥0∞ :=
  ⨅ (g : E2 → ℝ) (_ : IsWeakLaplacian Ω w g), eLpNorm g 2 (domainMeasure Ω)

/-- Classical solution of (10.1) with time derivative `ut`. -/
def IsSolution101 (Ω : Set E2) (f : ℝ → E2 → ℝ) (v : E2 → ℝ)
    (u ut : ℝ → E2 → ℝ) : Prop :=
  (∀ t > 0, ContDiffOn ℝ 2 (u t) Ω ∧ ContinuousOn (u t) (closure Ω)) ∧
  (∀ x ∈ Ω, ∀ t > 0, HasDerivAt (fun s => u s x) (ut t x) t) ∧
  (∀ t > 0, ∀ x ∈ Ω, ut t x - lapC (u t) x = f t x) ∧
  (∀ t > 0, ∀ x ∈ frontier Ω, u t x = 0) ∧
  (∀ x ∈ Ω, u 0 x = v x ∧ ContinuousWithinAt (fun s => u s x) (Set.Ici 0) 0)

/-- Crank–Nicolson scheme (10.25), `t_n = n k`. -/
def IsCN1025 (Ω : Set E2) (T : Finset (Fin 3 → E2)) (k : ℝ) (f : ℝ → E2 → ℝ)
    (vh : E2 → ℝ) (U : ℕ → E2 → ℝ) : Prop :=
  (∀ n, InSh T (U n)) ∧ U 0 = vh ∧
  ∀ n : ℕ, 1 ≤ n → ∀ χ, InSh T χ →
    (∫ x in Ω, (U n x - U (n - 1) x) / k * χ x) +
      aForm Ω (fun x => (U n x + U (n - 1) x) / 2) χ =
    ∫ x in Ω, f (((n : ℝ) - 1/2) * k) x * χ x

/-- Theorem 10.6 (Larsson–Thomée). `C` depends on `Ω`, the angle bound `θ₀`, the width
constant `cw` and the constant `C₀` of (10.10), not on `h`, `k`, `n` or the data. -/
theorem thm_10_6 (Ω : Set E2) (hΩo : IsOpen Ω) (hΩne : Ω.Nonempty)
    (hΩb : Bornology.IsBounded Ω) (hΩc : Convex ℝ Ω) (hΓ : SmoothBoundary Ω) :
    ∀ θ₀ : ℝ, 0 < θ₀ → ∀ cw C₀ : ℝ, ∃ C : ℝ,
    ∀ T : Finset (Fin 3 → E2), AdmissibleTriangulation Ω θ₀ cw T →
    ∀ k : ℝ, 0 < k →
    ∀ (v vh : E2 → ℝ) (f u ut utt uttt : ℝ → E2 → ℝ) (U : ℕ → E2 → ℝ),
      IsSolution101 Ω f v u ut →
      (∀ x ∈ Ω, ∀ t > 0, HasDerivAt (fun s => ut s x) (utt t x) t) →
      (∀ x ∈ Ω, ∀ t > 0, HasDerivAt (fun s => utt s x) (uttt t x) t) →
      eLpNorm (fun x => vh x - v x) 2 (domainMeasure Ω) ≤
        ENNReal.ofReal C₀ * ENNReal.ofReal (meshSize T ^ 2) * h2Norm Ω v →
      IsCN1025 Ω T k f vh U →
      ∀ n : ℕ,
        eLpNorm (fun x => U n x - u ((n : ℝ) * k) x) 2 (domainMeasure Ω) ≤
          ENNReal.ofReal C * ENNReal.ofReal (meshSize T ^ 2) *
            (h2Norm Ω v + ∫⁻ s in Set.Ioc 0 ((n : ℝ) * k), h2Norm Ω (ut s)) +
          ENNReal.ofReal C * ENNReal.ofReal (k ^ 2) *
            ∫⁻ s in Set.Ioc 0 ((n : ℝ) * k),
              (eLpNorm (uttt s) 2 (domainMeasure Ω) + lapNorm Ω (utt s)) := by
  sorry

end PDETextbook.LarssonThomee.Thm106