Integrating the Chebyshev interpolant gives Clenshaw–Curtis quadrature. Its weights and error estimates follow from interpolation, moments and positivity.
2.8.1 Integrating the interpolant
Clenshaw–Curtis quadrature integrates the polynomial reconstructed from the cosine-grid samples. This gives both a coefficient formula and a nodal-weight formula.
For \(N\ge 1\), define \(\mathcal Q_Nf=\sum _{k=0}^Nc_k(f)m_k\), where \(m_k=\int _{-1}^1T_k\). Separately, set \(\mathcal Q_0f=2f(1)\) for the singleton interpolant.
/-- For `N>0`, Clenshaw–Curtis quadrature via the finite cosine transform (ATAP 19.2). -/
noncomputable def clenshawCurtis (N : ℕ) (f : ℝ → ℝ) : ℝ :=
if h : N = 0 then 2 * f 1 else
∑ k : Fin (N + 1),
extremaCoefficient ⟨N, Nat.pos_of_ne_zero h⟩ f k * integralMoment k.val
For \(N\ge 1\), \(\mathcal Q_Nf=\int _{-1}^1P_Nf\), and this quadrature is exact for polynomials of degree at most \(N\). If \(f\) is interval integrable and \(|f-P_Nf|\le E\) on \([-1,1]\), its quadrature error is at most \(2E\). If \(\rho {\gt}1\), \(F\) is analytic on the open Bernstein ellipse \(E_\rho \), and \(|F|\le M\) there, then
theorem Polynomial.Chebyshev.clenshawCurtis_eval {N : ℕ} (hN : N ≠ 0) {p : Polynomial ℝ} (hp : p.degree ≤ ↑N) :
(Polynomial.Chebyshev.clenshawCurtis N fun x => Polynomial.eval x p) = ∫ (x : ℝ) in -1..1, Polynomial.eval x p
Lean proof
/-- The `(N+1)`-point rule is exact through degree `N` (ATAP 19.1). -/
theorem clenshawCurtis_eval {N : ℕ} (hN : N ≠ 0) {p : ℝ[X]} (hp : p.degree ≤ N) :
clenshawCurtis N (fun x ↦ p.eval x) = ∫ x in (-1 : ℝ)..1, p.eval x := by
rw [clenshawCurtis_eq_integral_dctInterpolant, dctInterpolant_eval hN hp]
theorem Polynomial.Chebyshev.abs_integral_sub_clenshawCurtis_le {N : ℕ} {f : ℝ → ℝ} {E : ℝ}
(hf : IntervalIntegrable f MeasureTheory.volume (-1) 1)
(he : ∀ x ∈ Set.Icc (-1) 1, |f x - Polynomial.eval x (Polynomial.Chebyshev.dctInterpolant N f)| ≤ E) :
|(∫ (x : ℝ) in -1..1, f x) - Polynomial.Chebyshev.clenshawCurtis N f| ≤ 2 * E
Lean proof
/-- A uniform interpolation error controls the quadrature error by the interval length. -/
theorem abs_integral_sub_clenshawCurtis_le {N : ℕ} {f : ℝ → ℝ} {E : ℝ}
(hf : IntervalIntegrable f volume (-1) 1)
(he : ∀ x ∈ Icc (-1 : ℝ) 1, |f x - (dctInterpolant N f).eval x| ≤ E) :
|(∫ x in (-1 : ℝ)..1, f x) - clenshawCurtis N f| ≤ 2 * E := by
rw [clenshawCurtis_eq_integral_dctInterpolant,
← intervalIntegral.integral_sub hf ((dctInterpolant N f).continuous.intervalIntegrable _ _)]
have h := intervalIntegral.norm_integral_le_of_norm_le_const
(a := (-1 : ℝ)) (b := 1) (f := fun x ↦ f x - (dctInterpolant N f).eval x)
(fun x hx ↦ by
rw [Real.norm_eq_abs]
apply he x
simpa only [uIcc_of_le (by norm_num : (-1 : ℝ) ≤ 1)] using uIoc_subset_uIcc hx)
simpa only [Real.norm_eq_abs, sub_neg_eq_add, one_add_one_eq_two,
abs_of_pos (by norm_num : (0 : ℝ) < 2), mul_comm] using h
theorem Polynomial.Chebyshev.abs_integral_sub_clenshawCurtis_le_of_analytic {F : ℂ → ℂ} {ρ M : ℝ} (hρ : 1 < ρ)
(hF : AnalyticOnNhd ℂ F (Complex.bernsteinEllipse ⟨ρ, hρ⟩)) (hM : ∀ z ∈ Complex.bernsteinEllipse ⟨ρ, hρ⟩, ‖F z‖ ≤ M)
{N : ℕ} (hN : N ≠ 0) :
|(∫ (x : ℝ) in -1..1, (F ↑x).re) - Polynomial.Chebyshev.clenshawCurtis N fun x => (F ↑x).re| ≤
8 * M / (ρ ^ N * (ρ - 1))
Lean proof
/-- Baseline geometric Clenshaw–Curtis error, from ATAP 8.2. Analyticity is only
required in the open ellipse. No positivity theorem for the weights is assumed. -/
theorem abs_integral_sub_clenshawCurtis_le_of_analytic {F : ℂ → ℂ} {ρ M : ℝ}
(hρ : 1 < ρ) (hF : AnalyticOnNhd ℂ F (bernsteinEllipse ⟨ρ, hρ⟩))
(hM : ∀ z ∈ bernsteinEllipse ⟨ρ, hρ⟩, ‖F z‖ ≤ M) {N : ℕ} (hN : N ≠ 0) :
|(∫ x in (-1 : ℝ)..1, (F x).re) - clenshawCurtis N (fun x ↦ (F x).re)| ≤
8 * M / (ρ ^ N * (ρ - 1)) := by
have hc : ContinuousOn (fun x : ℝ ↦ (F x).re) (Icc (-1) 1) :=
Complex.continuous_re.continuousOn.comp
(hF.continuousOn.comp Complex.continuous_ofReal.continuousOn
(fun x hx ↦ ofReal_mem_bernsteinEllipse hρ hx)) (fun _ _ ↦ mem_univ _)
have hi : IntervalIntegrable (fun x : ℝ ↦ (F x).re) volume (-1) 1 :=
ContinuousOn.intervalIntegrable (by simpa only [uIcc_of_le (by norm_num : (-1 : ℝ) ≤ 1)] using hc)
have h := abs_integral_sub_clenshawCurtis_le hi (E := 4 * M / (ρ ^ N * (ρ - 1)))
(fun x hx ↦ by
rw [eval_dctInterpolant_re hN]
exact (abs_re_le_norm (F x - (extremaInterpolant N (fun t ↦ F t)).eval (x : ℂ))).trans
(norm_sub_extremaInterpolant_le_of_analytic hρ hF hM hN hx))
convert h using 1
ring
The order-zero rule integrates its constant interpolant: \(\mathcal Q_0f=2f(1)\).
Lean Formalisation
theorem Polynomial.Chebyshev.clenshawCurtis_zero (f : ℝ → ℝ) : Polynomial.Chebyshev.clenshawCurtis 0 f = 2 * f 1
Lean proof
/-- The singleton interpolant integrates to twice its sole sample. -/
theorem clenshawCurtis_zero (f : ℝ → ℝ) : clenshawCurtis 0 f = 2 * f 1 := by
simp [clenshawCurtis]
For \(N\ge 1\) and \(0\le j\le N\), write \(x_j=\cos (j\pi /N)\), \(d_j=1/2\) at \(j=0,N\) and \(d_j=1\) otherwise, and let \(m_k\) be the classical moment from Definition 2.7.10. Define
For \(N\ge 1\) and arbitrary real nodal data, the coefficient-defined rule is \(Q_Nf=\sum _{j=0}^Nw_jf(x_j)\). Each weight is the integral of its Lagrange cardinal polynomial. In the explicit formula, replace \(T_k(x_j)\) by \(\cos (kj\pi /N)\) and \(m_k\) by \(2/(1-k^2)\) for even \(k\) and zero for odd \(k\). Thus the highest term contributes when \(N\) is even and vanishes when \(N\) is odd.
theorem Polynomial.Chebyshev.extremaSum_eq_sum_endpointWeight (N : { n // 0 < n }) (f : ℕ → ℝ) :
Polynomial.Chebyshev.extremaSum N f = ∑ j, Polynomial.Chebyshev.extremaEndpointWeight N j * f ↑j
Lean proof
/-- The half-weighted sum as an ordinary finite weighted sum. -/
theorem extremaSum_eq_sum_endpointWeight (N : {n : ℕ // 0 < n}) (f : ℕ → ℝ) :
extremaSum N f = ∑ j : Fin (N.val + 1), extremaEndpointWeight N j * f j := by
rcases N with ⟨N, hNp⟩
have hN : N ≠ 0 := Nat.ne_of_gt hNp
simp only [extremaEndpointWeight]
rw [Fin.sum_univ_eq_sum_range (fun j : ℕ => (if j = 0 ∨ j = N then (1 / 2 : ℝ) else 1) * f j)]
rw [Finset.sum_range_succ]
have he : ∑ j ∈ Finset.range N, (if j = 0 ∨ j = N then (1 / 2 : ℝ) else 1) * f j =
(∑ j ∈ Finset.range N, f j) - f 0 / 2 := by
have hj (j : ℕ) (hj : j ∈ Finset.range N) :
(if j = 0 ∨ j = N then (1 / 2 : ℝ) else 1) * f j = f j - (if j = 0 then f 0 / 2 else 0) := by
have hjN : j ≠ N := ne_of_lt (Finset.mem_range.mp hj)
by_cases hz : j = 0
· subst j; simp; ring
· simp [hz, hjN]
rw [Finset.sum_congr rfl hj, Finset.sum_sub_distrib]
simp [Nat.pos_of_ne_zero hN]
rw [he]
simp only [or_true, ite_true]
dsimp only [extremaSum]
ring
theorem Polynomial.Chebyshev.sum_clenshawCurtisWeight (N : { n // 0 < n }) :
∑ j, Polynomial.Chebyshev.clenshawCurtisWeight N j = 2
Lean proof
/-- The weights have the mass of the original interval. -/
theorem sum_clenshawCurtisWeight (N : {n : ℕ // 0 < n}) :
∑ j : Fin (N.val + 1), clenshawCurtisWeight N j = 2 := by
have he := clenshawCurtis_eval (Nat.ne_of_gt N.property) (p := (1 : ℝ[X])) (by simp)
rw [clenshawCurtis_eq_sum_weights N] at he
simpa only [eval_one, mul_one, intervalIntegral.integral_const, sub_neg_eq_add,
one_add_one_eq_two, smul_eq_mul] using he
theorem Polynomial.Chebyshev.clenshawCurtisWeight_last (N : { n // 0 < n }) :
Polynomial.Chebyshev.clenshawCurtisWeight N (Fin.last ↑N) = Polynomial.Chebyshev.clenshawCurtisWeight N 0
Lean proof
/-- Both endpoint weights agree. -/
theorem clenshawCurtisWeight_last (N : {n : ℕ // 0 < n}) :
clenshawCurtisWeight N (Fin.last N.val) = clenshawCurtisWeight N 0 := by
have hN : N.val ≠ 0 := Nat.ne_of_gt N.property
have he (k : ℕ) : integralMoment k * (-1 : ℝ) ^ k = integralMoment k := by
rcases Nat.even_or_odd k with hk | hk
· rw [hk.neg_one_pow, mul_one]
· simp [integralMoment_eq_ite, (Nat.not_even_iff_odd.mpr hk)]
simp only [clenshawCurtisWeight_eq_sum, classicalNode, Fin.val_last, Fin.val_zero, extremaEndpointWeight, or_true, true_or, ite_true,
node_eq_one, node_eq_neg_one hN, T_eval_one, mul_one, T_eval_neg_one,
Int.cast_negOnePow_natCast, mul_assoc, he]
The explicit weight formula uses the unweighted Chebyshev moments. For a positive grid order, the following code computes the weights and displays their positivity, sum and endpoint value.
m(k)=iseven(k)?2/(1-k^2):0.0# integral of T_kfunctionccweights(N)@assertN>=1[sum((k==0||k==N?1:2)/N*(j==0||j==N?1/2:1.0)*cos(k*j*pi/N)*m(k)forkin0:N)forjin0:N]endw=ccweights(8)(minimum(w)>0,sum(w),w[1])# positivity, total weight and endpoint weight
2.8.3 Error estimates and interval scaling
A polynomial competitor bounds the quadrature error. The sharper analytic estimate uses the even coefficient tail.
If \(|f(x_j)-g(x_j)|\le \varepsilon \) at all nodes, then \(|Q_Nf-Q_Ng|\le 2\varepsilon \). If \(f\) is interval-integrable and a real polynomial \(p\) of degree at most \(N\) satisfies \(|f(x)-p(x)|\le \varepsilon \) throughout \([-1,1]\), then
theorem Polynomial.Chebyshev.abs_integral_sub_clenshawCurtis_le_of_competitor {N : ℕ} (hN : N ≠ 0) {f : ℝ → ℝ}
(hf : IntervalIntegrable f MeasureTheory.volume (-1) 1) {p : Polynomial ℝ} (hp : p.degree ≤ ↑N) {ε : ℝ}
(hε : ∀ x ∈ Set.Icc (-1) 1, |f x - Polynomial.eval x p| ≤ ε) :
|(∫ (x : ℝ) in -1..1, f x) - Polynomial.Chebyshev.clenshawCurtis N f| ≤ 4 * ε
Lean proof
/-- Every degree-N uniform competitor gives a direct quadrature bound of 4ε.
Only interval integrability of the data is required. -/
theorem abs_integral_sub_clenshawCurtis_le_of_competitor {N : ℕ} (hN : N ≠ 0)
{f : ℝ → ℝ} (hf : IntervalIntegrable f volume (-1) 1) {p : ℝ[X]}
(hp : p.degree ≤ N) {ε : ℝ}
(hε : ∀ x ∈ Icc (-1 : ℝ) 1, |f x - p.eval x| ≤ ε) :
|(∫ x in (-1 : ℝ)..1, f x) - clenshawCurtis N f| ≤ 4 * ε := by
have hq := abs_clenshawCurtis_sub_le hN (f := f) (g := fun x ↦ p.eval x) (ε := ε)
(fun j (_ : j ≤ N) ↦ hε (node N j) node_mem_Icc)
have hi : |∫ x in (-1 : ℝ)..1, f x - p.eval x| ≤ 2 * ε := by
have he := intervalIntegral.norm_integral_le_of_norm_le_const
(a := (-1 : ℝ)) (b := 1) (f := fun x ↦ f x - p.eval x)
(fun x hx ↦ by
rw [Real.norm_eq_abs]
apply hε x
simpa only [uIcc_of_le (by norm_num : (-1 : ℝ) ≤ 1)] using uIoc_subset_uIcc hx)
simpa only [Real.norm_eq_abs, sub_neg_eq_add, one_add_one_eq_two,
abs_of_pos (by norm_num : (0 : ℝ) < 2), mul_comm] using he
have he : (∫ x in (-1 : ℝ)..1, f x) - clenshawCurtis N f =
(∫ x in (-1 : ℝ)..1, f x - p.eval x) -
(clenshawCurtis N f - clenshawCurtis N (fun x ↦ p.eval x)) := by
rw [intervalIntegral.integral_sub hf (p.continuous.intervalIntegrable _ _),
clenshawCurtis_eval hN hp]
ring
rw [he]
exact (abs_sub _ _).trans ((add_le_add hi hq).trans_eq (by ring))
If \(N\) is even, replace \(\rho ^{N-1}\) by \(\rho ^N\). More generally, any \(1\le m\) with \(2m\le N\) gives the denominator \(15\rho ^{2m}(\rho ^2-1)\). The underlying series bound uses \(|a_k|\le Cq^k\), \(0\le q{\lt}1\), to give \(\frac{32}{15}Cq^{2(m+1)}/(1-q^2)\).
For \(a{\lt}b\) and \(N\ge 1\), \(P_N(f\circ A)\circ B\) interpolates \(f\) at the mapped extrema \(A(x_j)\). Its integral is the scaled rule \(\frac{b-a}{2}\mathcal Q_N(f\circ A)\), which is exact for every physical polynomial of degree at most \(N\). If \(f\circ A\) is interval integrable and its canonical interpolation error is at most \(E\) on \([-1,1]\), the physical quadrature error is at most \((b-a)E\). In particular, let \(\rho {\gt}1\), let \(F\) be analytic on the open ellipse \(E_\rho \), and let \(|F|\le M\) there. For \(f(x)=\Re F(B(x))\) the error is at most
\[ \frac{4(b-a)M}{\rho ^N(\rho -1)}. \]
Reference: [ATAP, Chapters 8 and 19], with affine scaling.
Lean Formalisation
theorem Polynomial.Chebyshev.eval_dctInterpolant_interval_node {a b : ℝ} (hab : a < b) {N : ℕ} (hN : N ≠ 0) (f : ℝ → ℝ) {j : ℕ}
(hj : j ≤ N) :
Polynomial.eval
(Polynomial.Chebyshev.intervalCoord a ⟨b, ⋯⟩
(Polynomial.Chebyshev.intervalMap a b (Polynomial.Chebyshev.node N j)))
(Polynomial.Chebyshev.dctInterpolant N fun t => f (Polynomial.Chebyshev.intervalMap a b t)) =
f (Polynomial.Chebyshev.intervalMap a b (Polynomial.Chebyshev.node N j))
Lean proof
/-- The DCT polynomial pulled back to `[a,b]` matches all the mapped extrema. -/
theorem eval_dctInterpolant_interval_node {a b : ℝ} (hab : a < b)
{N : ℕ} (hN : N ≠ 0) (f : ℝ → ℝ) {j : ℕ} (hj : j ≤ N) :
(dctInterpolant N (fun t ↦ f (intervalMap a b t))).eval
(intervalCoord a ⟨b, ne_of_lt hab⟩ (intervalMap a b (node N j))) = f (intervalMap a b (node N j)) := by
rw [intervalCoord_intervalMap hab]
exact eval_dctInterpolant_node hN _ hj
theorem Polynomial.Chebyshev.integral_dctInterpolant_interval {a b : ℝ} (hab : a < b) (N : ℕ) (f : ℝ → ℝ) :
∫ (x : ℝ) in a..b,
Polynomial.eval (Polynomial.Chebyshev.intervalCoord a ⟨b, ⋯⟩ x)
(Polynomial.Chebyshev.dctInterpolant N fun t => f (Polynomial.Chebyshev.intervalMap a b t)) =
(b - a) / 2 * Polynomial.Chebyshev.clenshawCurtis N fun t => f (Polynomial.Chebyshev.intervalMap a b t)
Lean proof
/-- Integrating the transported DCT interpolant gives the scaled canonical rule. -/
theorem integral_dctInterpolant_interval {a b : ℝ} (hab : a < b) (N : ℕ) (f : ℝ → ℝ) :
(∫ x in a..b, (dctInterpolant N (fun t ↦ f (intervalMap a b t))).eval
(intervalCoord a ⟨b, ne_of_lt hab⟩ x)) =
(b - a) / 2 * clenshawCurtis N (fun t ↦ f (intervalMap a b t)) := by
have hi : IntervalIntegrable (fun x ↦
(dctInterpolant N (fun t ↦ f (intervalMap a b t))).eval
(intervalCoord a ⟨b, ne_of_lt hab⟩ x)) volume a b := by
apply Continuous.intervalIntegrable
exact (dctInterpolant N _).continuous.comp (by unfold intervalCoord; fun_prop)
rw [integral_comp_intervalMap a b ⟨_, hi⟩]
simp only [intervalCoord_intervalMap hab, clenshawCurtis_eq_integral_dctInterpolant]
theorem Polynomial.Chebyshev.clenshawCurtis_interval_eval {a b : ℝ} (hab : a < b) {N : ℕ} (hN : N ≠ 0) {p : Polynomial ℝ}
(hp : p.degree ≤ ↑N) :
((b - a) / 2 *
Polynomial.Chebyshev.clenshawCurtis N fun t => Polynomial.eval (Polynomial.Chebyshev.intervalMap a b t) p) =
∫ (x : ℝ) in a..b, Polynomial.eval x p
Lean proof
/-- Clenshaw--Curtis remains exact on physical polynomials of degree at most `N`. -/
theorem clenshawCurtis_interval_eval {a b : ℝ} (hab : a < b)
{N : ℕ} (hN : N ≠ 0) {p : ℝ[X]} (hp : p.degree ≤ N) :
(b - a) / 2 * clenshawCurtis N (fun t ↦ p.eval (intervalMap a b t)) =
∫ x in a..b, p.eval x := by
let q := p.comp (Polynomial.C ((b - a) / 2) * X + Polynomial.C ((a + b) / 2))
have hc : (b - a) / 2 ≠ 0 := ne_of_gt (half_pos (sub_pos.mpr hab))
have hd : (Polynomial.C ((b - a) / 2) * X + Polynomial.C ((a + b) / 2)).degree = 1 := by
rw [degree_add_eq_left_of_degree_lt (degree_C_lt_degree_C_mul_X hc), degree_C_mul_X hc]
have hq : q.degree ≤ N := by
dsimp only [q]
rwa [degree_comp (by rw [hd]; norm_num), hd, mul_one]
have he (t : ℝ) : q.eval t = p.eval (intervalMap a b t) := by simp [q, intervalMap]
rw [integral_comp_intervalMap a b ⟨_, p.continuous.intervalIntegrable _ _⟩]
congr 1
simpa only [he] using clenshawCurtis_eval hN hq
theorem Polynomial.Chebyshev.abs_integral_sub_clenshawCurtis_interval_le {a b E : ℝ} (hab : a < b) {N : ℕ} {f : ℝ → ℝ}
(hf : IntervalIntegrable (fun t => f (Polynomial.Chebyshev.intervalMap a b t)) MeasureTheory.volume (-1) 1)
(he :
∀ t ∈ Set.Icc (-1) 1,
|f (Polynomial.Chebyshev.intervalMap a b t) -
Polynomial.eval t
(Polynomial.Chebyshev.dctInterpolant N fun u => f (Polynomial.Chebyshev.intervalMap a b u))| ≤
E) :
|(∫ (x : ℝ) in a..b, f x) -
(b - a) / 2 * Polynomial.Chebyshev.clenshawCurtis N fun t => f (Polynomial.Chebyshev.intervalMap a b t)| ≤
(b - a) * E
Lean proof
/-- A canonical interpolation error `E` gives physical quadrature error at most
`(b-a)*E`, with the Jacobian accounted for exactly. -/
theorem abs_integral_sub_clenshawCurtis_interval_le {a b E : ℝ} (hab : a < b)
{N : ℕ} {f : ℝ → ℝ}
(hf : IntervalIntegrable (fun t ↦ f (intervalMap a b t)) volume (-1) 1)
(he : ∀ t ∈ Icc (-1 : ℝ) 1,
|f (intervalMap a b t) - (dctInterpolant N (fun u ↦ f (intervalMap a b u))).eval t| ≤ E) :
|(∫ x in a..b, f x) - (b - a) / 2 * clenshawCurtis N (fun t ↦ f (intervalMap a b t))| ≤
(b - a) * E := by
rw [integral_comp_intervalMap a b ⟨f, integrable_of_affine hab hf⟩, ← mul_sub, abs_mul,
abs_of_pos (half_pos (sub_pos.mpr hab))]
have h := mul_le_mul_of_nonneg_left (abs_integral_sub_clenshawCurtis_le hf he)
(half_pos (sub_pos.mpr hab)).le
convert h using 1
ring
theorem Polynomial.Chebyshev.abs_integral_sub_clenshawCurtis_interval_le_of_analytic {a b : ℝ} (hab : a < b) {F : ℂ → ℂ}
{ρ M : ℝ} (hρ : 1 < ρ) (hF : AnalyticOnNhd ℂ F (Complex.bernsteinEllipse ⟨ρ, hρ⟩))
(hM : ∀ z ∈ Complex.bernsteinEllipse ⟨ρ, hρ⟩, ‖F z‖ ≤ M) {N : ℕ} (hN : N ≠ 0) :
|(∫ (x : ℝ) in a..b, (F ↑(Polynomial.Chebyshev.intervalCoord a ⟨b, ⋯⟩ x)).re) -
(b - a) / 2 * Polynomial.Chebyshev.clenshawCurtis N fun t => (F ↑t).re| ≤
4 * (b - a) * M / (ρ ^ N * (ρ - 1))
Lean proof
/-- Geometric error on `[a,b]` when `F` is the analytic extension in canonical
coordinates. No boundary analyticity is assumed (ATAP, Theorem 8.2). -/
theorem abs_integral_sub_clenshawCurtis_interval_le_of_analytic
{a b : ℝ} (hab : a < b) {F : ℂ → ℂ} {ρ M : ℝ} (hρ : 1 < ρ)
(hF : AnalyticOnNhd ℂ F (bernsteinEllipse ⟨ρ, hρ⟩))
(hM : ∀ z ∈ bernsteinEllipse ⟨ρ, hρ⟩, ‖F z‖ ≤ M) {N : ℕ} (hN : N ≠ 0) :
|(∫ x in a..b, (F (intervalCoord a ⟨b, ne_of_lt hab⟩ x)).re) -
(b - a) / 2 * clenshawCurtis N (fun t ↦ (F t).re)| ≤
4 * (b - a) * M / (ρ ^ N * (ρ - 1)) := by
have hc : ContinuousOn (fun t : ℝ ↦ (F t).re) (Icc (-1) 1) :=
Complex.continuous_re.continuousOn.comp
(hF.continuousOn.comp Complex.continuous_ofReal.continuousOn
(fun t ht ↦ ofReal_mem_bernsteinEllipse hρ ht)) (fun _ _ ↦ mem_univ _)
have hi : IntervalIntegrable (fun t : ℝ ↦ (F t).re) volume (-1) 1 :=
ContinuousOn.intervalIntegrable (by simpa only [uIcc_of_le (by norm_num : (-1 : ℝ) ≤ 1)] using hc)
have hp : IntervalIntegrable (fun x ↦ (F (intervalCoord a ⟨b, ne_of_lt hab⟩ x)).re)
volume a b := integrable_of_affine hab (by simpa only [intervalCoord_intervalMap hab] using hi)
rw [integral_comp_intervalMap a b ⟨_, hp⟩]
simp only [intervalCoord_intervalMap hab]
rw [← mul_sub, abs_mul, abs_of_pos (half_pos (sub_pos.mpr hab))]
have h := mul_le_mul_of_nonneg_left
(abs_integral_sub_clenshawCurtis_le_of_analytic hρ hF hM hN)
(half_pos (sub_pos.mpr hab)).le
convert h using 1
ring
Suppose \(\sum _{k\ge 0}|a_k|{\lt}\infty \) and \(f(x)=\sum _{k\ge 0}a_kT_k(x)\) on \([-1,1]\). Integration and quadrature commute with this series. If \(N\ge 1\) and \(2m\le N\), then the error has the absolutely convergent expansion
theorem Polynomial.Chebyshev.hasSum_integralMoment {a : ℕ → ℝ} {f : ℝ → ℝ} (ha : Summable fun k => |a k|)
(hf : ∀ x ∈ Set.Icc (-1) 1, HasSum (fun k => a k * Polynomial.eval x (Polynomial.Chebyshev.T ℝ ↑k)) (f x)) :
HasSum (fun k => a k * Polynomial.Chebyshev.integralMoment k) (∫ (x : ℝ) in -1..1, f x)
Lean proof
/-- Absolute coefficient summability permits integration of a Chebyshev series. -/
theorem hasSum_integralMoment {a : ℕ → ℝ} {f : ℝ → ℝ}
(ha : Summable (fun k ↦ |a k|))
(hf : ∀ x ∈ Icc (-1 : ℝ) 1,
HasSum (fun k : ℕ ↦ a k * (T ℝ k).eval x) (f x)) :
HasSum (fun k ↦ a k * integralMoment k) (∫ x in (-1 : ℝ)..1, f x) := by
have h := intervalIntegral.hasSum_integral_of_dominated_convergence
(a := (-1 : ℝ)) (b := 1) (μ := volume)
(F := fun k x ↦ a k * (T ℝ k).eval x) (f := f) (fun k _ ↦ |a k|)
(fun k ↦ (continuous_const.mul (T ℝ k).continuous).aestronglyMeasurable)
(fun k ↦ Filter.Eventually.of_forall fun x hx ↦ by
have hx' : x ∈ Icc (-1 : ℝ) 1 := by
simpa only [uIcc_of_le (by norm_num : (-1 : ℝ) ≤ 1)] using uIoc_subset_uIcc hx
rw [Real.norm_eq_abs, abs_mul]
exact mul_le_of_le_one_right (abs_nonneg _) (abs_eval_T_real_le_one _ (abs_le.mpr hx')))
(Filter.Eventually.of_forall fun _ _ ↦ ha) intervalIntegrable_const
(Filter.Eventually.of_forall fun x hx ↦ hf x (by
simpa only [uIcc_of_le (by norm_num : (-1 : ℝ) ≤ 1)] using uIoc_subset_uIcc hx))
simpa only [intervalIntegral.integral_const_mul, ← integralMoment_eq_integral] using h
theorem Polynomial.Chebyshev.hasSum_clenshawCurtis {N : ℕ} (hN : N ≠ 0) {a : ℕ → ℝ} {f : ℝ → ℝ}
(hf :
∀ j ≤ N,
HasSum (fun k => a k * Polynomial.eval (Polynomial.Chebyshev.node N j) (Polynomial.Chebyshev.T ℝ ↑k))
(f (Polynomial.Chebyshev.node N j))) :
HasSum (fun k => a k * Polynomial.Chebyshev.clenshawCurtis N fun x => Polynomial.eval x (Polynomial.Chebyshev.T ℝ ↑k))
(Polynomial.Chebyshev.clenshawCurtis N f)
Lean proof
/-- Quadrature commutes with any series summable at its nodes. -/
theorem hasSum_clenshawCurtis {N : ℕ} (hN : N ≠ 0) {a : ℕ → ℝ} {f : ℝ → ℝ}
(hf : ∀ j ≤ N,
HasSum (fun k : ℕ ↦ a k * (T ℝ k).eval (node N j)) (f (node N j))) :
HasSum (fun k ↦ a k * clenshawCurtis N (fun x ↦ (T ℝ k).eval x))
(clenshawCurtis N f) := by
simp_rw [clenshawCurtis_eq_sum_weights ⟨N, Nat.pos_of_ne_zero hN⟩, Finset.mul_sum]
simpa only [mul_left_comm, classicalNode] using
hasSum_sum (s := (Finset.univ : Finset (Fin (N + 1))))
(fun j hj ↦ (hf j (Nat.le_of_lt_succ j.isLt)).mul_left (clenshawCurtisWeight ⟨N, Nat.pos_of_ne_zero hN⟩ j))
theorem Polynomial.Chebyshev.hasSum_clenshawCurtis_error_even_tail {N m : ℕ} (hN : N ≠ 0) (hm : 2 * m ≤ N) {a : ℕ → ℝ}
{f : ℝ → ℝ} (ha : Summable fun k => |a k|)
(hf : ∀ x ∈ Set.Icc (-1) 1, HasSum (fun k => a k * Polynomial.eval x (Polynomial.Chebyshev.T ℝ ↑k)) (f x)) :
HasSum
(fun k =>
a (2 * (k + (m + 1))) *
(Polynomial.Chebyshev.integralMoment (2 * (k + (m + 1))) -
Polynomial.Chebyshev.clenshawCurtis N fun x =>
Polynomial.eval x (Polynomial.Chebyshev.T ℝ ↑(2 * (k + (m + 1))))))
((∫ (x : ℝ) in -1..1, f x) - Polynomial.Chebyshev.clenshawCurtis N f)
Lean proof
/-- Only even coefficients above the exactness degree contribute to the error.
The integer `m` may be any cutoff with `2m ≤ N`. -/
theorem hasSum_clenshawCurtis_error_even_tail {N m : ℕ} (hN : N ≠ 0) (hm : 2 * m ≤ N)
{a : ℕ → ℝ} {f : ℝ → ℝ} (ha : Summable (fun k ↦ |a k|))
(hf : ∀ x ∈ Icc (-1 : ℝ) 1,
HasSum (fun k : ℕ ↦ a k * (T ℝ k).eval x) (f x)) :
HasSum (fun k ↦ a (2 * (k + (m + 1))) *
(integralMoment (2 * (k + (m + 1))) -
clenshawCurtis N (fun x ↦ (T ℝ (↑(2 * (k + (m + 1)) : ℕ))).eval x)))
((∫ x in (-1 : ℝ)..1, f x) - clenshawCurtis N f) := by
let d (k : ℕ) := a k *
(integralMoment k - clenshawCurtis N (fun x ↦ (T ℝ k).eval x))
have hs : HasSum d ((∫ x in (-1 : ℝ)..1, f x) - clenshawCurtis N f) := by
simpa only [d, mul_sub] using
(hasSum_integralMoment ha hf).sub (hasSum_clenshawCurtis hN (fun j _ ↦ hf _ node_mem_Icc))
have he : Summable (fun k ↦ d (2 * k)) :=
hs.summable.comp_injective (fun _ _ h ↦ by omega)
have ho : ∀ k, d (2 * k + 1) = 0 := by
intro k
simp only [d, integralMoment_odd, clenshawCurtis_eval_T_odd hN, sub_self, mul_zero]
have heq : (∑' k, d (2 * k)) = (∫ x in (-1 : ℝ)..1, f x) - clenshawCurtis N f := by
have ht := tsum_even_add_odd he (hs.summable.comp_injective (fun _ _ h ↦ by omega))
simpa only [ho, tsum_zero, add_zero, hs.tsum_eq] using ht
have hz : ∑ k ∈ Finset.range (m + 1), d (2 * k) = 0 := by
apply Finset.sum_eq_zero
intro k hk
have hkN : 2 * k ≤ N := by simp only [Finset.mem_range] at hk; omega
have hp : (T ℝ (↑(2 * k : ℕ))).degree ≤ N := by
simpa only [degree_T, Int.natAbs_natCast, Nat.cast_le] using hkN
simp only [d, clenshawCurtis_eval hN hp, ← integralMoment_eq_integral, sub_self, mul_zero]
have ht := he.sum_add_tsum_nat_add (m + 1)
rw [hz, zero_add, heq] at ht
exact ht ▸ ((summable_nat_add_iff (m + 1)).2 he).hasSum
In particular, every finite absolute moment sum is strictly below \(3\). For \(N\ge 1\), the endpoint-mode sum satisfies \(\sum _{k=0}^Nb_km_k=2\sum _{k=0}^Nm_k-2-m_N\).
theorem Polynomial.Chebyshev.sum_integralMoment (N : ℕ) :
∑ k ∈ Finset.range (N + 1), Polynomial.Chebyshev.integralMoment k =
4 - ∑ k ∈ Finset.range (N + 1), |Polynomial.Chebyshev.integralMoment k|
Lean proof
/-- The signed moment sum is the complementary absolute moment sum. -/
theorem sum_integralMoment (N : ℕ) :
∑ k ∈ Finset.range (N + 1), integralMoment k =
4 - ∑ k ∈ Finset.range (N + 1), |integralMoment k| := by
have ht (k : ℕ) : integralMoment k =
(if k = 0 then (4 : ℝ) else 0) - |integralMoment k| := by
by_cases hk : k = 0
· subst k; norm_num [integralMoment_eq_ite]
· rw [ite_eq_right hk, zero_sub, abs_of_nonpos (integralMoment_nonpos (Nat.pos_of_ne_zero hk))]
ring
conv_lhs => rw [Finset.sum_congr rfl (fun k _ ↦ ht k)]
simp only [Finset.sum_sub_distrib, Finset.sum_ite_eq', Finset.mem_range,
Nat.zero_lt_succ, ite_true]
theorem Polynomial.Chebyshev.sum_endpoint_modes {N : ℕ} (hN : N ≠ 0) :
∑ k ∈ Finset.range (N + 1), (if k = 0 ∨ k = N then 1 else 2) * Polynomial.Chebyshev.integralMoment k =
2 * ∑ k ∈ Finset.range (N + 1), Polynomial.Chebyshev.integralMoment k - 2 - Polynomial.Chebyshev.integralMoment N
Lean proof
/-- Removing half of the constant and highest frequencies gives the endpoint kernel. -/
theorem sum_endpoint_modes {N : ℕ} (hN : N ≠ 0) :
(∑ k ∈ Finset.range (N + 1),
(if k = 0 ∨ k = N then (1 : ℝ) else 2) * integralMoment k) =
2 * (∑ k ∈ Finset.range (N + 1), integralMoment k) - 2 - integralMoment N := by
have he (k : ℕ) : (if k = 0 ∨ k = N then (1 : ℝ) else 2) * integralMoment k =
2 * integralMoment k - (if k = 0 then integralMoment 0 else 0) -
(if k = N then integralMoment N else 0) := by
by_cases h0 : k = 0
· subst k; norm_num [Ne.symm hN, integralMoment_eq_ite]
· by_cases hk : k = N
· subst k; simp [hN]; ring
· simp [h0, hk]
rw [Finset.sum_congr rfl (fun k _ ↦ he k)]
norm_num [Finset.sum_sub_distrib, ← Finset.mul_sum, integralMoment_eq_ite]
simp only [Finset.mul_sum, mul_ite, mul_zero]
For \(N\ge 1\), the folded index \(A_N(k)\) has the same parity as \(k\), and \(Q_NT_k=m_{A_N(k)}\). In particular, \(Q_NT_{2k+1}=0\) for every \(k\ge 0\), even beyond the polynomial exactness degree.
theorem Polynomial.Chebyshev.aliasIndex_mod_two {N : ℕ} (hN : N ≠ 0) (k : ℕ) : Polynomial.Chebyshev.aliasIndex N k % 2 = k % 2
Lean proof
/-- Folding a frequency at the extrema preserves its parity (ATAP, Chapter 4). -/
theorem aliasIndex_mod_two {N : ℕ} (hN : N ≠ 0) (k : ℕ) :
aliasIndex N k % 2 = k % 2 := by
have hm := Nat.mod_lt k (show 0 < 2 * N by omega)
have hp : (k % (2 * N)) % 2 = k % 2 :=
Nat.mod_mod_of_dvd k (dvd_mul_right 2 N)
simp only [aliasIndex, ite_eq_right hN]
rcases le_total (k % (2 * N)) (2 * N - k % (2 * N)) with h | h
· rw [min_eq_left h]
exact hp
· rw [min_eq_right h]
omega
theorem Polynomial.Chebyshev.clenshawCurtis_eval_T_alias {N : ℕ} (hN : N ≠ 0) (k : ℕ) :
(Polynomial.Chebyshev.clenshawCurtis N fun x => Polynomial.eval x (Polynomial.Chebyshev.T ℝ ↑k)) =
Polynomial.Chebyshev.integralMoment (Polynomial.Chebyshev.aliasIndex N k)
Lean proof
/-- The quadrature of a Chebyshev mode is the moment of its folded frequency. -/
theorem clenshawCurtis_eval_T_alias {N : ℕ} (hN : N ≠ 0) (k : ℕ) :
clenshawCurtis N (fun x ↦ (T ℝ k).eval x) = integralMoment (aliasIndex N k) := by
simp only [clenshawCurtis, dite_eq_right hN]
have he (j : Fin (N + 1)) :=
extremaCoefficient_eval_T_alias hN (Nat.le_of_lt_succ j.isLt) k
simp_rw [he]
rw [Fin.sum_univ_eq_sum_range
(fun j ↦ (if aliasIndex N k = j then (1 : ℝ) else 0) * integralMoment j)]
simp [ite_mul, Finset.sum_ite_eq, Finset.mem_range,
Nat.lt_succ_iff.mpr (aliasIndex_le hN k)]
theorem Polynomial.Chebyshev.clenshawCurtis_eval_T_odd {N : ℕ} (hN : N ≠ 0) (k : ℕ) :
(Polynomial.Chebyshev.clenshawCurtis N fun x => Polynomial.eval x (Polynomial.Chebyshev.T ℝ ↑(2 * k + 1))) = 0
Lean proof
/-- Every odd mode has zero quadrature, regardless of its degree (ATAP 19.3). -/
theorem clenshawCurtis_eval_T_odd {N : ℕ} (hN : N ≠ 0) (k : ℕ) :
clenshawCurtis N (fun x ↦ (T ℝ (↑(2 * k + 1 : ℕ))).eval x) = 0 := by
rw [clenshawCurtis_eval_T_alias hN]
have hp := aliasIndex_mod_two hN (2 * k + 1)
have he : aliasIndex N (2 * k + 1) = 2 * (aliasIndex N (2 * k + 1) / 2) + 1 := by omega
rw [he, integralMoment_odd]
theorem Polynomial.Chebyshev.abs_clenshawCurtis_eval_T_le_two {N : ℕ} (hN : N ≠ 0) (k : ℕ) :
|Polynomial.Chebyshev.clenshawCurtis N fun x => Polynomial.eval x (Polynomial.Chebyshev.T ℝ ↑k)| ≤ 2
Lean proof
/-- Positive weights of total mass two bound the quadrature of every mode. -/
theorem abs_clenshawCurtis_eval_T_le_two {N : ℕ} (hN : N ≠ 0) (k : ℕ) :
|clenshawCurtis N (fun x ↦ (T ℝ k).eval x)| ≤ 2 := by
have h := abs_clenshawCurtis_sub_le hN (f := fun x ↦ (T ℝ k).eval x)
(g := fun _ ↦ 0) (ε := 1) (fun j _ ↦ by
simpa using abs_eval_T_real_le_one (k : ℤ) (abs_le.mpr (node_mem_Icc (n := N) (i := j))))
simpa [clenshawCurtis_eq_sum_weights ⟨N, Nat.pos_of_ne_zero hN⟩] using h
theorem Polynomial.Chebyshev.abs_integralMoment_sub_clenshawCurtis_le {N k : ℕ} (hN : N ≠ 0) (hk : 4 ≤ k) :
|Polynomial.Chebyshev.integralMoment k -
Polynomial.Chebyshev.clenshawCurtis N fun x => Polynomial.eval x (Polynomial.Chebyshev.T ℝ ↑k)| ≤
32 / 15
Lean proof
/-- The modewise constant used in the sharper analytic quadrature bound. -/
theorem abs_integralMoment_sub_clenshawCurtis_le {N k : ℕ}
(hN : N ≠ 0) (hk : 4 ≤ k) :
|integralMoment k - clenshawCurtis N (fun x ↦ (T ℝ k).eval x)| ≤ 32 / 15 := by
exact (abs_sub _ _).trans ((add_le_add (abs_integralMoment_le_two_fifteenths hk)
(abs_clenshawCurtis_eval_T_le_two hN k)).trans_eq (by norm_num))