#include <boost/math/special_functions/chebyshev.hpp>
namespace boost{ namespace math{ template<class Real1, class Real2, class Real3> calculated-result-type chebyshev_next(Real1 const & x, Real2 const & Tn, Real3 const & Tn_1); template<class Real> calculated-result-type chebyshev_t(unsigned n, Real const & x); template<class Real, class Policy> calculated-result-type chebyshev_t(unsigned n, Real const & x, const Policy&); template<class Real> calculated-result-type chebyshev_u(unsigned n, Real const & x); template<class Real, class Policy> calculated-result-type chebyshev_u(unsigned n, Real const & x, const Policy&); template<class Real> calculated-result-type chebyshev_t_prime(unsigned n, Real const & x); template<class Real1, class Real2> calculated-result-type chebyshev_clenshaw_recurrence(const Real1* const c, size_t length, Real2 x); template<typename Real> Real chebyshev_clenshaw_recurrence(const Real* const c, size_t length, Real a, Real b, Real x); }} // namespaces
“实分析学家离不开傅里叶分析,复分析学家离不开洛朗分析,数值分析学家离不开切比雪夫分析” --Lloyd N. Trefethen
第一类切比雪夫多项式由以下递推关系定义:Tn+1(x) := 2xTn(x) - Tn-1(x), n > 0,其中 T0(x) := 1 且 T1(x) := x。这些可以使用以下简单的 Boost 代码计算
double x = 0.5; double T12 = boost::math::chebyshev_t(12, x);
导数的计算也很直接
double T12_prime = boost::math::chebyshev_t_prime(12, x);
通过这些函数评估第 n 阶切比雪夫多项式的复杂度是线性的。因此,它们不适合用于计算(例如)切比雪夫级数,因为线性缩放函数的总和呈二次方缩放。尽管有非常复杂的切比雪夫级数评估算法,但下面介绍了一种线性时间算法
double x = 0.5; std::vector<double> c{14.2, -13.7, 82.3, 96}; double T0 = 1; double T1 = x; double f = c[0]*T0/2; unsigned l = 1; while(l < c.size()) { f += c[l]*T1; std::swap(T0, T1); T1 = boost::math::chebyshev_next(x, T0, T1); ++l; }
这使用 chebyshev_next
函数在恒定时间内评估切比雪夫级数的每一项。然而,当 x 接近 1 时,这种朴素的算法会造成灾难性的精度损失。Clenshaw 给出了一种缓解这种情况的方法,并在 boost 中实现为
double x = 0.5; std::vector<double> c{14.2, -13.7, 82.3, 96}; double f = chebyshev_clenshaw_recurrence(c.data(), c.size(), x);
注意:在切比雪夫级数中,我们对第一个系数的定义与 Clenshaw 的原始工作存在 2 倍的差异。这是因为切比雪夫级数展开的符号存在两种传统,
和
boost math 始终使用第二种约定,第一个系数的因子为 1/2。
另一种常见的用例是必须在某个区间 [a, b] 上评估多项式时。转换为区间 [-1, 1] 会导致一些精度问题,同时也为我们提供了一些机会。对于这种情况,我们使用 Reinch 对 Clenshaw 递推的修改,这在 此处 也有讨论。用法如下
double x = 9; double a = 7; double b = 12; std::vector<double> c{14.2, -13.7, 82.3, 96}; double f = chebyshev_clenshaw_recurrence(c.data(), c.size(), a, b, x);
第二类切比雪夫多项式可以通过 chebyshev_u
评估
double x = -0.23; double U1 = boost::math::chebyshev_u(1, x);
已知通过三项递推关系评估切比雪夫多项式在 x ∊ [-1, 1] 时是 混合向前向后稳定的。然而,作者不知道 x 在 [-1, 1] 之外是否有类似的结果。因此,强烈建议不要在 [-1, 1] 之外评估切比雪夫多项式。也就是说,计算过程中的小舍入误差通常会导致这种情况,并且由于这些小问题而终止计算是非常令人沮丧的。因此,chebyshev_t
和 chebyshev_u
具有用于 x > 1 和 x < -1 的代码路径,这些代码路径不使用三项递推关系。这些代码路径 慢得多,应尽可能避免使用。
评估切比雪夫级数相对简单。真正的挑战是切比雪夫级数的 生成。为此,boost 提供了一个 切比雪夫变换,这是一种投影算子,可将函数投影到切比雪夫多项式的有限维跨度上。但在我们讨论 API 之前,让我们分析一下为什么我们可能想要将函数投影到切比雪夫多项式的跨度上。
API 如下所示。
#include <boost/math/special_functions/chebyshev_transform.hpp>
namespace boost{ namespace math{ template<class Real> class chebyshev_transform { public: template<class F> chebyshev_transform(const F& f, Real a, Real b, Real tol=500*std::numeric_limits<Real>::epsilon()); Real operator()(Real x) const Real integrate() const const std::vector<Real>& coefficients() const Real prime(Real x) const }; }}// end namespaces
切比雪夫变换接受一个函数 f,并返回 f 的 近极小极大 近似值,以切比雪夫多项式表示。通过 近极小极大,我们的意思是生成的切比雪夫多项式“非常接近”多项式 pn,后者最小化 f - pn 的一致范数。“非常接近”的概念可以变得严谨;有关详细信息,请参阅 Trefethen 的《逼近理论与逼近实践》。
切比雪夫变换的工作原理是:通过在切比雪夫点评估输入函数来创建值向量,然后在结果向量上执行离散余弦变换。为了高效地执行此操作,我们使用了 FFTW3。因此,要编译,您必须安装 FFTW3
,并链接 -lfftw3
以获得双精度,-lfftw3f
以获得单精度,-lfftw3l
以获得长双精度,以及 -lfftw3q
以获得四倍 (__float128
) 精度。在已知切比雪夫级数的系数后,例程会重新遍历它们,并滤除所有绝对比率小于构造函数中请求的容差的最大系数的系数。