Second-order L2 error for interval finite elements

import PDETextbook.Base

open MeasureTheory
open scoped BigOperators ENNReal

noncomputable section

namespace PDETextbook.LarssonThomee.Thm52

/-- The coordinate of a point of the one-dimensional Euclidean space. -/
def coord (x : Euclidean 1) : ℝ := (EuclideanSpace.proj (0 : Fin 1) : Euclidean 1 →L[ℝ] ℝ) x

/-- The point of `Euclidean 1` with coordinate `t`. -/
def pt (t : ℝ) : Euclidean 1 := EuclideanSpace.single (0 : Fin 1) t

/-- `Ω = (0,1)`, realized in `Euclidean 1`. -/
def Ω : Set (Euclidean 1) := {x | 0 < coord x ∧ coord x < 1}

/-- `Ω̄ = [0,1]`. -/
def Ωbar : Set (Euclidean 1) := {x | 0 ≤ coord x ∧ coord x ≤ 1}

/-- A partition `0 = x_0 < x_1 < ⋯ < x_M = 1` of `Ω`. -/
structure Partition where
  M : ℕ
  x : Fin (M + 1) → ℝ
  strictMono : StrictMono x
  first : x 0 = 0
  last : x (Fin.last M) = 1

/-- `h = max_j h_j`, `h_j = x_j - x_{j-1}`, `j = 1,…,M`. -/
def meshSize (P : Partition) : ℝ :=
  ⨆ j : Fin P.M, (P.x j.succ - P.x j.castSucc)

/-- `v ∈ S_h`: `v ∈ C(Ω̄)`, `v` linear on each `K_j = [x_{j-1}, x_j]`, `v(0) = v(1) = 0`. -/
def InSh (P : Partition) (v : Euclidean 1 → ℝ) : Prop :=
  ContinuousOn v Ωbar ∧
  (∀ j : Fin P.M, ∃ α β : ℝ, ∀ y : Euclidean 1,
      P.x j.castSucc ≤ coord y → coord y ≤ P.x j.succ → v y = α * coord y + β) ∧
  v (pt 0) = 0 ∧ v (pt 1) = 0

/-- `a(v,w) = ∫_Ω (a v' w' + c v w) dx`, with `v'`, `w'` the (weak) derivatives. -/
def bilinForm (a c : ℝ → ℝ) (v v' w w' : Euclidean 1 → ℝ) : ℝ :=
  ∫ x in Ω, (a (coord x) * v' x * w' x + c (coord x) * v x * w x)

/-- `(f, v) = ∫_Ω f v dx`. -/
def l2Pair (f v : Euclidean 1 → ℝ) : ℝ := ∫ x in Ω, f x * v x

/-- `u ∈ H_0^1` solves (5.2): `a(u, φ) = (f, φ)` for all `φ ∈ H_0^1`. -/
def IsVariationalSolution (a c f u : Euclidean 1 → ℝ) : Prop := False

end PDETextbook.LarssonThomee.Thm52

namespace PDETextbook.LarssonThomee.Thm52'

open PDETextbook.LarssonThomee.Thm52

/-- `u ∈ H_0^1` solves (5.2): `a(u, φ) = (f, φ)` for all `φ ∈ H_0^1`.
Weak derivatives are a.e. unique, so the choice of representative `u'`, `φ'` is immaterial. -/
def IsVariationalSolution (a c : ℝ → ℝ) (f u : Euclidean 1 → ℝ) : Prop :=
  InH10Closure Ω u ∧
  ∃ u' : Euclidean 1 → ℝ, HasWeakPartial Ω u (0 : Fin 1) u' ∧ InL2 Ω u' ∧
    ∀ φ φ' : Euclidean 1 → ℝ, InH10Closure Ω φ → HasWeakPartial Ω φ (0 : Fin 1) φ' →
      bilinForm a c u u' φ φ' = l2Pair f φ

/-- `u_h ∈ S_h` solves (5.3): `a(u_h, χ) = (f, χ)` for all `χ ∈ S_h`. -/
def IsGalerkinSolution (P : Partition) (a c : ℝ → ℝ) (f uh : Euclidean 1 → ℝ) : Prop :=
  InSh P uh ∧
  ∃ uh' : Euclidean 1 → ℝ, HasWeakPartial Ω uh (0 : Fin 1) uh' ∧
    ∀ χ χ' : Euclidean 1 → ℝ, InSh P χ → HasWeakPartial Ω χ (0 : Fin 1) χ' →
      bilinForm a c uh uh' χ χ' = l2Pair f χ

/-- `‖v‖_2 = (‖v‖² + ‖v'‖² + ‖v''‖²)^{1/2}` on `Ω`, given weak derivatives `v'`, `v''`. -/
def h2Norm (v v' v'' : Euclidean 1 → ℝ) : ℝ :=
  Real.sqrt ((eLpNorm v 2 (domainMeasure Ω)).toReal ^ 2 +
    (eLpNorm v' 2 (domainMeasure Ω)).toReal ^ 2 +
    (eLpNorm v'' 2 (domainMeasure Ω)).toReal ^ 2)

/-- Theorem 5.2 (Larsson–Thomée): for `-(a u')' + c u = f` on `(0,1)`, `u(0)=u(1)=0`,
with `a, c` smooth on `[0,1]`, `a ≥ a₀ > 0`, `c ≥ 0`, `f ∈ L₂`, the piecewise linear
finite element solution satisfies `‖u_h - u‖ ≤ C h² ‖u‖₂`, with `C` independent of the
partition, `f` and `u`. -/
theorem thm_5_2 (a c : ℝ → ℝ)
    (ha : ∀ n : ℕ, ContDiffOn ℝ n a (Set.Icc (0 : ℝ) 1))
    (hc : ∀ n : ℕ, ContDiffOn ℝ n c (Set.Icc (0 : ℝ) 1))
    (ha0 : ∃ a₀ : ℝ, 0 < a₀ ∧ ∀ t ∈ Set.Icc (0 : ℝ) 1, a₀ ≤ a t)
    (hc0 : ∀ t ∈ Set.Icc (0 : ℝ) 1, 0 ≤ c t) :
    ∃ C : ℝ, ∀ (P : Partition) (f u uh u' u'' : Euclidean 1 → ℝ),
      InL2 Ω f →
      IsVariationalSolution a c f u →
      IsGalerkinSolution P a c f uh →
      HasWeakPartial Ω u (0 : Fin 1) u' → InL2 Ω u' →
      HasWeakPartial Ω u' (0 : Fin 1) u'' → InL2 Ω u'' →
      (eLpNorm (fun x => uh x - u x) 2 (domainMeasure Ω)).toReal ≤
        C * meshSize P ^ 2 * h2Norm u u' u'' := by
  sorry

end PDETextbook.LarssonThomee.Thm52'