Stability for a boundary-adapted difference scheme

import PDETextbook.Base

open MeasureTheory
open scoped BigOperators ENNReal

noncomputable section

namespace PDETextbook.LarssonThomee.Lem45

abbrev E2 := Euclidean 2

/-- Mesh-width `h = 1/M`. -/
def meshWidth (M : ℕ) : ℝ := 1 / (M : ℝ)

/-- Mesh-point `x_j = j h`, `j ∈ ℤ²`. -/
def meshPt (M : ℕ) (j : ℤ × ℤ) : E2 :=
  (EuclideanSpace.equiv (Fin 2) ℝ).symm
    ![(j.1 : ℝ) * meshWidth M, (j.2 : ℝ) * meshWidth M]

def e1 : ℤ × ℤ := (1, 0)
def e2 : ℤ × ℤ := (0, 1)

/-- `x_i` is one of the four neighbors `x_{j ± e_l}` of `x_j`. -/
def IsMeshNeighbor (i j : ℤ × ℤ) : Prop :=
  i = j + e1 ∨ i = j - e1 ∨ i = j + e2 ∨ i = j - e2

/-- Smooth boundary: near each boundary point `Ω = {φ < 0}` for a `C^∞` `φ` with `dφ(p) ≠ 0`. -/
def HasSmoothBoundary (Ω : Set E2) : Prop :=
  ∀ p ∈ frontier Ω, ∃ V ∈ nhds p, ∃ φ : E2 → ℝ,
    ContDiff ℝ (⊤ : ℕ∞) φ ∧ fderiv ℝ φ p ≠ 0 ∧ Ω ∩ V = {x | x ∈ V ∧ φ x < 0}

/-- `Ω_h`: interior mesh-points all four of whose neighbors lie in `Ω̄`. -/
def OmegaH (Ω : Set E2) (M : ℕ) : Set (ℤ × ℤ) :=
  {j | meshPt M j ∈ Ω ∧ ∀ k, IsMeshNeighbor k j → meshPt M k ∈ closure Ω}

/-- `ω_h`: mesh-points in `Ω` not in `Ω_h`. -/
def omegaH (Ω : Set E2) (M : ℕ) : Set (ℤ × ℤ) :=
  {j | meshPt M j ∈ Ω ∧ j ∉ OmegaH Ω M}

/-- Selection preceding (4.14): neighbor `x_i = x_{nb j} ∈ Ω_h ∪ ω_h` of `x_j ∈ ω_h` whose
grid line through `x_j` cuts `Γ` at the non-mesh-point `x̄_j = xbar j`, with `x_j` between
`x_i` and `x̄_j` and `|x_j - x̄_j| < h` (so `α_j ≤ 1/2` as asserted in (4.14)). -/
def IsAdmissibleSelection (Ω : Set E2) (M : ℕ) (nb : ℤ × ℤ → ℤ × ℤ)
    (xbar : ℤ × ℤ → E2) : Prop :=
  ∀ j ∈ omegaH Ω M,
    IsMeshNeighbor (nb j) j ∧ nb j ∈ OmegaH Ω M ∪ omegaH Ω M ∧
    xbar j ∈ frontier Ω ∧
    (∃ s : ℝ, 0 < s ∧ s < 1 ∧
      xbar j = meshPt M j + s • (meshPt M j - meshPt M (nb j))) ∧
    ∀ k : ℤ × ℤ, xbar j ≠ meshPt M k

/-- `Γ_h`: points of `Γ` that are neighbors of points in `Ω_h` or points `x̄_j`, `x_j ∈ ω_h`. -/
def GammaH (Ω : Set E2) (M : ℕ) (xbar : ℤ × ℤ → E2) : Set E2 :=
  {p | p ∈ frontier Ω ∧
    ((∃ j ∈ OmegaH Ω M, ∃ k, IsMeshNeighbor k j ∧ p = meshPt M k) ∨
      ∃ j ∈ omegaH Ω M, p = xbar j)}

/-- `Δ_h U_j = ∂₁∂̄₁U_j + ∂₂∂̄₂U_j`, (4.10)-(4.11). -/
def discreteLaplacian (M : ℕ) (U : E2 → ℝ) (j : ℤ × ℤ) : ℝ :=
  (U (meshPt M (j + e1)) - 2 * U (meshPt M j) + U (meshPt M (j - e1))) / (meshWidth M) ^ 2 +
  (U (meshPt M (j + e2)) - 2 * U (meshPt M j) + U (meshPt M (j - e2))) / (meshWidth M) ^ 2

/-- `α_j = |x_j - x̄_j| / (h + |x_j - x̄_j|)`, (4.14). -/
def alphaCoeff (M : ℕ) (xbar : ℤ × ℤ → E2) (j : ℤ × ℤ) : ℝ :=
  dist (meshPt M j) (xbar j) / (meshWidth M + dist (meshPt M j) (xbar j))

/-- `ℓ_h U_j = U_j - α_j U_i - (1-α_j) U(x̄_j)`, (4.14). -/
def ellH (M : ℕ) (nb : ℤ × ℤ → ℤ × ℤ) (xbar : ℤ × ℤ → E2) (U : E2 → ℝ)
    (j : ℤ × ℤ) : ℝ :=
  U (meshPt M j) - alphaCoeff M xbar j * U (meshPt M (nb j))
    - (1 - alphaCoeff M xbar j) * U (xbar j)

/-- (4.7) `|V|_S = max_{x_j ∈ S} |V_j|` over mesh indices (finite `S` for bounded `Ω`). -/
def meshMaxNorm (S : Set (ℤ × ℤ)) (V : ℤ × ℤ → ℝ) : ℝ :=
  sSup ((fun j => |V j|) '' S)

/-- (4.7) over a point set such as `Γ_h` (which contains the non-mesh points `x̄_j`). -/
def pointMaxNorm (S : Set E2) (U : E2 → ℝ) : ℝ :=
  sSup ((fun p => |U p|) '' S)

/-- Lemma 4.5: `|U|_{Ω_h ∪ ω_h} ≤ 2(|U|_{Γ_h} + |ℓ_h U|_{ω_h} + C|Δ_h U|_{Ω_h})` for any
mesh-function `U`, in the setting of a convex plane domain with smooth boundary. -/
theorem lem_4_5 (Ω : Set E2) (hopen : IsOpen Ω) (hconn : IsConnected Ω)
    (hconv : Convex ℝ Ω) (hbdd : Bornology.IsBounded Ω) (hsmooth : HasSmoothBoundary Ω) :
    ∃ C : ℝ, ∀ M : ℕ, 0 < M →
      ∀ (nb : ℤ × ℤ → ℤ × ℤ) (xbar : ℤ × ℤ → E2), IsAdmissibleSelection Ω M nb xbar →
      ∀ U : E2 → ℝ,
        meshMaxNorm (OmegaH Ω M ∪ omegaH Ω M) (fun j => U (meshPt M j)) ≤
          2 * (pointMaxNorm (GammaH Ω M xbar) U +
            meshMaxNorm (omegaH Ω M) (ellH M nb xbar U) +
            C * meshMaxNorm (OmegaH Ω M) (discreteLaplacian M U)) := by
  sorry

end PDETextbook.LarssonThomee.Lem45