import PDETextbook.Base

open MeasureTheory
open scoped BigOperators ENNReal

noncomputable section

namespace PDETextbook.LarssonThomee.Thm96

open Set

/-- Mesh width `h = 1/M` of the grid `x_j = j h`, `j = 0,…,M`, on `Ω̄ = [0,1]`. -/
def meshWidth (M : ℕ) : ℝ := 1 / (M : ℝ)

/-- Discrete maximum norm `‖W‖_{∞,h} = max_{0 ≤ j ≤ M} |W_j|`. -/
def discreteMaxNorm (M : ℕ) (W : ℕ → ℝ) : ℝ :=
  ⨆ j : Fin (M + 1), |W (j : ℕ)|

/-- `|w|_{C^k(Ω̄)} = max_{|α| = k} ‖D^α w‖_{C(Ω̄)}`, here `sup_{x ∈ [0,1]} |w^{(k)}(x)|`. -/
def cSeminorm (k : ℕ) (w : ℝ → ℝ) : ℝ :=
  ⨆ x : Icc (0 : ℝ) 1, |iteratedDerivWithin k w (Icc 0 1) (x : ℝ)|

/-- Classical solution of (9.19): `u_t = u_xx` in `(0,1) × (0,∞)`,
`u(0,t) = u(1,t) = 0` for `t > 0`, `u(·,0) = v` in `(0,1)`,
`u` continuous on `[0,1] × [0,∞)`. -/
def IsHeatDirichletSolution (u : ℝ → ℝ → ℝ) (v : ℝ → ℝ) : Prop :=
  ContinuousOn (fun p : ℝ × ℝ => u p.1 p.2) (Icc (0 : ℝ) 1 ×ˢ Ici (0 : ℝ)) ∧
  (∀ x ∈ Ioo (0 : ℝ) 1, ∀ t : ℝ, 0 < t →
    DifferentiableAt ℝ (fun s => u x s) t ∧
    ContDiffAt ℝ 2 (fun y => u y t) x ∧
    deriv (fun s => u x s) t = iteratedDeriv 2 (fun y => u y t) x) ∧
  (∀ t : ℝ, 0 < t → u 0 t = 0 ∧ u 1 t = 0) ∧
  (∀ x ∈ Ioo (0 : ℝ) 1, u x 0 = v x)

/-- Backward Euler scheme (9.21): `∂̄_t U_j^{n+1} = ∂_x ∂̄_x U_j^{n+1}` (j = 1..M-1, n ≥ 0),
`U_0^{n+1} = U_M^{n+1} = 0`, `U_j^0 = v(x_j)` (j = 0..M), with
`∂̄_t W^{n+1} = (W^{n+1}-W^n)/k`, `∂_x W_j = (W_{j+1}-W_j)/h`, `∂̄_x W_j = (W_j-W_{j-1})/h`. -/
def IsBackwardEulerSolution (M : ℕ) (k : ℝ) (v : ℝ → ℝ) (U : ℕ → ℕ → ℝ) : Prop :=
  (∀ n j : ℕ, 1 ≤ j → j ≤ M - 1 →
    (U (n + 1) j - U n j) / k =
      ((U (n + 1) (j + 1) - U (n + 1) j) / meshWidth M -
        (U (n + 1) j - U (n + 1) (j - 1)) / meshWidth M) / meshWidth M) ∧
  (∀ n : ℕ, U (n + 1) 0 = 0 ∧ U (n + 1) M = 0) ∧
  (∀ j : ℕ, j ≤ M → U 0 j = v ((j : ℝ) * meshWidth M))

/-- Theorem 9.6: `‖U^n - u^n‖_{∞,h} ≤ C t_n (h^2 + k) max_{t ≤ t_n} |u(·,t)|_{C^4}`, `t_n ≥ 0`.
The maximum is encoded via an arbitrary upper bound `K` of `|u(·,s)|_{C^4}`, `0 ≤ s ≤ t_n`
(equivalent when finite; vacuous, as in the source, when infinite). -/
theorem thm_9_6 :
    ∃ C : ℝ, ∀ (M : ℕ), 1 ≤ M → ∀ (k : ℝ), 0 < k →
      ∀ (v : ℝ → ℝ) (u : ℝ → ℝ → ℝ) (U : ℕ → ℕ → ℝ),
        IsHeatDirichletSolution u v →
        ContinuousOn v (Icc (0 : ℝ) 1) →
        IsBackwardEulerSolution M k v U →
        ∀ n : ℕ,
          (∀ s ∈ Icc (0 : ℝ) ((n : ℝ) * k),
            ContDiffOn ℝ 4 (fun y => u y s) (Icc (0 : ℝ) 1)) →
          ∀ K : ℝ,
            (∀ s ∈ Icc (0 : ℝ) ((n : ℝ) * k), cSeminorm 4 (fun y => u y s) ≤ K) →
            discreteMaxNorm M
                (fun j => U n j - u ((j : ℝ) * meshWidth M) ((n : ℝ) * k)) ≤
              C * ((n : ℝ) * k) * (meshWidth M ^ 2 + k) * K := by
  sorry

end PDETextbook.LarssonThomee.Thm96