import PDETextbook.Base

open MeasureTheory
open scoped BigOperators ENNReal

noncomputable section

namespace PDETextbook.LarssonThomee.Thm13_3

/-- Partition `0 = x_0 < x_1 < ... < x_M = 1` of `Ω = (0,1)`. -/
structure Mesh where
  M : ℕ
  x : Fin (M + 1) → ℝ
  x_zero : x 0 = 0
  x_last : x (Fin.last M) = 1
  strictMono : StrictMono x

/-- `h = max_j h_j` with `h_j = x_j - x_{j-1}`. -/
def Mesh.h (m : Mesh) : ℝ := ⨆ j : Fin m.M, (m.x j.succ - m.x j.castSucc)

/-- Membership in `S_h^-` of (13.7): continuous on `[0,1]`, linear on each
`K_j = [x_{j-1}, x_j]`, and `χ(0) = 0`. -/
def InShMinus (m : Mesh) (χ : ℝ → ℝ) : Prop :=
  ContinuousOn χ (Set.Icc 0 1) ∧
  (∀ j : Fin m.M, ∃ a b : ℝ,
    ∀ y ∈ Set.Icc (m.x j.castSucc) (m.x j.succ), χ y = a * y + b) ∧
  χ 0 = 0

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

/-- Transport of a function of `x ∈ ℝ` to `Euclidean 1`. -/
def lift (w : ℝ → ℝ) : Euclidean 1 → ℝ := fun p => w (p 0)

/-- `‖w‖ = ‖w‖_{L₂(Ω)}` (value `∞` if not square integrable). -/
def l2Norm (w : ℝ → ℝ) : ℝ≥0∞ := eLpNorm (lift w) 2 (domainMeasure Ω1)

/-- `‖w‖_1 = (‖w‖² + ‖w'‖²)^{1/2}` with weak derivative `w'`; `∞` if none exists. -/
def h1Norm (w : ℝ → ℝ) : ℝ≥0∞ :=
  ⨅ (g : ℝ → ℝ) (_ : HasWeakPartial Ω1 (lift w) 0 (lift g)),
    (l2Norm w ^ 2 + l2Norm g ^ 2) ^ (1 / 2 : ℝ)

/-- `‖w‖_2 = (‖w‖² + ‖w'‖² + ‖w''‖²)^{1/2}` with weak derivatives; `∞` if none exist. -/
def h2Norm (w : ℝ → ℝ) : ℝ≥0∞ :=
  ⨅ (g : ℝ → ℝ) (g2 : ℝ → ℝ)
    (_ : HasWeakPartial Ω1 (lift w) 0 (lift g) ∧ HasWeakPartial Ω1 (lift g) 0 (lift g2)),
    (l2Norm w ^ 2 + l2Norm g ^ 2 + l2Norm g2 ^ 2) ^ (1 / 2 : ℝ)

/-- `u` (written `u x t`) solves (13.6): `u_t + u_x = f` in `Ω` for `t > 0`,
`u(0,t) = 0` for `t > 0`, `u(·,0) = v` in `Ω`; classical pointwise reading,
continuous on `Ω̄ × [0,∞)`. -/
def IsSolution136 (f : ℝ → ℝ → ℝ) (v : ℝ → ℝ) (u : ℝ → ℝ → ℝ) : Prop :=
  ContinuousOn (Function.uncurry u) (Set.Icc 0 1 ×ˢ Set.Ici 0) ∧
  (∀ x ∈ Set.Ioo (0 : ℝ) 1, ∀ t : ℝ, 0 < t →
    DifferentiableAt ℝ (fun r => u x r) t ∧ DifferentiableAt ℝ (fun y => u y t) x ∧
    deriv (fun r => u x r) t + deriv (fun y => u y t) x = f x t) ∧
  (∀ t : ℝ, 0 < t → u 0 t = 0) ∧
  (∀ x ∈ Set.Ioo (0 : ℝ) 1, u x 0 = v x)

/-- `uh` (written `uh t x`) solves (13.8): `u_h(t) ∈ S_h^-` for `t ≥ 0`,
`(u_{h,t} + u_{h,x}, χ) = (f, χ)` for all `χ ∈ S_h^-`, `t > 0`, and `u_h(0) = v_h`. -/
def IsSolution138 (m : Mesh) (f : ℝ → ℝ → ℝ) (vh : ℝ → ℝ) (uh : ℝ → ℝ → ℝ) : Prop :=
  (∀ t : ℝ, 0 ≤ t → InShMinus m (uh t)) ∧
  (∀ x ∈ Set.Icc (0 : ℝ) 1, ContinuousWithinAt (fun s => uh s x) (Set.Ici 0) 0) ∧
  (∀ t : ℝ, 0 < t → ∀ x ∈ Set.Icc (0 : ℝ) 1, DifferentiableAt ℝ (fun s => uh s x) t) ∧
  (∀ t : ℝ, 0 < t → ∀ χ : ℝ → ℝ, InShMinus m χ →
    (∫ x in Set.Ioo (0 : ℝ) 1, (deriv (fun s => uh s x) t + deriv (uh t) x) * χ x) =
      ∫ x in Set.Ioo (0 : ℝ) 1, f x t * χ x) ∧
  (∀ x ∈ Set.Icc (0 : ℝ) 1, uh 0 x = vh x)

/-- Theorem 13.3: with `v_h` suitably chosen,
`‖u_h(t) - u(t)‖ ≤ C h (‖v‖_1 + ∫_0^t (‖u‖_2 + ‖u_t‖_1) ds)` for `t ≥ 0`. -/
theorem thm_13_3 :
    ∃ C : ℝ, 0 < C ∧
      ∀ (m : Mesh) (f : ℝ → ℝ → ℝ) (v : ℝ → ℝ) (u : ℝ → ℝ → ℝ),
        IsSolution136 f v u →
        ∃ vh : ℝ → ℝ, InShMinus m vh ∧
          ∀ uh : ℝ → ℝ → ℝ, IsSolution138 m f vh uh →
            ∀ t : ℝ, 0 ≤ t →
              l2Norm (fun x => uh t x - u x t) ≤
                ENNReal.ofReal (C * m.h) *
                  (h1Norm v + ∫⁻ s in Set.Ioo 0 t,
                    (h2Norm (fun x => u x s) +
                      h1Norm (fun x => deriv (fun r => u x r) s))) := by
  sorry

end PDETextbook.LarssonThomee.Thm13_3