Second-order approximation of the first interval eigenvalue

import PDETextbook.Base

open MeasureTheory
open scoped BigOperators ENNReal

noncomputable section

namespace PDETextbook.LarssonThomee.Thm66

/-- Classical eigenpair of (6.26): `-(a φ')' + c φ = lam φ` in `Ω = (0,1)`, with
`φ(0) = φ(1) = 0`. Here `φ ∈ C²(Ω) ∩ C(Ω̄)` and `φ` is not identically zero. -/
def IsContEigenpair (a c : ℝ → ℝ) (lam : ℝ) (φ : ℝ → ℝ) : Prop :=
  ContinuousOn φ (Set.Icc 0 1) ∧ ContDiffOn ℝ 2 φ (Set.Ioo 0 1) ∧
  φ 0 = 0 ∧ φ 1 = 0 ∧ (∃ x ∈ Set.Ioo (0 : ℝ) 1, φ x ≠ 0) ∧
  ∀ x ∈ Set.Ioo (0 : ℝ) 1,
    -(deriv (fun y => a y * deriv φ y) x) + c x * φ x = lam * φ x

/-- The difference operator of (6.27) at the mesh point `x_j = j h`, `h = 1/M`:
`A_h U_j = -∂̄(a_{j+1/2} ∂U_j) + c_j U_j`, with `∂U_j = (U_{j+1}-U_j)/h`,
`∂̄V_j = (V_j - V_{j-1})/h`, `a_{j+1/2} = a(x_j + h/2)`, `c_j = c(x_j)`. -/
def meshOp (a c : ℝ → ℝ) (M : ℕ) (U : ℕ → ℝ) (j : ℕ) : ℝ :=
  let h : ℝ := 1 / (M : ℝ);
  let x : ℝ := (j : ℝ) * h;
  -(a (x + h / 2) * ((U (j + 1) - U j) / h) - a (x - h / 2) * ((U j - U (j - 1)) / h)) / h
    + c x * U j

/-- `Λ` is an eigenvalue of the finite-dimensional problem (6.27): there is a mesh function
`U` with `U_0 = U_M = 0`, nonzero at some interior mesh point, satisfying
`A_h U_j = Λ U_j` for `j = 1, …, M-1` (equivalently `A Ū = Λ Ū` with `Ū ≠ 0`). -/
def IsDiscreteEigenvalue (a c : ℝ → ℝ) (M : ℕ) (Λ : ℝ) : Prop :=
  ∃ U : ℕ → ℝ, U 0 = 0 ∧ U M = 0 ∧ (∃ j, 1 ≤ j ∧ j ≤ M - 1 ∧ U j ≠ 0) ∧
    ∀ j, 1 ≤ j → j ≤ M - 1 → meshOp a c M U j = Λ * U j

/-- Theorem 6.6: if `Λ₁` and `λ₁` are the smallest eigenvalues of (6.27) and (6.26),
then `|Λ₁ - λ₁| ≤ C h²`, where `h = 1/M` and `C` depends only on the data `a, c`
(not on `h`). -/
theorem thm_6_6 (a c : ℝ → ℝ) (a₀ : ℝ) (ha₀ : 0 < a₀)
    (ha : ContDiffOn ℝ (⊤ : ℕ∞) a (Set.Icc 0 1))
    (hc : ContDiffOn ℝ (⊤ : ℕ∞) c (Set.Icc 0 1))
    (haa₀ : ∀ x ∈ Set.Icc (0 : ℝ) 1, a₀ ≤ a x)
    (hc0 : ∀ x ∈ Set.Icc (0 : ℝ) 1, 0 ≤ c x) :
    ∃ C : ℝ, ∀ M : ℕ, 2 ≤ M → ∀ Λ₁ lam₁ : ℝ,
      IsLeast {Λ | IsDiscreteEigenvalue a c M Λ} Λ₁ →
      IsLeast {lam | ∃ φ, IsContEigenpair a c lam φ} lam₁ →
      |Λ₁ - lam₁| ≤ C * (1 / (M : ℝ)) ^ 2 := by
  sorry

end PDETextbook.LarssonThomee.Thm66