Second-order convergence of the five-point Poisson scheme

import PDETextbook.Base

open MeasureTheory
open scoped BigOperators ENNReal

noncomputable section

namespace PDETextbook.LarssonThomee.Thm4_2

/-- The square `Ω = (0,1) × (0,1)` of Section 4.2 (standing assumption under (4.9)). -/
def unitSquare : Set (Euclidean 2) :=
  {x | ∀ l : Fin 2, 0 < x l ∧ x l < 1}

/-- Mesh point `x_j = j h` with `h = 1/M`, `j = (j₁, j₂) ∈ ℤ²`. -/
def meshPoint (M : ℕ) (j : ℤ × ℤ) : Euclidean 2 :=
  EuclideanSpace.single (0 : Fin 2) ((j.1 : ℝ) / M) +
    EuclideanSpace.single (1 : Fin 2) ((j.2 : ℝ) / M)

/-- Unit index shifts `e₁ = (1,0)`, `e₂ = (0,1)`. -/
def unitShift (l : Fin 2) : ℤ × ℤ := if l = 0 then (1, 0) else (0, 1)

/-- `∂_l ∂̄_l U_j = (U_{j+e_l} - 2U_j + U_{j-e_l}) / h²`, as in (4.10), with `h = 1/M`. -/
def secondDiff (M : ℕ) (U : ℤ × ℤ → ℝ) (l : Fin 2) (j : ℤ × ℤ) : ℝ :=
  (U (j + unitShift l) - 2 * U j + U (j - unitShift l)) / (1 / (M : ℝ)) ^ 2

/-- `Δ_h U_j = ∂₁∂̄₁ U_j + ∂₂∂̄₂ U_j`, as in (4.11). -/
def discreteLaplacian (M : ℕ) (U : ℤ × ℤ → ℝ) (j : ℤ × ℤ) : ℝ :=
  secondDiff M U 0 j + secondDiff M U 1 j

/-- `U` solves the finite difference problem (4.11): `-Δ_h U_j = f_j` for `x_j ∈ Ω`,
`U_j = 0` for `x_j ∈ Γ`. -/
def IsDiscreteSolution (M : ℕ) (f : Euclidean 2 → ℝ) (U : ℤ × ℤ → ℝ) : Prop :=
  (∀ j : ℤ × ℤ, meshPoint M j ∈ unitSquare →
      -discreteLaplacian M U j = f (meshPoint M j)) ∧
  (∀ j : ℤ × ℤ, meshPoint M j ∈ frontier unitSquare → U j = 0)

/-- Classical Laplacian `Δu = ∂₁²u + ∂₂²u`. -/
def laplacian (u : Euclidean 2 → ℝ) (x : Euclidean 2) : ℝ :=
  ∑ i : Fin 2, partialDeriv (partialDeriv u i) i x

/-- `u` is a (classical) solution of the Dirichlet problem (4.9) on the square:
`-Δu = f` in `Ω`, `u = 0` on `Γ`. -/
def IsPoissonSolution (f u : Euclidean 2 → ℝ) : Prop :=
  ContDiffOn ℝ 2 u unitSquare ∧ ContinuousOn u (closure unitSquare) ∧
  (∀ x ∈ unitSquare, -laplacian u x = f x) ∧
  (∀ x ∈ frontier unitSquare, u x = 0)

/-- Discrete maximum norm (4.7): `|V|_S = max_{x_j ∈ S} |V_j|`. -/
def meshMaxNorm (M : ℕ) (S : Set (Euclidean 2)) (V : ℤ × ℤ → ℝ) : ℝ :=
  sSup ((fun j => |V j|) '' {j : ℤ × ℤ | meshPoint M j ∈ S})

/-- `C⁴` seminorm `|u|_{C⁴}`: supremum over `Ω` of the size of the fourth-order
derivatives (realized via the norm of the fourth Fréchet derivative, equivalent to
`max_{|α|=4} sup |D^α u|` up to a fixed dimensional factor). -/
def c4Seminorm (u : Euclidean 2 → ℝ) : ℝ :=
  sSup ((fun x => ‖iteratedFDeriv ℝ 4 u x‖) '' unitSquare)

/-- Theorem 4.2 (Larsson–Thomée): for the solutions `U` of (4.11) and `u` of (4.9) on the
unit square, `|U - u|_Ω ≤ C h² |u|_{C⁴}` with `h = 1/M` and `C` independent of `h` and `u`.
The finiteness of `|u|_{C⁴}` implicit in the statement is expressed as `u ∈ C⁴(Ω̄)`. -/
theorem thm_4_2 :
    ∃ C : ℝ, ∀ (M : ℕ), 0 < M →
      ∀ (f u : Euclidean 2 → ℝ) (U : ℤ × ℤ → ℝ),
        IsPoissonSolution f u →
        ContDiffOn ℝ 4 u (closure unitSquare) →
        IsDiscreteSolution M f U →
        meshMaxNorm M unitSquare (fun j => U j - u (meshPoint M j)) ≤
          C * (1 / (M : ℝ)) ^ 2 * c4Seminorm u := by
  sorry

end PDETextbook.LarssonThomee.Thm4_2