Weak Neumann solvability with positive reaction

import PDETextbook.Base

open MeasureTheory
open scoped BigOperators ENNReal

noncomputable section

namespace PDETextbook.LarssonThomee.Thm38

/-- `Ω` is a domain (open and connected) with smooth boundary `Γ = ∂Ω`. This is realized by
a global smooth defining function `ρ` with `Ω = {ρ < 0}` and nonvanishing differential on `{ρ = 0}`. -/
def IsSmoothDomain {d : ℕ} (Ω : Set (Euclidean d)) : Prop :=
  IsOpen Ω ∧ IsConnected Ω ∧
  ∃ ρ : Euclidean d → ℝ, ContDiff ℝ (⊤ : ℕ∞) ρ ∧ Ω = {x | ρ x < 0} ∧
    ∀ x, ρ x = 0 → fderiv ℝ ρ x ≠ 0

/-- Smooth coefficient on `Ω̄` in the `C^k(Ω̄)` sense: `C^∞` on `Ω` with every derivative bounded on `Ω`. -/
def IsSmoothCoeff {d : ℕ} (Ω : Set (Euclidean d)) (f : Euclidean d → ℝ) : Prop :=
  ContDiffOn ℝ (⊤ : ℕ∞) f Ω ∧
  ∀ n : ℕ, ∃ M : ℝ, ∀ x ∈ Ω, ‖iteratedFDerivWithin ℝ n f Ω x‖ ≤ M

/-- `a(u,v) = ∫_Ω (a ∇u·∇v + c u v) dx`, i.e. (3.25) with `b = 0`, as in (3.31).
Here `G` and `H` are the weak gradients of `u` and `v`. -/
def neumannForm {d : ℕ} (Ω : Set (Euclidean d)) (a c : Euclidean d → ℝ)
    (u : Euclidean d → ℝ) (G : Fin d → Euclidean d → ℝ)
    (v : Euclidean d → ℝ) (H : Fin d → Euclidean d → ℝ) : ℝ :=
  ∫ x in Ω, (a x * ∑ i, G i x * H i x + c x * (u x * v x))

/-- `L(v) = (f, v) = ∫_Ω f v dx`. -/
def loadFunctional {d : ℕ} (Ω : Set (Euclidean d)) (f v : Euclidean d → ℝ) : ℝ :=
  ∫ x in Ω, f x * v x

/-- Weak solution of the Neumann problem (3.30): `u ∈ H¹` (with weak gradient `G`) and
`a(u,v) = L(v)` for all `v ∈ H¹`, i.e. (3.33). -/
def IsNeumannWeakSolution {d : ℕ} (Ω : Set (Euclidean d)) (a c f : Euclidean d → ℝ)
    (u : Euclidean d → ℝ) (G : Fin d → Euclidean d → ℝ) : Prop :=
  HasH1Gradient Ω u G ∧
  ∀ (v : Euclidean d → ℝ) (H : Fin d → Euclidean d → ℝ), HasH1Gradient Ω v H →
    neumannForm Ω a c u G v H = loadFunctional Ω f v

/-- `‖u‖₁ = (∑_{|α|≤1} ‖D^α u‖²)^{1/2}`, with `G` the weak gradient of `u`. -/
def h1Norm {d : ℕ} (Ω : Set (Euclidean d)) (u : Euclidean d → ℝ)
    (G : Fin d → Euclidean d → ℝ) : ℝ :=
  Real.sqrt ((∫ x in Ω, u x ^ 2) + ∑ i, ∫ x in Ω, G i x ^ 2)

/-- `‖f‖ = (∫_Ω f² dx)^{1/2}`. -/
def l2Norm {d : ℕ} (Ω : Set (Euclidean d)) (f : Euclidean d → ℝ) : ℝ :=
  Real.sqrt (∫ x in Ω, f x ^ 2)

/-- Theorem 3.8 (Larsson–Thomée). Hypotheses (3.2) are `a ≥ a₀ > 0` on `Ω` with smooth
coefficients. For the Neumann problem, the source additionally requires `c ≥ c₀ > 0` on `Ω`.
For `f ∈ L₂` there is a unique (up to a.e. equality on `Ω`) `u ∈ H¹` satisfying (3.33).
Moreover `‖u‖₁ ≤ C‖f‖`, with `C` independent of `f`. -/
theorem neumann_weak_solution_exists_unique {d : ℕ} (Ω : Set (Euclidean d))
    (hΩ : IsSmoothDomain Ω)
    (a c : Euclidean d → ℝ) (ha : IsSmoothCoeff Ω a) (hc : IsSmoothCoeff Ω c)
    (a₀ c₀ : ℝ) (ha₀ : 0 < a₀) (hc₀ : 0 < c₀)
    (ha_lower : ∀ x ∈ Ω, a₀ ≤ a x) (hc_lower : ∀ x ∈ Ω, c₀ ≤ c x) :
    ∃ C : ℝ, ∀ f : Euclidean d → ℝ, InL2 Ω f →
      ∃ (u : Euclidean d → ℝ) (G : Fin d → Euclidean d → ℝ),
        IsNeumannWeakSolution Ω a c f u G ∧
        (∀ (u' : Euclidean d → ℝ) (G' : Fin d → Euclidean d → ℝ),
          IsNeumannWeakSolution Ω a c f u' G' → u' =ᵐ[domainMeasure Ω] u) ∧
        h1Norm Ω u G ≤ C * l2Norm Ω f := by
  sorry

end PDETextbook.LarssonThomee.Thm38