多项式插值

已知函数 `f` 在 `n+1` 个不同的点处的值 `f(x_k) = y_k`, `quad k = 0, 1, cdots, n`, 求一个次数不超过 `n` 的多项式 `p(x)`, 满足 `p(x_k) = y_k`, `quad k = 0, 1, cdots, n`. `p` 称为 `f` 的 `n` 次插值多项式.

插值多项式的存在唯一性.

使用待定系数法. 设 `p` 的 `j` 次项系数为 `p_j`, 它们适合下面的线性方程组: `[1, x_0, cdots, x_0^n; 1, x_1, cdots, x_1^n; vdots, vdots, , vdots; 1, x_n, cdots, x_n^n] [p_0; p_1; vdots; p_n] = [y_0; y_1; vdots; y_n]`. 方程组的系数矩阵是 `n+1` 阶 Vandermonde 矩阵, 由假设 `x_k` 两两不同, 故系数矩阵非奇异, 方程组有唯一解.

下面将要介绍各种形式的插值多项式. 虽然形式不同, 由插值多项式的唯一性知道, 这些插值多项式都是相等的.

Lagrange 插值

线性插值 考虑两点 `(x_0, y_0)` 和 `(x_1, y_1)`, 它们的一次插值多项式显然就是连接这两点的直线, 方程为 `(y-y_0)/(x-x_0) = (y_1-y_0)/(x_1-x_0)`, 或者写成 `y = (x-x_1)/(x_0-x_1) y_0 + (x-x_0)/(x_1-x_0) y_1`. 插值多项式是两个基函数 `f_0(x) = (x-x_1)/(x_0-x_1)` 和 `f_1(x) = (x-x_0)/(x_1-x_0)` 的线性组合, 系数恰为 `y_0, y_1`. 容易验证 `f_0(x_0) = f_1(x_1) = 1`, `f_0(x_1) = f_1(x_0) = 0`.

Lagrange 插值 构造 Lagrange 基函数 `l_i(x) = prod_(j!=i) (x-x_j)/(x_i-x_j)`, 它满足 `l_i(x_k) = delta_(i k) = { 1, if k = i; 0, otherwise :}`. 换言之, `l_i` 是函数 `delta_(i k)` 的插值多项式. 一般地, 将函数写成 `delta_(i k)` 的线性组合: `f(x_k) = sum_k y_k delta_(i k)` `= sum_k y_k l_i(x_k)`. 故 `L(x) = sum_k y_k l_i(x)` 就是 `f` 的插值多项式, 称为 Lagrange 插值多项式.

Neville 插值 记函数在 `x_i, cdots, x_j` 上的插值多项式为 `p_(i j)`, 则有递推公式 `p_(i j)(x) = y_i`,
`p_(i j)(x) = (x-x_j)/(x_i-x_j) p_(i,j-1)(x) + (x-x_i)/(x_j-x_i) p_(i+1,j)(x)`. 换言之, `p_(i j)` 是两个低一次的多项式 `p_(i,j-1)` 和 `p_(i+1,j)` 线性插值的结果.

对插值多项式的次数 `j-i` 进行归纳. 次数为 0 时, 插值多项式是常数 `y_i`; 次数为 1 时, 即为线性插值. 现在假设 `p_(i,j-1)` 和 `p_(i+1,j)` 是相应点集上的插值多项式, 验证可知 `p_(i j)(x_k)` `= (x_k-x_j)/(x_i-x_j) p_(i,j-1)(x_k) + (x_k-x_i)/(x_j-x_i) p_(i+1,j)(x_k)` `= (x_k-x_j)/(x_i-x_j) y_k + (x_k-x_i)/(x_j-x_i) y_k` `= y_k`.

Newton 插值

差商 设 `x_i le x_(i+1) le cdots le x_j`, 函数 `f` 在该点集上的 `j-i` 阶差商定义为 `f[x_i] = f(x_i)`,
`f[x_i, cdots, x_j] = (f[x_(i+1), cdots, x_j] - f[x_i, cdots, x_(j-1)])/(x_j-x_i)`. 差商是离散版本的导数. 事实上, `lim_(x_j-x_i to 0) f[x_i, cdots, x_j]` `= (f^((j-i))(x_i))/(j-i)!`. 因此定义 `f overset(n+1)overbrace([x_i, cdots, x_i]) = (f^((n))(x_i))/n!`, 我们称这里 `x_i` 的重复度为 `n+1`.

设 `f(x) = (a-x)^-1`, 证明: `f[x_0, x_1, cdots, x_n] = prod_(k=0)^n (a-x_k)^-1`,

对 `n` 进行归纳证明. `n = 0` 时， `f[x_0] = (a-x_0)^-1`, 公式成立. 设公式对非负整数 `n` 成立, 则对 `n+1` 有 `f[x_0, cdots, x_(n+1)]`
`= (f[x_1, cdots, x_(n+1)] - f[x_0, cdots, x_n]) / (x_(n+1) - x_0)`
`= 1/(x_(n+1) - x_0) (prod_(k=1)^(n+1) (a-x_k)^-1 - prod_(k=0)^n (a-x_k)^-1)`
`= 1/(x_(n+1) - x_0) ((a-x_0) - (a-x_(n+1))) / (prod_(k=0)^(n+1) (a-x_k))`
`= prod_(k=0)^(n+1) (a-x_k)^-1`.

下面证明 Newton 插值公式, 这是 Taylor 公式的离散版本.

Newton 插值 `f` 在 `x_0, cdots, x_n` 上写为 `f(x) = N(x) + R_n(x)`, 其中 `N(x) = sum_(k=0)^n f[x_0, cdots, x_k] prod_(j=0)^(k-1) (x-x_j)`
`= f[x_0] + f[x_0, x_1](x-x_0) + cdots + f[x_0, cdots, x_n](x-x_0)cdots(x-x_(n-1))` 是 `f` 的 `n` 次插值多项式, `R_n(x) = f[x, x_0, cdots, x_n](x-x_0)cdots(x-x_n)` 是余项.

从余项入手, `R_k(x) = f[x, x_0, cdots, x_k](x-x_0)cdots(x-x_k)`
`= (f[x, x_0, cdots, x_(k-1)] - f[x_0, cdots, x_k]) (x-x_0) cdots(x-x_(k-1))`
`= R_(k-1)(x) - f[x_0, cdots, x_k] prod_(j=0)^(k-1) (x-x_j)`, 求和, `R_0(x) - R_n(x)` `= sum_(k=1)^n [R_(k-1)(x) - R_k(x)]` `= sum_(k=1)^n f[x_0, cdots, x_k] prod_(j=0)^(k-1) (x-x_j)`, 再代入 `R_0(x) = f(x) - f(x_0)` 即得结论. 最后, 注意到对 `k = 0, cdots, n` 有 `R_n(x_k) = 0`, 所以 `f(x_k) = N(x_k)`.

Newton 插值的特点是, 只有最后一项含有 `x_n`, 而前 `n` 项都与 `x_n` 无关. 因此我们可以逐点逐项地计算插值多项式. 要增加、删除最后一个节点, 只需增加、删除最后一项即可, 这是 Newton 插值带来的便利.

Hermite 插值

假设 `f(x) = ax^3 + bx^2 + cx + d`, 于是 `f'(x) = 3ax^2 + 2bx + c`. 规定 `f` 及其导数在给定点处的值, 运用待定系数法可以算出 `f` 的表达式:

`f(0)` `f(1)` `f'(0)` `f'(1)` `a` `b` `c` `d`
1 0 0 0 2 -3 0 1
0 1 0 0 -2 3 0 0
0 0 1 0 1 -2 1 0
0 0 0 1 1 -1 0 0

一般地, 若 `f(x_1) = y_1`, `f(x_2) = y_2`, `f'(x_1) = k_1`, `f'(x_2) = k_2`, 则通过线性叠加得到 `f(x) = y_1 f_1(t) + y_2 f_2(t) + k_1 f_3(t) + k_2 f_4(t)`, 其中, `x = a + (b-a) t`, `f_1` 到 `f_4` 是上面 4 个问题的解.

样条曲线

Bezier 曲线

区分：Bezier 曲线, Bessel 函数, Basel 问题.

Bezier 曲线 `n` 次 Bezier 曲线定义为 Bernstein 基函数 关于 `n+1` 个控制点 `P_0, P_1, cdots, P_n` 的线性组合 `P(t) = sum_k B_k^n(t) P_k`, `quad t in [0, 1]`. 其中 `B_k^n(t) = (n;k) t^k (1-t)^(n-k)`, `quad t in [0, 1]`. 容易看出 `sum_(k=0)^n B_k^n(t) = 1`. 且 `P(0) = P_0`, `P(1) = P_n`. 根据二项式系数的 Pascal 公式, Bernstein 基函数的递推关系为 `B_k^n(t) = (1-t) B_k^(n-1)(t) + t B_(k-1)^(n-1)(t)`, `quad 0 lt k lt n`.

一次 Bezier 曲线即线段 `P_0 P_1`: `P_(0, 1)(t)` `= (1-t) P_0 + t P_1` `= [t, 1][-1, 1; 1, 0][P_0; P_1]`. 二次 Bezier 曲线是 `P_0, P_2` 间的一条抛物线: `P_(0, 1, 2)(t)` `= (1-t)^2 P_0 + 2(1-t)t P_1 + t^2 P_2` `= [t^2, t, 1][1, -2, 1; -2, 2, 0; 1, 0, 0][P_0; P_1; P_2]`. 二次 Bezier 曲线还可以看作一次 Bezier 曲线的线性组合: `P_(0, 1, 2)(t) = (1-t)P_(0,1)(t) + t P_(1,2)(t)`. 一般地, `n` 次 Bezier 曲线是两条 `n-1` 次 Bezier 曲线的线性组合. 工程中最常用的是 2 次和 3 次 Bezier 曲线.

为什么二次 Bezier 曲线是抛物线? 考虑参数方程 `x = a t^2 + b t + c`, `quad y = d t^2 + e t + f` 注意到抛物线在仿射变换下的像仍为抛物线 (或退化为直线), 令 `u = d x - a y`, 则 `u` 是 `t` 的一次函数, 从而 `(x, t)` 是抛物线 `rArr (x, u)` 是抛物线 `rArr (x, y)` 是抛物线.

B 样条曲线

[来自 知乎@书剑飘零]

样条 (spline) 本来是指一些细小的有弹性的木条或者钢片, 用这些木条来设计模型结构. Bezier 曲线的不足之处是, 一个控制点的变动会影响整条曲线, "牵一发而动全身". B 样条曲线对这一点作了改进.

当第一个节点 `t_0` 和最后一个节点 `t_m` 的重复度为 `k+1` 时, `k` 次样条曲线的首尾分别与控制点的线段相切. 这种情形称为 clamped B-spline.

绘制样条曲线, 固定次数 `k = 3` 和节点列表 `[0, 0, 0, 0, 1, 2, 3, 3, 3, 3]`, 6 个控制点的位置可以用鼠标拖拽. 注意到拖拽一个控制点时, 只有曲线的一部分形状受到影响, 而每个控制点的影响范围取决于基函数的支撑集. 这就避免了 Bezier 曲线 "牵一发而动全身" 的问题.

假如次数 `k = 3`, 而节点列表改为 `[0, 0, 0, 0, 1, 1, 1, 1]`, 控制点的数目变成 4 个, 则上例退化为 3 次 Bezier 曲线. 因此 Bezier 曲线是 B 样条的一个特例.

样条曲线的计算

`P(t) = sum_(i=0)^n N_i^k(t) P_i` `= sum_(i=0)^n (t-t_i)/(t_(j-1)-t_i) N_i^(k-1)(t) P_i` `+ sum_(i=0)^n (t_j-t)/(t_j-t_(i+1)) N_(i+1)^(k-1)(t) P_i`, 限定 `t in [t_k, t_(n+1))` 时, 上述和式的第一项与最后一项为零, 得到 `P(t)` `= sum_(i=1)^n (t-t_i)/(t_(j-1)-t_i) N_i^(k-1)(t) P_i` `+ sum_(i=0)^(n-1) (t_j-t)/(t_j-t_(i+1)) N_(i+1)^(k-1)(t) P_i`
`= sum_(i=1)^n [(t-t_i)/(t_(j-1)-t_i) P_i + (t_(j-1)-t)/(t_(j-1)-t_i) P_(i-1)] N_i^(k-1)(t)`. 若定义 `tau_i^l = (t-t_i)/(t_(j-l)-t_i)`,
`P_i^l(t) = { P_i, if l = 0; tau_i^l P_i^(l-1)(t) + (1-tau_i^l) P_(i-1)^(l-1)(t), if l gt 0 :}`, 则 `P(t) = sum_(i=1)^n (tau_i^1 P_i^0(t) + (1-tau_i^1) P_(i-1)^0(t)) N_i^(k-1)(t)`
`= sum_(i=1)^n P_i^1(t) N_i^(k-1)(t)` `= sum_(i=2)^n P_i^2(t) N_i^(k-2)(t)` `= cdots` `= sum_(i=k)^n P_i^k(t) N_i^0(t)`. 上式右端是分段函数, `t in [t_i, t_(i+1))` 时, 有 `P(t) = P_i^k(t)`.

NURBS

[来自 知乎@EC果酱技术团队]

多项式逼近

Legendre 多项式

    Legendre 多项式
  1. 计算 `p_n(x) = "d"^n/dx^n (x^2-1)^n` 的 `n` 次项系数.
  2. 验证 `p_n(1) = 2^n n!`, `p_n(-1) = (-2)^n n!`.
  3. 验证 `{p_n}` 的正交性: `(p_n, p_m) = int_(-1)^1 p_n(x) p_m(x) dx` `= {0, if m != n; 2^(2n+1)/(2n+1) n!^2, otherwise:}`.
    4. 令 `L_n(x) = 1/(2^n n!) p_n(x)`, 则 `L_n(1) = 1`, `L_n(-1) = (-1)^n`, `(L_n, L_n) = 2//(2n+1)`. 称 `L_n(x)` 为 `n` 次 Legendre 多项式.
  4. 递推公式: `L_(n+1)(x) = (2n+1)/(n+1) x L_n(x) - n/(n+1) L_(n-1)(x)`.
`n` `L_n(x)`
0 `1`
1 `x`
2 `(3x^2 - 1)//2`
3 `(5x^3 - 3x)//2`
4 `(35x^4 - 30x^2 + 3)//8`
5 `(63x^5 - 70x^3 + 15x)//8`
6 `(231x^6 - 315x^4 + 105x^2 - 5)//16`
  1. `(x^2-1)^n` 是 `2n` 次多项式, 首项系数为 1, 因此它的 `n` 阶导数的首项系数是 `(2n) * (2n-1) cdots (n+1) = (2n)!//n!`.
  2. 以 `p_n(1)` 为例, `p_n(1) = {:"d"^n/dx^n|_(x=1) (x^2-1)^n` `= {:"d"^n/dy^n|_(y=0) (y(y+2))^n` `= (y(y+2))^n` 的 `n` 次项系数乘以 `n!` `= 2^n n!`.
  3. 不妨设 `m le n`, 分部积分得 `(p_n, p_m)` `= int_(-1)^1 "d"^n/dx^n (x^2-1)^n "d"^m/dx^m (x^2-1)^m dx` `= (-1)^n int_(-1)^1 (x^2-1)^n "d"^(m+n)/dx^(m+n) (x^2-1)^m dx`, 注意, 由于 `k lt n` 时, `{:"d"^k/dx^k|_(x=1) (x^2-1)^n` `= (y(y+2))^n` 的 `k` 次项系数乘以 `k!` `= 0`,
    `{:"d"^k/dx^k|_(x=-1) (x^2-1)^n` `= (y(y-2))^n` 的 `k` 次项系数乘以 `k!` `= 0`,     所以分部积分的中间项全部为零. 现在回到分部积分的结果, `(x^2-1)^m` 是 `2m` 次多项式, 我们有 `"d"^(m+n)/dx^(m+n) (x^2 - 1)^m` `= { (2n)!, if m = n; 0, if m lt n :}` 于是 `m != n` 时 `(p_n, p_m) = 0`. `m = n` 时 `(p_n, p_n)` `= (-1)^n (2n)! int_(-1)^1 (x^2 - 1)^n dx` `= 2 (2n)! int_0^1 (1-x^2)^n dx` `= 2 (2n)! int_0^(pi//2) cos^(2n+1) theta "d"theta` `= 2 (2n)! ((2n)!!)/((2n+1)!!)` `= 2^(2n+1)/(2n+1) n!^2`.

Chebyshev 多项式

下面主要讨论 Chebyshev 多项式. 它来源于三角函数的 `n` 倍角问题, 在数值分析中作为一类 "最佳多项式" 有着重要应用. 使用和差化积公式, 便得到求 `n` 倍角正余弦值的 Chebyshev 算法: `cos n theta = 2 cos theta cos(n-1)theta - cos(n-2)theta`,
`sin n theta = 2 cos theta sin(n-1)theta - sin(n-2)theta`. 而初始条件为 `cos 0 theta = 1`, `quad cos 1 theta = cos theta`,
`sin 0 theta = 0`, `quad sin 1 theta = sin theta`.

容易看出 `cos n theta` 可以被表为 `cos theta` 的多项式, 称为 (第一类) Chebyshev 多项式: `cos n theta = T_n(cos theta)`. 类似地, `U_n` 称为 第二类 Chebyshev 多项式: `(sin n theta)/(sin theta) = U_(n-1)(cos theta)`.

第一类 Chebyshev 多项式的递归定义: `T_n(x) = { 1, if n = 0; x, if n = 1; 2x T_(n-1)(x) - T_(n-2)(x), if n ge 2; :}` 第二类 Chebyshev 多项式的递归定义: `U_(n-1)(x) = { 0, if n = 0; 1, if n = 1; 2x U_(n-2)(x) - U_(n-3)(x), if n ge 2; :}`

显式表达式: `T_n = [(x + sqrt(1-x^2)"i")^n + (x - sqrt(1-x^2))^n]//2`, `quad x in [-1, 1]`. 代入 `x = cos theta` 即可验证.

Chebyshev 多项式
`n` `T_n(x)` `U_(n-1)(x)`
`0` `1` `0`
`1` `x` `1`
`2` `2x^2-1` `2x`
`3` `4x^3 - 3x` `4x^2-1`
`4` `8x^4 - 8x^2 + 1` `8x^3 - 4x`
`5` `16x^5 - 20x^3 + 5x` `16x^4 - 12x^2 + 1`

零点与最值 由公式 `cos n theta = T_n(cos theta)` 知道, `T_n(x)` 在 `x in [-1, 1]` 上有 `n` 个零点: `cos((2k-1)/(2n) pi)`, `quad k = 1, cdots, n`. 且 `T_n(x)` 在 `x in [-1, 1]` 上的最值在 `cos(k/n pi)`, `quad k = 0, cdots, n` 处取得, 其中 `k` 为偶数时取得最大值 `1`, 奇数时取得最小值 `-1`.

    最佳一致逼近 在 `[-1, 1]` 上全体首 1 的 `n` 次实系数多项式中, `T_n(x) * 2^(1-n)` 是对 `0` 的最佳一致逼近. 换言之, 对任意首 1 的 `n` 次实系数多项式 `P_n(x)`, 都有 `max_(x in [-1, 1]) |P_n(x)|` `ge max_(x in [-1, 1]) |T_n(x)| * 2^(1-n)` `= 2^(1-n)`.

反设 `max_(x in [-1, 1]) |P_n(x)| lt 2^(1-n)`. 考虑 `T_n(x)` 的最值点 `x_k = cos(k/n pi)`, `quad k = 0, cdots, n`, 定义 `g(x) = P_n(x) - T_n(x) * 2^(1-n)`, 于是 `g(x_k) { gt 0, if k " odd"; lt 0, if k " even" :}` `g(x)` 的值在 `x_k`, `k = 0, cdots, n` 共 `n+1` 个点处交错, 故 `g(x)` 在 `[-1, 1]` 上有 `n` 个零点. 但 `g(x)` 是次数不超过 `n` 的多项式, 因此必有 `g(x) -= 0`, 即 `P_n(x) = T_n(x) * 2^(1-n)`, 一个矛盾.

最大系数性质 若 `f(x)` 是 `n` 次实系数多项式, 且 `max_(x in [-1, 1])|f(x)| le 1`, 则 `f` 的 `n-2k` 次项系数的绝对值不超过 `T_n(x)` 的对应系数的绝对值.
例如, 取 `n = 2`, 则 `T_n(x) = 2x^2 - 1` 在 `x = cos(k/2 pi)` 处取得最值, 即 `T_n(1) = T_n(-1) = 1`, `T_n(0) = -1`. 对于任意二次实系数多项式 `f(x) = a x^2 + b x + c`, 由于 `max_(x in [-1,1])|f(x)| le 1`, 特别地 `f(1) = a + b + c le 1`,
`f(-1) = a - b + c le 1`,
`f(0) = c ge -1`. 前两式相加得 `a + c le 1`, 再由 `c ge -1` 得 `a le 2`. 类似可得 `|c| le 1`, `|a| le 2`.

??