Home 计算机 数值计算

函数近似主要想解决的问题是，如果我们只知道一个函数的有限值，如何尽量精确地补齐这个函数的其他值
本节将主要介绍两种方法

  • Lagrange 插值
  • Spline 插值

暂时省略

  • Ermit 插值
  • Chebyshev 插值

# Lagrange 插值

作为数学事实，如果我们能唯一确定相异的 n+1n+1 个点 x0,x1,…,xnx_0, x_1, \ldots, x_n 上的函数值 f(x0),f(x1),…,f(xn)f(x_0), f(x_1), \ldots, f(x_n)，就可以唯一确定一个 nn 次多项式 PP 使得 P(xk)=f(xk)P(x_k) = f(x_k) 对于 k=0,1,…,nk = 0, 1, \ldots, n 都成立
该多项式称为函数 ff 在点 x0,x1,…,xnx_0, x_1, \ldots, x_n 处的 Lagrange 插值多项式

令多项式的形式为

p(x):=∑k=0nakxkp(x) := \sum_{k=0}^n a_k x^k

我们想要求解出系数 aka_k，根据条件，应该有 n+1n+1 个方程

{a0+a1x0+⋯+anx0n=f(x0)a0+a1x1+⋯+anx1n=f(x1)⋯a0+a1xn+⋯+anxnn=f(xn)\begin{cases} a_0 + a_1 x_0 + \cdots + a_n x_0^n = f(x_0) \\ a_0 + a_1 x_1 + \cdots + a_n x_1^n = f(x_1) \\ \cdots \\ a_0 + a_1 x_n + \cdots + a_n x_n^n = f(x_n) \end{cases}

将其改写为矩阵形式

(1x0⋯x0n1x1⋯x1n⋮⋮⋱⋮1xn⋯xnn)(a0a1⋮an)=(f(x0)f(x1)⋮f(xn))\begin{pmatrix} 1 & x_0 & \cdots & x_0^n \\ 1 & x_1 & \cdots & x_1^n \\ \vdots & \vdots & \ddots & \vdots\\ 1 & x_n & \cdots & x_n^n \end{pmatrix} \begin{pmatrix} a_0 \\ a_1 \\ \vdots \\ a_n \end{pmatrix} = \begin{pmatrix} f(x_0) \\ f(x_1) \\ \vdots \\ f(x_n) \end{pmatrix}

注意看系数矩阵，其行列式为 Vandermonde 行列式，因此

∣1x0⋯x0n1x1⋯x1n⋮⋮⋱⋮1xn⋯xnn∣=∏0≤j<k≤n(xk−xj)≠0\begin{vmatrix} 1 & x_0 & \cdots & x_0^n \\ 1 & x_1 & \cdots & x_1^n \\ \vdots & \vdots & \ddots & \vdots\\ 1 & x_n & \cdots & x_n^n \end{vmatrix} = \prod_{0 \leq j < k \leq n} (x_k - x_j) \neq 0

这说明系数矩阵可逆，因此方程组有唯一解，进而多项式 PP 也唯一存在

实际上，也可以直接写出多项式 PP 的表达式

p(x)=∑k=0nf(xk)∏j=0,j≠knx−xjxk−xjp(x) = \sum_{k=0}^n f(x_k) \prod_{j=0, j \neq k}^n \frac{x - x_j}{x_k - x_j}

接下来考察 Lagrange 插值的误差，令误差函数

E(x):=f(x)−p(x)E(x) := f(x) - p(x)

选取并固定一个点 xx，令

F(z):=E(z)−E(x)∏k=0n(z−xk)∏k=0n(x−xk)F(z) := E(z) - E(x) \frac{\prod_{k=0}^n (z - x_k)}{\prod_{k=0}^n (x - x_k)}

由于对于点 x0,x1,…,xnx_0, x_1, \ldots, x_n，Lagrange 插值多项式 pp 都与函数 ff 的值一致，因此 E(xk)=0E(x_k) = 0，并且对于所固定的 xx，F(x)=0F(x) = 0，也就是说 F(z)F(z) 在 n+2n+2 个点上取值为 00
对其重复应用 n+1n+1 次 Rolle 定理，可以得到 F(n+1)(ξ)=0F^{(n+1)}(\xi) = 0
再注意到 p(x)p(x) 是 nn 次多项式，因此 p(n+1)(x)=0p^{(n+1)}(x) = 0

综上所述

F(n+1)(ξ)=f(n+1)(ξ)−E(x)n!∏k=0n(x−xk)=0F^{(n+1)}(\xi) = f^{(n+1)}(\xi) - E(x) \frac{n!}{\prod_{k=0}^n (x - x_k)} = 0

移项得到

E(x)=f(n+1)(ξ)n!∏k=0n(x−xk)E(x) = \frac{f^{(n+1)}(\xi)}{n!} \prod_{k=0}^n (x - x_k)

也就是说，通过差值区间上函数的 n+1n+1 阶导数的最大值与差值区间长度的 n+1n+1 次幂的乘积，可以估计出 Lagrange 插值的误差

接下来考虑如何计算 Lagrange 插值多项式 pp 的值，直接套用其定义式会有 O(n2)O(n^2) 的时间复杂度，我们将引入 Neville 迭代法

我们想要通过 n+1n+1 个点来近似一个 nn 次多项式 p(x)p(x)
该算法的核心思想是:

  • 先移除最末尾的一个点，转为用 nn 个点近似 n−1n-1 次多项式的问题，并假设已经得到了这么一个多项式 q(x)q(x)
  • 再移除最前面的一个点，转为用 nn 个点近似 n−1n-1 次多项式的问题，并假设已经得到了这么一个多项式 r(x)r(x)

通过

p(x)=(x−x0)q(x)−(x−xn)r(x)xn−x0p(x) = \frac{(x - x_0) q(x) - (x - x_n) r(x)}{x_n - x_0}

就可以简单地得到 p(x)p(x) 的值

以下开始算法说明，给定 n+1n+1 个点 x0,x1,…,xnx_0, x_1, \ldots, x_n 与函数值 f(x0),f(x1),…,f(xn)f(x_0), f(x_1), \ldots, f(x_n)，以及一个点 α\alpha，我们想要计算 p(α)p(\alpha) 的值

  • 将 x0,x1,…,xnx_0, x_1, \ldots, x_n 按照靠近 α\alpha 的程度进行排序并重新编号
  • 对于每个 m=0,1,…,nm = 0, 1, \ldots, n，定义起始值 Pm,0:=f(xm)P_{m, 0} := f(x_m)
  • 逐一开始迭代，选取点 xmx_m，并在 ℓ=0,1,…,m\ell = 0, 1, \ldots, m 下循环计算

Pm,ℓ:=(α−xm−ℓ)Pm,ℓ−1−(α−xm)Pm−1,ℓ−1xm−xm−ℓP_{m, \ell} := \frac{(\alpha - x_{m-\ell}) P_{m, \ell-1} - (\alpha - x_m) P_{m-1, \ell-1}}{x_m - x_{m-\ell}}

  • 最后 Pn,nP_{n, n} 就是 p(α)p(\alpha) 的值

# Spline 插值

在点数较多的情况下，Lagrange 插值的误差会变得非常大，因此我们需要引入 Spline 插值来解决问题