Simple root of second kind Chebyshev

Stream: mathlib4

[Nick_adfor]

I think this theorem is a good one for mathlib. Now I've finished the construction, but some details need to be filled.

Theorem: Prove that all roots of the second kind Chebyshev polynomials (Un(x)) are real and simple.
Proof:

The second kind Chebyshev polynomial can be expressed as:
Un(x) = sin((n+1)θ) / sinθ, where x = cosθ
Its roots are given by:
x_k = cos(kπ / (n+1)), for k = 1, 2, ..., n
Since when k = 1, 2, ..., n, we have 0 < kπ / (n+1) < π, all roots are real and distinct.

To prove that the roots are simple, we need to verify that for all k, Un'(x_k) ≠ 0. Using the chain rule to differentiate:
dUn/dx = (dUn/dθ) · (dθ/dx) = d/dθ [sin((n+1)θ) / sinθ] · dθ/dx
Since x = cosθ, we have dx/dθ = -sinθ, thus dθ/dx = -1/sinθ.
At θ = θ_k = kπ / (n+1), sin((n+1)θ_k) = 0. Calculating the derivative:
[d/dθ (sin((n+1)θ) / sinθ)]_θ_k = [(n+1)cos((n+1)θ_k)sinθ_k - sin((n+1)θ_k)cosθ_k] / sin²θ_k
= (n+1)cos(kπ)sinθ_k / sin²θ_k
= (n+1)(-1)^k / sinθ_k
Therefore:
Un'(x_k) = [(n+1)(-1)^k / sinθ_k] · [-1/sinθ_k] = -(n+1)(-1)^k / sin²θ_k
Since θ_k = kπ / (n+1) ∈ (0,π), we have sinθ_k ≠ 0, thus Un'(x_k) ≠ 0. Therefore, all roots are real and simple.

import Mathlib

open Polynomial Real
open Polynomial
open Chebyshev

set_option maxHeartbeats 1000000

/-- Recursion relation for Chebyshev polynomials of the second kind.
For any natural number n ≥ 2, Uₙ(x) = 2xUₙ₋₁(x) - Uₙ₋₂(x). -/
lemma Chebyshev_U_recursion :
  ∀ n : ℕ, n ≥ 2 → Chebyshev.U ℝ n = 2 * X * Chebyshev.U ℝ (n-1) - Chebyshev.U ℝ (n-2) := by
  -- Introduce n and the hypothesis n ≥ 2
  intro n hn
  -- Apply the standard recursion relation for Chebyshev polynomials
  exact Chebyshev.U_eq ℝ ↑n

/-- Evaluation of Chebyshev polynomials at x = 1.
For any natural number n, Uₙ(1) = n + 1. -/
lemma Chebyshev_U_eval_one (n : ℕ) : (Chebyshev.U ℝ n).eval (1 : ℝ) = (n + 1 : ℝ) := by
  -- Use induction on n
  induction' n with n IH
  · -- Base case: n = 0
    simp
  · -- Inductive step
    cases' n with n
    · -- Case n = 1
      simp
      ring_nf
    · -- General case
      simp

/-- Non-vanishing at endpoints for scaled Chebyshev polynomials.
For any natural number n, Uₙ(1) ≠ 0 and Uₙ(-1) ≠ 0. -/
lemma endpoints_nonzero (n : ℕ) :
    ((Chebyshev.U ℝ n).comp (Polynomial.C (1/2) * X)).eval 2 ≠ 0 ∧
    ((Chebyshev.U ℝ n).comp (Polynomial.C (1/2) * X)).eval (-2) ≠ 0 := by
  -- Split the conjunction into two goals
  constructor
  · -- First goal: evaluation at 2 is nonzero
    simp
    -- Use linear arithmetic with the fact that n+1 > 0
    linarith [show (0 : ℝ) < n + 1 from by exact_mod_cast Nat.succ_pos n]
  · -- Second goal: evaluation at -2 is nonzero
    simp
    -- n+1 is nonzero as a real number
    exact Nat.cast_add_one_ne_zero n

/-- Evaluation of Chebyshev polynomials at x = -1.
For any natural number n, Uₙ(-1) = (-1)ⁿ(n + 1). -/
lemma Chebyshev.U_eval_neg_one (n : ℕ) :
    (Chebyshev.U ℝ n).eval (-1) = (-1 : ℝ) ^ n * (n + 1 : ℝ) := by
  -- Use strong induction on n
  induction' n using Nat.strong_induction_on with n IH
  -- Case analysis on n
  rcases n with (rfl|rfl|m)
  · -- Case n = 0
    simp
  · -- Case n = 1
    simp
    ring_nf
  · -- Case n = m + 2 where m is a natural number
    -- Rewrite m+2 as (m+1)+1 for clarity
    have : m + 2 = (m + 1) + 1 := by omega
    rw [<-this]
    -- Inductive hypotheses for m
    have hm1 : m < m + 2 := by omega
    -- Inductive hypotheses for m+1
    have hm2 : m + 1 < m + 2 := by omega
    -- Get the recursion relation for Chebyshev polynomials
    have rec_eq := Chebyshev.U_add_two ℝ (m)
    simp only at rec_eq ⊢
    -- Rewrite the evaluation using the recursion relation
    have eval_rec_eq : (Chebyshev.U ℝ (m + 2)).eval (-1) =
                      (2 * X * Chebyshev.U ℝ (m + 1) - Chebyshev.U ℝ m).eval (-1) := by
      rw [rec_eq]
    -- Type casting lemmas
    have cast_eq1 : (U ℝ (m + 2) : Polynomial ℝ) = U ℝ (↑(m + 2) : ℕ) := by simp
    -- Type casting lemmas
    have cast_eq2 : (U ℝ (m + 1) : Polynomial ℝ) = U ℝ (↑(m + 1) : ℕ) := by simp
    -- Type casting lemmas
    have cast_eq3 : (U ℝ m : Polynomial ℝ) = U ℝ (↑m : ℕ) := by simp
    -- Begin the main calculation
    rw [← cast_eq1]
    rw [rec_eq]
    simp only [eval_sub, eval_mul, eval_ofNat, eval_X]
    -- Apply inductive hypotheses
    rw [IH m hm1]
    rw [cast_eq2]
    rw [IH (m + 1) hm2]
    simp
    ring_nf

/-- Parity property of Chebyshev polynomials.
For any natural number n and real number x, Uₙ(-x) = (-1)ⁿUₙ(x). -/
lemma Chebyshev.U_neg (n : ℕ) (x : ℝ) :
    (U ℝ n).eval (-x) = (-1 : ℝ)^n * (U ℝ n).eval x := by
  -- Use strong induction on n
  induction' n using Nat.strong_induction_on with n ih
  -- Case analysis on n
  rcases n with (rfl|rfl|m)
  · -- Case n = 0
    simp
  · -- Case n = 1
    simp
  · -- Case n = m + 2
    -- Express evaluation using recursion relation
    have eval_rec : (U ℝ (m + 2)).eval (-x) =
                    2 * (-x) * (U ℝ (m + 1)).eval (-x) - (U ℝ m).eval (-x) := by
                      rw [ih m (Nat.le.step Nat.le.refl)]
                      simp
                      -- Calculation
                      have h2 := ih m (by omega)
                      rw [h2]
    -- Main calculation
    calc
      (U ℝ (m + 2)).eval (-x)
          = 2 * (-x) * (U ℝ (m + 1)).eval (-x) - (U ℝ m).eval (-x) := by rw [eval_rec]
      _ = 2 * (-x) * ((-1 : ℝ)^(m + 1) * (U ℝ (m + 1)).eval x)
            - ((-1 : ℝ)^m * (U ℝ m).eval x) := by
            norm_num
            ring_nf
            rw [ih m (Nat.le.step Nat.le.refl)]
            ring_nf
            simp
            -- Handle the inner evaluation term
            calc
              -(x * eval (-x) (U ℝ (1 + m)) * 2)
                  = -(x * ((-1)^(m+1) * eval x (U ℝ (m+1))) * 2) := by
                    congr 1
                    congr 1
                    congr 1
                    ring_nf
                    -- Prove the equality for the inner evaluation
                    have h_eq : eval (-x) (U ℝ (1 + ↑m)) = -(eval x (U ℝ (1 + ↑m)) * (-1)^m) := by
                      calc
                        eval (-x) (U ℝ (1 + ↑m))
                            = eval (-x) (U ℝ ↑(m+1)) := by simp; ring_nf
                        _   = (-1)^(m+1) * eval x (U ℝ ↑(m+1)) := by rw [ih (m + 1) Nat.le.refl]
                        _   = -(eval x (U ℝ (1+↑m)) * (-1)^m) := by ring_nf; simp
                    exact h_eq
              _ = x * eval x (U ℝ (m+1)) * (-1)^m * 2 := by
                    rw [pow_succ]
                    ring
            rw [Int.add_comm]
      _ = (-1 : ℝ)^m * (2 * x * (U ℝ (m + 1)).eval x - (U ℝ m).eval x) := by ring
      _ = (-1 : ℝ)^m * (U ℝ (m + 2)).eval x := by simp
      _ = (-1 : ℝ)^(m + 2) * (U ℝ (m + 2)).eval x := by
        norm_num
        ring_nf
        exact (true_or (x * eval x (U ℝ (1 + ↑m)) * 2 - eval x (U ℝ ↑m) = 0)).mpr trivial

/-- Cosine representation theorem for real numbers in a symmetric interval.
For any a > 0 and real number x with |x| ≤ a, there exists θ such that x = a cos θ. -/
lemma exists_cos_representation_general (x a : ℝ) (ha : a > 0) (h : |x| ≤ a) :
    ∃ (θ : ℝ), x = a * Real.cos θ := by
  -- Show that x/a is in the domain of arccos [-1, 1]
  have bound : -1 ≤ x/a ∧ x/a ≤ 1 := by
    constructor
    · -- Prove -1 ≤ x/a
      have h_left : -a ≤ x := abs_le.mp h |>.left
      calc
        -1 = (-a)/a := by field_simp [ne_of_gt ha]
        _ ≤ x/a := by exact (div_le_div_iff_of_pos_right ha).mpr h_left
    · -- Prove x/a ≤ 1
      have h_right : x ≤ a := abs_le.mp h |>.right
      calc
        x/a ≤ a/a := by exact (div_le_div_iff_of_pos_right ha).mpr h_right
        _ = 1 := by field_simp [ne_of_gt ha]
  -- Use θ = arccos(x/a)
  exact ⟨Real.arccos (x/a), by rw [Real.cos_arccos bound.1 bound.2]; field_simp [ne_of_gt ha]⟩

lemma exists_cos_representation (x : ℝ) (h : x^2 < 4) : ∃ (θ : ℝ), x = 2 * Real.cos θ := by
  -- Show that |x| < 2
  have : |x| < 2 := by
    rw [abs_lt]
    constructor <;> nlinarith
  -- Apply the general cosine representation theorem
  exact exists_cos_representation_general x 2 (by norm_num) (by linarith)

[Nick_adfor]

Now this two part haven't been finished. I think it needs a lot of calculation, which I'm not so familiar with. By the way, these two code should be pasted together so that they can work.

lemma Chebyshev.U_pos_of_ge_one {n : ℕ} {y : ℝ} (hy : 1 ≤ y) :
    0 < (U ℝ n).eval y := by
  induction' n using Nat.strong_induction_on with n IH
  rcases n with (rfl | rfl | m)
  · -- Case n = 0
    simp
  · -- Case n = 1
    simp
    linarith [hy]
  · -- Case n = m + 2
    -- Get recursion relation for Chebyshev polynomials
    have rec := Chebyshev_U_recursion (m+2) (by omega)
    simp at rec
    -- Inductive hypotheses for m and m+1
    have IH_m   := IH m (by omega)
    -- Inductive hypotheses for m and m+1
    have IH_m1  := IH (m + 1) (by omega)
    calc
      (U ℝ (m+2)).eval y
          = (2 * y) * (U ℝ (m+1)).eval y - (U ℝ m).eval y := by
          rw [rec]
          simp
          ring_nf
          exact (true_or (y = 0)).mpr trivial
      _ > (U ℝ (m+1)).eval y - (U ℝ m).eval y := by
          -- Calculation
          have h2y : 2 * y ≥ 2 := by linarith
          simp
          -- Calculation
          have pos : 0 < eval y (U ℝ (↑m + 1)) := by
            exact IH_m1
          nlinarith [pos, h2y]
      _ > 0 := by
          -- Calculation
          have : (U ℝ (m+1)).eval y > (U ℝ m).eval y := by
            rcases m with _ | m
            · -- Case m = 0
              simp
              linarith [hy]
            · -- Case m = m+1
              -- Calculation
                have h : m + 2 ≥ 2 := by
                  have : m ≥ 0 := Nat.zero_le m
                  linarith
                have rec' := Chebyshev_U_recursion (m + 2) h
                rw [Int.natCast_add] at rec'
                have rec'' : U ℝ (m + 2) = 2 * X * U ℝ (m + 1) - U ℝ m := by exact U_add_two ℝ ↑m
                have : U ℝ (↑(m + 1) + 1) = U ℝ (↑m + ↑2) := by exact rfl
                rw [this]
                rw [rec'']
                simp
                -- As before
                have h2y : 2 * y ≥ 2 := by linarith
                -- As before
                have pos : 0 < eval y (U ℝ (↑m + 1)) := by
                  exact IH_m
                -- Calculation
                have h : 2 * y - 1 ≥ 1 := by
                  nlinarith [hy]
                -- Calculation
                have IH_m := IH m (by omega)
                have IH_m_sub1 : 0 < (U ℝ (m - 1)).eval y := by
                  have : m - 1 < m + 2 := by omega
                  sorry
                have : 2 * y > 1 := by nlinarith [hy]
                sorry
          nlinarith

theorem Chebyshev_U_simple_roots (n : ℕ) (x : ℝ) :
    ((Chebyshev.U ℝ n).comp (Polynomial.C (1/2) * X)).eval x = 0 →
    ((Chebyshev.U ℝ n).comp (Polynomial.C (1/2) * X)).derivative.eval x ≠ 0 := by
  intro hx
  -- Evaluation lemmas at x = 1
  have U_at_1 := Chebyshev_U_eval_one n
  -- Evaluation lemmas at x = -1
  have U_at_neg1 := Chebyshev.U_eval_neg_one n
  -- Non-vanishing at endpoints for positive n
  have endpoints_nonzero : n > 0 →
      ((Chebyshev.U ℝ n).comp (Polynomial.C (1/2) * X)).eval 2 ≠ 0 ∧
      ((Chebyshev.U ℝ n).comp (Polynomial.C (1/2) * X)).eval (-2) ≠ 0 := by
    intro hn
    exact endpoints_nonzero n
  -- Show that root must be in open interval (-2, 2)
  have root_in_open_interval : x^2 < 4 := by
    by_cases h : n = 0
    · -- Case n = 0: polynomial is constant 1, no roots
      subst h
      simp at hx ⊢
    · -- Case n > 0
      -- Calculation
      have hn : n > 0 := Nat.pos_of_ne_zero h
      -- Calculation
      have ⟨h1, h2⟩ := endpoints_nonzero hn
      by_contra! H
      -- Calculation
      have : x ≤ -2 ∨ x ≥ 2 := by
        by_cases h : x ≥ 0
        · right
          nlinarith [H]
        · left
          -- Calculation
          have : x ≤ 0 := by linarith
          nlinarith [H]
      rcases this with (hx' | hx')
      · -- Case x ≤ -2
        -- Calculation
        have : (1/2 : ℝ) * x ≤ -1 := by linarith
        -- Calculation
        have mono : ∀ (y : ℝ), y ≤ -1 → (Chebyshev.U ℝ n).eval y ≠ 0 := by
          intro y hy
          rw [show y = -(-y) by ring] at hy ⊢
          -- Calculation
          have : -y ≥ 1 := by linarith
          -- Calculation
          have U_neg : (Chebyshev.U ℝ n).eval (-y) = (-1)^n * (Chebyshev.U ℝ n).eval y := by exact Chebyshev.U_neg n y
          -- Calculation
          have : (Chebyshev.U ℝ n).eval 1 = n + 1 := U_at_1
          -- Calculation
          have pos : 0 < (Chebyshev.U ℝ n).eval 1 := by linarith
          intro H
          simp at H
          simp [H] at U_neg
          sorry
        simp [Polynomial.eval_comp, eval_mul, eval_C, eval_X] at hx
        sorry
      · -- Case x ≥ 2
        have : ∀ (y : ℝ), y ≥ 1 → (Chebyshev.U ℝ n).eval y ≠ 0 := by
          intro y hy
          -- Calculation
          have pos : 0 < (Chebyshev.U ℝ n).eval y := by exact Chebyshev.U_pos_of_ge_one hy
          exact ne_of_gt pos
        simp at hx
        exact this ((2⁻¹ : ℝ) * x) (by linarith) hx
  -- Express x in cosine form: x = 2 cos θ
  have : ∃ (θ : ℝ), x = 2 * Real.cos θ := exists_cos_representation x root_in_open_interval

  rcases this with ⟨θ, hx_eq⟩
  rw [hx_eq]

  -- Derivative formula for scaled Chebyshev polynomial
  have deriv_formula : ((Chebyshev.U ℝ n).comp (Polynomial.C (1/2) * X)).derivative =
                     Polynomial.C ((n + 1 : ℝ) / 2) * (Chebyshev.U ℝ (n - 1)).comp (Polynomial.C (1/2) * X) := by
    sorry

  rw [deriv_formula, Polynomial.eval_mul, Polynomial.eval_C]

  refine mul_ne_zero (by norm_num; exact Nat.cast_add_one_ne_zero n) ?_
  simp [Polynomial.eval_comp, eval_mul, eval_C, eval_X]

  -- Show that U_{n-1}(cos θ) ≠ 0
  have : (Chebyshev.U ℝ (n - 1)).eval (cos θ) ≠ 0 := by
    intro H
    -- Evaluation lemmas for U_{n-1} at 1 and -1
    have U_prev_at_1 : (Chebyshev.U ℝ (n - 1)).eval (1 : ℝ) = n := by
      by_cases h : n = 0
      · subst h; simp
        -- Calculation
      · have := Chebyshev_U_eval_one (n - 1)
        simp at this ⊢
    -- Calculation
    have U_prev_at_neg1 : (Chebyshev.U ℝ (n - 1)).eval (-1 : ℝ) = (-1)^(n - 1) * n := by
      by_cases h : n = 0
      · subst h; simp
        -- Calculation
      · have := Chebyshev.U_eval_neg_one (n - 1)
        simp at this ⊢
        -- Calculation
        have n_ge_one : n ≥ 1 := Nat.one_le_iff_ne_zero.mpr h
        simpa [n_ge_one]
    -- We know U_n(cos θ) = 0 from hypothesis
    have h₀ : (Chebyshev.U ℝ n).eval (cos θ) = 0 := by
      simp [hx_eq] at hx
      exact hx
    -- And U_{n-1}(cos θ) = 0 from assumption
    have h₁ : (Chebyshev.U ℝ (n - 1)).eval (cos θ) = 0 := H
    -- But this leads to a contradiction
    have : (Chebyshev.U ℝ (n - 1)).eval (cos θ) ≠ 0 := by
      intro h₁
      -- Use trigonometric formulas for Chebyshev polynomials
      have U_formula := Polynomial.Chebyshev.U_real_cos θ (n : ℤ)
      -- Calculation
      have U_prev_formula := Polynomial.Chebyshev.U_real_cos θ ((n - 1 : ℕ) : ℤ)

      -- From h₀ we get sin((n+1)θ) = 0
      have hsin : Real.sin ((n+1) * θ) = 0 := by
        -- Calculation
        have U_formula' :
            (U ℝ n).eval (Real.cos θ) * Real.sin θ = Real.sin ((n+1) * θ) := by simp
        rw [h₀] at U_formula'
        simp only [zero_mul] at U_formula'
        exact id (Eq.symm U_formula')
      -- From h₁ we get sin(nθ) = 0
      have hsin' : Real.sin (n * θ) = 0 := by
        -- Calculation
        have U_prev_formula' :
            (U ℝ (n-1)).eval (Real.cos θ) * Real.sin θ = Real.sin (n * θ) := by simp
        rw [h₁] at U_prev_formula'
        simp only [zero_mul] at U_prev_formula'
        exact id (Eq.symm U_prev_formula')
      -- Express zeros using periodicity of sine
      obtain ⟨k, hk⟩ := Real.sin_eq_zero_iff.1 hsin
      obtain ⟨l, hl⟩ := Real.sin_eq_zero_iff.1 hsin'
      -- Calculation
      have : n*k = l*(n+1) := by
        -- From sin((n+1)θ) = 0 we get θ = kπ/(n+1)
        have hsin_eq1 : (n + 1) * θ = k * π := by
          rcases Real.sin_eq_zero_iff.mp hsin with ⟨k', hk'⟩
          -- Calculation
          have : θ = k' * π / (n + 1) := by
            field_simp
            nlinarith
          rw [this] at hk'
          field_simp at hk' ⊢
          nlinarith

        -- From sin(nθ) = 0 we get θ = lπ/n
        have hsin_eq2 : n * θ = l * π := by
          rcases Real.sin_eq_zero_iff.mp hsin' with ⟨l', hl'⟩
          -- Calculation
          have : θ = l' * π / n := by
            sorry
          rw [this] at hl'
          field_simp at hl' ⊢
          nlinarith

        -- Now we have nθ = lπ and (n+1)θ = kπ
        -- Subtract to get θ = (k - l)π
        have θ_eq : θ = (k - l) * π := by
          -- Calculation
          have : (n + 1) * θ - n * θ = k * π - l * π := by
            rw [hsin_eq1, hsin_eq2]
          ring_nf at this
          nlinarith

        -- Substitute back into nθ = lπ
        rw [θ_eq] at hsin_eq2
        -- Get n(k - l)π = lπ
        have : n * (k - l) * π = l * π := by
          nlinarith

        -- Cancel π (since π ≠ 0)
        have π_ne_zero : π ≠ 0 := by exact Real.pi_pos.ne.symm
        field_simp [π_ne_zero] at this ⊢
        -- Get n(k - l) = l
        sorry
      sorry
    exact this H

  exact this

[Nick_adfor]

Mathlib.Analysis.SpecialFunctions.Trigonometric.Chebyshev