跳转至

Chap 3: Interpolation and Polynomial Approximation⚓︎

5001 个字 21 行代码 预计阅读时间 25 分钟

如果函数 \(y = f(x)\) 的计算过于复杂,或者甚至是未知的,一种近似求解的办法是首先在一组点序列 \(x_0, \dots, x_n\) 上获取一组函数值 \(y_0 = f(x_0), \dots, y_n = f(x_n)\),然后根据这些值构造一个相对简单的近似函数 \(g(x) \approx f(x)\)

如果 \(g(x)\) 满足 \(\forall\ i = 0, \dots, n, g(x_i) = f(x_i)\),我们称 \(g(x)\) \(f(x)\) 插值函数(interpolating function)。最常用的插值函数是代数多项式(algebraic polynomials)。

Interpolation and the Lagrange Polynomial⚓︎

目标:找到 \(n\) 阶多项式 \(P_n(x) = a_0 + a_1 x + \dots + a_n x^n\),使得 \(\forall\ i = 0, \dots, n, P_n(x_i) = y_i\)

注:对任何 \(i \ne j\),必须满足 \(x_i \ne x_j\)

\(n = 1\) 时:给定 \(x_0, x_1; y_0, y_1\)。找到 \(P_1(x) = a_0 + a_1 x\),使得 \(P_1(x_0) = y_0, P_1(x_1) = y_1\)。此时 \(P_1(x)\) 是一个经过给定两点 \((x_0, y_0), (x_1, y_1)\) 线函数(line function),即:

\[ \begin{align} P_1(x) & = y_0 + \dfrac{y_1 - y_0}{x_1 - x_0}(x - x_0) \notag \\ & = \underbrace{\Big(\dfrac{x - x_1}{x_0 - x_1}\Big)}_{\textcolor{red}{L_{1, 0}}(x)} y_0 + \underbrace{\Big(\dfrac{x - x_0}{x_1 - x_0}\Big)}_{\textcolor{red}{L_{1, 1}}(x)} y_0 = \sum\limits_{i=0}^1 \textcolor{red}{L_{1, i}}(x)y_i \notag \end{align} \]

其中标红的项被称为拉格朗日基(Lagrange Basis),它满足:\(L_{1, i}(x_j) = \delta_{ij}\)(称为 Kronecker 增量)


\(n \ge 1\) 时:我们要寻找 \(L_{n, i}(x)\ (i = 0, \dots, n)\),使得 \(L_{n, i} (x_j) = \delta_{ij}\)。然后令 \(P_n(x) = \sum\limits_{i=0}^n L_{n, i}(x) y_i\)。因此 \(P_n(x_i) = y_i\)

每个 \(L_{n, i}\) 都有 \(n\) 个根 \(x_0 \dots \widehat{x_i} \dots x_n\)。可以得到:

\[ L_{n, j}(x) = C_i (x - x_0) \dots (\widehat{x - x_i}) \dots (x - x_n) = C_i \prod\limits_{\substack{j \ne i \\ j = 0}}^n (x - x_j) \]

\(L_{n, i}(x_i) = 1\) 时,\(C_i = \prod\limits_{j \ne i} \dfrac{1}{x_i - x_j}\)

\[ L_{n, i}(x) = \prod\limits_{\substack{j \ne i \\ j = 0}}^n \dfrac{(x - x_j)}{(x_i - x_j)} \Rightarrow P_n(x) = \sum\limits_{i=0}^n L_{n, i}(x) y_i \]

这里的 \(P_n(x)\) 就是 n 拉格朗日插值多项式(n-th Lagrange interpolating polynomial)。

定理

如果 \(x_0, x_1, \dots, x_n\) \(n + 1\) 个不同的数,且用函数 \(f\) 得到这 \(n + 1\) 数对应的函数值,那么第 n 拉格朗日插值多项式就是唯一的

证明

(反证法)如果不是唯一的没那么存在两个多项式 \(P_n(x)\) \(Q_n(x)\),它们都满足插值条件。这样我们可以得到一个多项式 \(D(x) = P_n(x) - Q_n(x)\),它的阶数 \(\le n\)。但 \(D(x)\) \(n + 1\) 个不同的根 \(x_0, x_1, \dots, x_n\),这和它的阶数条件矛盾。因此假设不成立,这样就证明了多项式的唯一性。

插值多项式不是唯一的,除非它的阶数被约束在不超过 \(n\) 的范围内。

相应的反例就是 \(P(x) = L_n(x) + p(x) \prod\limits_{i=0}^n (x - x_i)\),其中 \(p(x)\) 可以是任意阶数的多项式。


下面我们来分析余项 (remainder):假如 \(a \le x_0 < x_1 < \dots < x_n \le b\) \(f \in C^{n+1} [a, b]\)。考虑截断误差 \(R_n(x) = f(x) - P_n(x)\)

罗尔定理

如果 \(\varphi(x)\) 足够平滑,且 \(\varphi(x_0) = \varphi(x_1) = 0\),那么 \(\exists \xi \in (x_0, x_1)\),使得 \(\varphi'(\xi) = 0\)

一般来说,如果 \(\varphi(x_0) = \varphi(x_1) \varphi(x_2) = 0\),使得 \(\varphi'(\xi_0) = \varphi'(\xi_1) = 0\),那么 \(\exists \xi \in (\xi_0, \xi_1)\),使得 \(\varphi''(\xi) = 0\)

同理,\(\varphi(x_0) = \dots = \varphi(x_n) = 0 \Rightarrow \exists \xi \in (a, b)\),使得 \(\xi^{(n)}(\xi) = 0\)

\(R_n(x)\) 至少有 \(n + 1\) 个根 \(\Rightarrow R_n(x) = K(x) \prod\limits_{i=0}^n (x - x_i)\)

修正任何 \(x \ne x_i\ (i = 0, \dots, n)\)。定义为 \(t \in [a, b]\) 定义函数 \(g\) 为:

\[ g(t) = R_n(t) - K(x) \prod\limits_{i=0}^n (t - x_i) \]

\(g(x)\) \(n + 2\) 个不同的根 \(x_0, \dots, x_n, x\ \Rightarrow g^{(n+1)}(\xi_x) = 0, \xi_x \in (a, b)\)

\[ \begin{align} & f^{(n+1)}(\xi_x) - \cancel{P_n^{(n+1)}(\xi_x)} - K(x)(n+1)! = R_n^{(n+1)}(\xi_x) - K(x)(n+1)! \notag \\ & \Rightarrow K(x) = \dfrac{f^{(n+1)}(\xi_x)}{(n+1)!} \notag \end{align} \]

因此 \(R(x) = \dfrac{f^{(n+1)}(\xi_x)}{(n+1)!} \prod\limits_{i=0}^n (x - x_i)\)

  • 因为大多数情况下 \(\xi_x\) 无法确定,因此我们会获取 \(f^{(n+1)}\) 的上界。也就是说,得到一个估计值 \(M_{n+1}\),满足 \(\forall\ x \in (a, b),\ |f^{(n+1)}(x)| \le M_{n+1}\),并且取 \(\dfrac{M_{n+1}}{(n+1)!} \prod\limits_{i=0}^n |x - x_i|\) 为总误差的上界。
  • 对于任意阶数 \(\le n\) 的多项式函数 \(f\),插值多项式是精确的,因为 \(f^{(n+1)}(x) \equiv 0\)
Quiz

A

讨论

假设为 \(f(x) = e^x\ (x \in [0, 1])\) 做一个表格,表格中的每一项精确到 8 位小数,并且步幅为 \(h\)。对于绝对误差不超过 \(10^{-6}\) 的线性插值,\(h\) 应该是多少?

注意

答案由 Gemini 2.5 Pro 生成。

我们需要确定步长 \(h\) 的值,使得对于函数 \(f(x) = e^x\) 在区间 \([0, 1]\) 上的线性插值,其绝对误差不超过 \(10^{-6}\)。表格条目的精度为小数点后 8 位。

1. 线性插值总误差

线性插值的总误差 \(E_{total}\) 由两部分组成:截断误差 \(E_{trunc}\) 和由于表格值舍入误差引起的插值误差 \(E_{roundoff}\)

\(|E_{total}| \le |E_{trunc}| + |E_{roundoff}|\)

2. 截断误差

线性插值的截断误差公式为:

\(|E_{trunc}(x)| = \left| \frac{f''(\xi)}{2!} (x-x_i)(x-x_{i+1}) \right|\)

其中 \(x_i \le x \le x_{i+1}\)\(x_{i+1} - x_i = h\),且 \(\xi \in (x_i, x_{i+1})\)

该误差项中的 \(|(x-x_i)(x-x_{i+1})|\) 在区间中点 \(x = x_i + h/2\) 达到最大值 \(h^2/4\)

所以,\(|E_{trunc}| \le \frac{h^2}{8} \max_{\xi \in [0,1]} |f''(\xi)|\)

对于 \(f(x) = e^x\)

\(f'(x) = e^x\)

\(f''(x) = e^x\)

在区间 \([0, 1]\) 上,\(f''(x) = e^x\) 是增函数,其最大值为 \(f''(1) = e^1 = e\)

因此,\(|E_{trunc}| \le \frac{h^2 e}{8}\)

3. 舍入误差

表格条目的精度为小数点后 8 位,这意味着每个列表值 \(f(x_i)\) 的误差 \(|\epsilon_i| \le 0.5 \times 10^{-8}\)

线性插值公式为 \(P_1(x) = L_i(x)f(x_i) + L_{i+1}(x)f(x_{i+1})\),其中 \(L_i(x) = \frac{x_{i+1}-x}{h}\) \(L_{i+1}(x) = \frac{x-x_i}{h}\) 是拉格朗日基函数。

由于 \(L_i(x) \ge 0\), \(L_{i+1}(x) \ge 0\) \(L_i(x) + L_{i+1}(x) = 1\) 对于 \(x \in [x_i, x_{i+1}]\), 由表格值不精确引起的插值误差 \(|E_{roundoff}|\) 的上界为:

\(|E_{roundoff}| \le L_i(x)|\epsilon_i| + L_{i+1}(x)|\epsilon_{i+1}| \le (L_i(x) + L_{i+1}(x)) \times (0.5 \times 10^{-8}) = 1 \times 0.5 \times 10^{-8} = 0.5 \times 10^{-8}\)

4. 确定步长 \(h\)

总误差必须满足:

\(\frac{h^2 e}{8} + 0.5 \times 10^{-8} \le 10^{-6}\)

\(\frac{h^2 e}{8} \le 10^{-6} - 0.5 \times 10^{-8}\)

\(\frac{h^2 e}{8} \le 10^{-6} - 0.005 \times 10^{-6}\)

\(\frac{h^2 e}{8} \le 0.995 \times 10^{-6}\)

\(h^2 \le \frac{8 \times 0.995 \times 10^{-6}}{e}\)

使用 \(e \approx 2.71828\)

\(h^2 \le \frac{7.96}{e} \times 10^{-6} \approx \frac{7.96}{2.71828} \times 10^{-6} \approx 2.92826 \times 10^{-6}\)

\(h \le \sqrt{2.92826 \times 10^{-6}} \approx 1.7112158 \times 10^{-3}\)

所以,步长 \(h\) 应满足:

\(h \le 0.0017112\)

为了确保线性插值的绝对误差最多为 \(10^{-6}\),步长 \(h\) 不应超过 \(0.0017112\)

Neville's Method⚓︎

例子

给定 \(\sin \dfrac{\pi}{6} = \dfrac{1}{2}, \sin \dfrac{\pi}{4} = \dfrac{1}{\sqrt{2}}, \sin \dfrac{\pi}{3} = \dfrac{\sqrt{3}}{2}\)。使用关于 \(\sin x\) 的线形和二次拉格朗日多项式,计算 \(\sin 50 \degree\) 并评估误差。

(已知 \(\sin 50 \degree = 0.7660444...\)

  • 先使用 \(x_0, x_1\) \(x_1, x_2\) 计算线形插值。

    • 使用 \(x_0 = \dfrac{\pi}{6}, x_1 = \dfrac{\pi}{4}\)
      • \(P_1(x) = \dfrac{x - \frac{\pi}{4}}{\frac{\pi}{6} - \frac{\pi}{4}} \times \dfrac{1}{2} + \dfrac{x - \frac{\pi}{6}}{\frac{\pi}{4} - \frac{\pi}{6}} \times \dfrac{1}{\sqrt{2}}\)
      • \(50 \degree = \dfrac{5\pi}{18}\)
      • \(\sin 50 \degree \approx P_1(\dfrac{5 \pi}{18}) \approx 0.77614\)
      • \(f(x) = \sin x, f^{(2)} = - sin \xi_x, \xi_x \in (\dfrac{\pi}{6}), \dfrac{\pi}{3}\),且 \(\dfrac{1}{2} < \sin \xi_x < \dfrac{\sqrt{3}}{2}\)
      • \(R_1 (x) = \dfrac{f^{(2)(\xi_x)}}{2!}(x - \dfrac{\pi}{6})(x - \dfrac{\pi}{4})\),得到 \(-0.01319 < R_1(\dfrac{5\pi}{18}) < -0.00762\),因此外推误差 \(\approx -0.01001\)
    • 使用 \(x_1 = \dfrac{\pi}{4}, x_2 = \dfrac{\pi}{3}\)
      • 计算得到 \(\sin 50 \degree \approx 0.76008, 0.00538 < \widetilde{R_1}(\dfrac{5\pi}{18}) < 0.00660\)
      • 因此插值误差 \(\approx 0.00596\)

        一般情况下,插值比外推效果更好。

  • 再使用 \(x_0, x_1, x_2\) 计算二次插值。

    • \(P_2(x) = \frac{(x - \frac{\pi}{4})(x - \frac{\pi}{3})}{(\frac{\pi}{6} - \frac{\pi}{4})(\frac{\pi}{6} - \frac{\pi}{2})} \times \dfrac{1}{2} + \frac{(x - \frac{\pi}{6})(x - \frac{\pi}{3})}{(\frac{\pi}{4} - \frac{\pi}{6})(\frac{\pi}{4} - \frac{\pi}{3})} \times \dfrac{1}{\sqrt{2}} + \frac{(x - \frac{\pi}{6})(x - \frac{\pi}{4})}{(\frac{\pi}{3} - \frac{\pi}{6})(\frac{\pi}{3} - \frac{\pi}{4})} \times \dfrac{\sqrt{3}}{2}\)
    • \(\sin 50 \degree \approx P_2(\dfrac{5\pi}{18}) \approx 0.76543\)
    • \(R_2(x) = \dfrac{- \cos \xi_x}{3!}(x - \dfrac{\pi}{6})(x - \dfrac{\pi}{4})(x - \dfrac{\pi}{3}),\ \dfrac{1}{2} < \cos \xi_x < \dfrac{\sqrt{3}}{2}\)
    • \(0.00044 < R_2(\dfrac{5 \pi}{18}) < 0.00077\),所以二次插值的误差 \(\approx 0.00061\)

更高次的插值法通常会带来更好的结果,但并不总是如此。

定义

\(f\) 是关于 \(x_0, x_1, \dots, x_n\) 的函数,并假设 \(m_1, \dots, m_k\) \(k\) 个不同的整数且满足 \(\forall i, 0 \le m_i \le n\)。拉格朗日多项式在 \(k\) 个点 \(x_{m_1}, \dots, x_{m_k}\) 上与 \(f(x)\) 具有相同值时,记作 \(P_{m_1, \dots, m_k}(x)\)

定理

\(f\) 是关于 \(x_0, x_1, \dots, x_k\) 的函数,并令 \(x_i, x_j\) 为其中两个不相等的数,那么:

\[ P(x) = \dfrac{(x - x_j)P_{0, 1, \dots, j-1, j+1, \dots, k}(x) - (x - x_i)P_{0, 1, \dots, i-1, i+1, \dots, k}(x)}{x_i - x_j} \]

描述了在 \(k+1\) 个点 \(x_0, x_1, \dots, x_k\) 上向 \(f\) 插值的第 \(k\) 个拉格朗日多项式。

证明
  • 对任意 \(0 \le r \le k\) \(r \ne i \text{ and } j\),两个在分子上插值多项式等于 \(f(x_r)\),因此 \(P(x_r) = f(x_r)\)
  • 第一个在分子上的多项式等于 \(f(x_i)\),且第二项为 0,所以 \(P(x_i) = f(x_i)\)。同理,\(P(x_j) = f(x_j)\)
  • 因此,在 \(k+1\) 个点 \(x_0, x_1, \dots, x_k\) 上向 \(f\) 插值的第 \(k\) 个拉格朗日多项式是唯一的

上述定理表明插值多项式可以递归生成。比如,它们可以以下表所示的方式,一行行地生成插值多项式:

\[ \begin{matrix} x_0 & P_0 & & & & \\ x_1 & P_1 & P_{0, 1} & & & \\ x_2 & P_2 & P_{1, 2} & P_{0, 1, 2} & & \\ x_3 & P_3 & P_{2, 3} & P_{1, 2, 3} & P_{0, 1, 2, 3} \\ x_4 & P_4 & P_{3, 4} & P_{2, 3, 4} & P_{1, 2, 3, 4} & P_{0, 1, 2, 3, 4} \end{matrix} \]

上述过程被称为 Neville 。但 \(P\) 的记号显得过于笨重(一堆下标表示参与到多项式中的插值点。观察发现,只需要两个下标就行了——我们用新的记号 \(Q_{i, j}(x)\ (0 \le j \le i)\) 来表示阶数为 \(j\) 的,在 \((j+1)\) 个数 \(x_{i-j}, x_{i-j+1}, \dots, x_{i-1}, x_i\) 上的插值多项式,即:\(Q_{i, j} = P_{i-j, i-j+1, \dots, i-1, i}\)

算法:Neville 迭代插值

求解对于数 \(x\),在 \(n+1\) 个不同的数 \(x_0, \dots, x_n\) 上的函数 \(f\) 的插值多项式 \(P\)

  • 输入:数 \(x_0, x_1, \dots, x_n\);值 \(f(x_0), f(x_1), \dots, f(x_n)\),分别作为 \(Q\) 的第 1 \(Q_{0, 0}, Q_{1, 0}, \dots, Q_{n, 0}\) 上的值。
  • 输出:表 \(Q\),其中 \(P(x) = Q_{n, n}\)
Step 1  for i = 1, 2, ..., n:
            for j = 1, 2, ..., i:
                set Q[i][j] = ((x - x[i-j]) * Q[i][j-1] - (x - x[i]) * Q[i-1][j-1]) / (x[i] - x[i-j]);
Step 2  Output(Q);
        STOP;

Divided Difference⚓︎

Divided differences is a recursive division process. Given a sequence of data points \((x_0,y_0),\dots,(x_n,y_n)\), the method calculates the coefficients of the interpolation polynomial of these points in the Newton form. -- Wikipedia

  • 1 阶差商\(f[x_i, x_j] = \dfrac{f(x_i) - f(x_j)}{x_i - x_j} (i \ne j, x_i \ne x_j)\)
  • 2 阶差商\(f[x_i, x_j, x_k] = \dfrac{f[x_i, x_j] - f[x_j, x_k]}{x_i - x_k} (i \ne k)\)
  • \(k+1\) 阶差商

    \[ \begin{align} f[x_0, \dots, x_{k+1}] & = \dfrac{f[\textcolor{cornflowerblue}{x_0}, x_1, \dots, x_k] - f[x_1, \dots, x_k, \textcolor{cornflowerblue}{x_{k+1}}]}{\textcolor{cornflowerblue}{x_0 - x_{k+1}}} \notag \\ & = \dfrac{f[x_0, \dots, x_{k-1}, \textcolor{cornflowerblue}{x_k}] - f[x_0, \dots, x_{k-1}, \textcolor{cornflowerblue}{x_{k+1}}]}{\textcolor{cornflowerblue}{x_k - x_{k+1}}} \notag \end{align} \]

事实上,\(f[x_0, \dots, x_k] = \sum\limits_{i=0}^k \dfrac{f(x_i)}{\omega_{k+1}' (x_i)}\),其中 \(\omega_{k+1}(x) = \prod\limits_{i=0}^k (x - x_i), \omega_{k+1}'(x_i) = \prod\limits_{\substack{j = 0 \\ j \ne i}}^k (x_i - x_j)\)。这个公式的要点在于:\(f[x_0, \dots, x_k]\) 的值和 \(x_0, \dots, x_k\) 的顺序无关。

Newton's Interpolation⚓︎

目标:得到 \(N_n(x) = a_0 + a_1(x - x_0) + a_2(x - x_0)(x - x_1) + \dots + a_n(x - x_0) \dots (x - x_{n-1})\)

\[ \begin{cases} f(x) = f(x_0) + (x - x_0)f[x, x_0] & (1)\\ f[x, x_0] = f[x_0, x_1] + (x - x_1)f[x, x_0, x_1] & (2)\\ \dots\ \dots\ \dots\\ f[x, x_0, \dots, x_{n-1}] = f[x_0, \dots, x_n] + (x - x_n) f[x, x_0, \dots, x_n] & (n-1) \end{cases} \]

计算 \((1) + (x - x_0) \times (2) + \dots + (x - x_0) \dots (x - x_{n-1}) \times (n-1)\),得到:

\[ \begin{align} f(x) = & \textcolor{red}{f(x_0) + f[x_0, x_1](x - x_0) + f[x_0, x_1, x_2](x - x_0)(x - x_1) + \dots} \notag \\ & \textcolor{red}{+ f[x_0, \dots, x_n](x - x_0) \dots (x - x_{n-1})} \notag \\ & \textcolor{green}{+ f[x, x_0, \dots, x_n](x - x_0) \dots (x - x_{n-1})(x - x_n)} \notag \end{align} \]

其中红色部分就是我们要求的 \(N_n(x)\) ,而绿色部分是 \(R_n(x)\)。另外,\(a_i = f[x_0, \dots, x_i]\)

算法:牛顿插值差商公式

求得对于数 \(x\),在 \(n+1\) 个不同的数 \(x_0, \dots, x_n\) 上的函数 \(f\) 的插值多项式 \(P\) 的差商系数。

  • 输入:数 \(x_0, x_1, \dots, x_n\);值 \(f(x_0), f(x_1), \dots, f(x_n)\),分别记作 \(F_{0, 0}, F_{1, 0}, \dots, F_{n,0}\)
  • 输出:数 \(F_{0, 0}, F_{1, 1}, F_{n, n}\),其中 \(P(x) = \sum\limits_{i=0}^n F_{i, i} \prod_{j=0}^{i-1} (x - x_j)\)
Step 1  for i = 1, 2, ..., n:
            for j = 1, 2, ..., i:
                set F[i][j] = (F[i][j-1] - F[i-1][j-1]) / (x[i] - x[i-j]);
Step 2  Output(F[0][0], F[1][1], ..., F[n][n]);  // F[i][i] is f[x[0], x[1], ..., x[i]]
        STOP;

  • 因为第 n 个插值多项式是唯一的,所以 \(N_n(x) \equiv P_n(x)\)
  • 它们必须有相同的截断误差,即:

    \[ \begin{align} & f[x, x_0, \dots, x_n] \omega_{k+1} (x) = \dfrac{f^{(n+1)}(\xi_x)}{(n+1)!} \omega_{k+1}(x) \notag \\ & \Rightarrow f[x_0, \dots, x_k] = \dfrac{f^{(k)}(\xi)}{k!}, \xi \in (x_{\text{min}}, x_{\text{max}}) \notag \end{align} \]
  • 该过程和 Neville 法类似:

    \[ \begin{matrix} f(x_0) & & & & & \notag \\ f(x_1) & f[x_0, x_1] & & & & \notag \\ f(x_2) & f[x_1, x_2] & f[x_0, x_1, x_2] & & & \notag \\ \dots & \dots & \dots & & & \notag \\ f(x_{n-1}) & \dots & \dots & & & \notag \\ f(x_n) & f[x_{n-1}, x_n] & f[x_{n-2}, x_{n-1}, x_n] & & f[x_0, \dots, x_n] \notag \\ f(x_{n+1}) & f[x_n, x_{n+1}] & f[x_{n-1}, x_n x_{n+1}] & \dots & f[x_1, \dots, x_{n+1}] & f[x_0, \dots, x_{n+1}] \notag \end{matrix} \]

Formulae with Equal Spacing⚓︎

如果这些点是等间距的,即 \(x_i = x_0 + ih\ (i = 0, \dots, n)\),那么:

  • 前向差(forward difference):\(\Delta f_i = f_{i+1} - f_i, \Delta^k f_i = \Delta(\Delta^{k-1} f_i) = \Delta^{k-1} f_{i+1} - \Delta^{k-1} f_i\)
  • 后向差(backward difference)::\(\nabla f_i = f_i - f_{i-1}, \nabla^k f_i = \nabla(\Delta^{k-1} f_i) = \nabla^{k-1} f_i - \nabla^{k-1} f_{i-1}\)
  • 中心差(centered difference):\(\delta^k f_i = \delta^{k-1} f_{i+\frac{1}{2}} - \delta^{k-1} f_{i - \frac{1}{2}}\),其中 \(f_{i \pm \frac{1}{2}} = f(x_i \pm \dfrac{h}{2})\)

Some Important Properties⚓︎

  • 线性:\(\Delta(a \cdot f(x) + b \cdot g(x)) = a \Delta f + b \Delta g\)
  • 如果 \(f(x)\) 是一个 \(m\) 阶多项式,那么 \(\Delta^k f(x)\ (0 \le k \le m)\) 是一个 \(m - k\) 阶多项式且 \(\Delta^k f(x) = 0\ (k > m)\)
  • 差值还能从以下函数中得到:

    • \(\Delta^n f_k = \sum\limits_{j=0}^n (-1)^j \left( \begin{array}{cccc}n \\ j\end{array}\right) f_{n+k-j}\)
    • \(\nabla^n f_k = \sum\limits_{j=0}^n (-1)^{n-j} \left( \begin{array}{cccc}n \\ j\end{array}\right) f_{k+j-n}\)
  • 反之亦然:\(f_{n+k} = \sum\limits_{j=0}^n \left( \begin{array}{cccc}n \\ j\end{array}\right) \Delta^j f_k\)

  • \(f[x_0, \dots, x_k] = \dfrac{\Delta^k f_0}{k! h^k}, f[x_n, x_{n-1}, \dots, x_{n-k}] = \dfrac{\nabla^k f_n}{k!h^k}\)。从 \(R_n\) 可以得到:\(f^{(k)}(\xi) = \dfrac{\Delta^k f_0}{h^k}\)

总结

  • 牛顿前向差公式(Newton forward-difference formula):令 \(x = x_0 + th\),那么

    \[ \begin{align} N_n(x) & = N_n(x_0 + th) = \sum\limits_{k=0}^n \left( \begin{array}{cccc}t \\ k\end{array}\right) \Delta^k f(x_0), \notag \\ R_n(x) & = \dfrac{f^{(n+1)} (\xi)}{(n+1)!} t(t-1) \dots (t-n)h^{n+1}, \xi \in (x_0, x_n) \notag \end{align} \]
  • 牛顿后向差公式(Newton backward-difference formula):颠倒点的顺序,即计算 \(N_n(x) = f(x_n) + f[x_n, x_{n-1}](x - x_n) + dots + f[x_n, \dots, x_0](x - x_n) \dots (x - x_1)\)。令 \(x = x_n + th\),那么

    \[ N_n(x) = N_n(x_n + th) = \sum\limits_{k=0}^n (-1)^k \left( \begin{array}{cccc}-t \\ k\end{array}\right) \nabla^k f(x_n) \]

Hermite Interpolation⚓︎

密切多项式

\(x_0, x_1, \dots, x_n\) 为在 \([a, b]\) 上的 \(n+1\) 个不同的数,\(m_i\) 是和 \(x_i\) 关联的非负整数(\(i = 0, \dots, n\)。假设 \(f \in C^m[a, b]\),其中 \(m = \max\limits_{0 \le i \le n}m_i\),那么用于近似 \(f\) 密切多项式(osculating polynomial) 为满足以下条件的阶数最小的多项式 \(P(x)\)

\[ \dfrac{d^k P(x_i)}{dx^k} = \dfrac{d^k f(x_i)}{dx^k},\ \text{for each } i = 0, \dots, n \text{ and } k = 0, \dots, m \]

目标:找到一个密切多项式 \(P(x)\),使得 \(\forall i = 0, 1, \dots, n, P(x_i) = f(x_i), P'(x_i) = f'(x_i), \dots, P^{(m_i)}(x_i) = f^{(m_i)}(x_i)\)

  • 给定 \(N\) 个条件(即有 \(N\) 个方程\(N - 1\) 阶多项式就能确定下来
  • \(f\) 以及所有在一个点 \(x_0\) 上的 \(\le m_0\) 阶的导数吻合的密切多项式就是一个泰勒多项式:

    \[ P(x) = f(x_0) + f'(x_0)(x - x_0) + \dots + \dfrac{f^{(m_0)}(x_0)}{m_0!}(x - x_0)^{m_0} \]

    且余项 \(R(x) = f(x) - \varphi(x) = \dfrac{f^{(m_0 + 1)}(\xi)}{(m_0 + 1)!}(x - x_0)^{(m_0 + 1)}\)

  • \(\forall i = 0, 1, \dots, n,\ m_i = 1\) 时,此时的多项式为埃尔米特多项式(Hermite polynomials)

例子

假设 \(x_0 \ne x_1 \ne x_2\)。给定 \(f(x_0), f(x_1), f(x_2)\) \(f'(x_1)\),寻找多项式 \(P(x)\),满足 \(P(x_i) = f(x_i),\ i = 0, 1, 2\),且 \(P'(x_1) = f'(x_1)\)。并分析误差。

首先,\(P(x)\) 的阶必须 \(\le 3\)

与拉格朗日多项式类似,我们要待定一个埃尔米特多项式:\(P_3(x) = \sum\limits_{i=0}^2 f(x_i) h_i(x) + f'(x_1) \hat{h_1}(x)\),其中 \(h_i(x_j) = \delta_{ij}, h_i'(x_1) = 0, \hat{h_1}(x_i) = 0, \hat{h_1}'(x_1) = 1\)

  • \(h_0(x)\):有根 \(x_1, x_2\),且 \(h_0'(x_1) = 0 \quad \Rightarrow \quad x_1\) 是一个重根
    • \(\begin{cases}h_0(x) = C_0(x - x_1)^2(x - x_2) \\ h_0(x_0) = 1 \Rightarrow C_0\end{cases} \quad \Rightarrow \quad h_0(x) = \dfrac{(x - x_1)^2(x - x_2)}{(x_0 - x_1)^2(x_0 - x_2)}\)
  • \(h_2(x)\):与 \(h_0(x)\) 类似
  • \(h_1(x)\):有根 \(x_0, x_2 \Rightarrow h_1(x) = (Ax + B)(x - x_0)(x - x_2)\)\(A, B\) 可通过 \(h_1(x_1) = 0\) \(h_1'(x_1) = 0\) 求解
  • \(\hat{h_1}(x)\)::有根 \(x_0, x_1, x_2 \Rightarrow \hat{h_1}(x) = C_1(x - x_0)(x - x_1)(x - x_2)\)\(h_1(x_1) = 1 \Rightarrow C_1\) 能被求解

一般情况下,给定 \(x_0, \dots, x_n; y_0, \dots, y_n\) 以及 \(y_0', \dots, y_n'\),埃尔米特多项式 \(H_{2n+1}(x)\) 满足对于所有的 \(i\)\(H_{2n+1}(x_i) = y_i\) \(H_{2n+1}'(x_i) = y_i'\)

求解过程

\(H_{2n+1}(x) = \sum\limits_{i=0}^n y_i h_i(x) + \sum\limits_{i=0}^n y_i' \hat{h_i}(x)\),其中 \(h_i(x_j) = \delta_{ij}, h_i'(x_j) = 0, \hat{h_i}(x_j) = 0, \hat{h_i}'(x_j) = \delta_{ij}\)

  • \(h_i(x)\)
    • \(x_0, \dots, \hat{x_i}, \dots, x_n\) 是重数为 2 的根 \(\Rightarrow\ h_i(x) = (A_i x + B_i) L_{n, i}^2(x)\)
    • \(A_i, B_i\) 能通过 \(h_i(x_i) = 1, h_i'(x_i) = 0\) 求解
    • \(h_i(x) = [1 - 2L_{n, i}'(x_i)(x - x_i)L_{n, i}^2(x)]\)
  • \(\hat{x_i}(x)\)
    • 除了 \(x_i\) 外,所有的根 \(x_0, \dots, x_n\) 的重数均为 2,得到:
    • \(\begin{cases}\hat{h_i}(x) = C_i(x - x_i) L_{n, i}^2(x) \\ \hat{h_i}'(x_i) = 1 \Rightarrow C_i = 1\end{cases} \quad \Rightarrow \quad \hat{h_i}(x) = (x - x_i) L_{n, i}^2(x)\)

如果 \(a = x_0 < x_1 < \dots < x_n = b, f \in C^{2n}[a, b]\),那么 \(R_n(x) = \dfrac{f^{(2n+2)}(\xi_x)}{(2n+2)!}\Big[\prod\limits_{i=0}^n (x - x_i) \Big]^2\)

定理

如果 \(f \in C^1[a, b]\) \(x_0, \dots, x_n \in [a, b]\) 是不同的数,那么在函数 \(f\) 及其导数 \(f'\) 一致的最小次数唯一多项式即为次数不超过 \(2n+1\) 的埃尔米特多项式:

\[ H_{2n+1}(x) = \sum\limits_{j=0}^n f(x_j) H_{n,j}(x) + \sum\limits_{j=0}^n f'(x_j) \hat{H}_{n,j}(x) \]

其中 \(H_{n,j} = [1 - 2(x - x_j) L'_{n,j}(x_j)]L_{n,j}^2(x)\)\(\hat{H}_{n,j}(x) = (x - x_j) L^2_{n,j}(x)\)。在这里,\(L_{n,j}(x)\) 指代的是 \(n\) 阶多项式中第 \(j\) 个拉格朗日系数。

另外,若 \(f \in C^{2n-2}[a, b]\),那么

\[ f(x) = H_{2n+1}(x) + \dfrac{(x - x_0)^2 \dots (x - x_n)^2}{(2n+2)!}f^{(2n+2)}(\xi) \]

\(\xi\) 为某个满足 \(a < \xi < b\) 的数。

思考

给定 \(x_i = i + 1, i = 0, 1, 2, 3, 4, 5\),哪一个是 \(\hat{h_2}(x)\)

算法:埃尔米特插值

求得在 \(n+1\) 个不同的数 \(x_0, \dots, x_n\) 上的函数 \(f\) 的埃尔米特插值多项式 \(H(x)\) 的系数。

  • 输入:数 \(x_0, x_1, \dots, x_n\);值 \(f(x_0), f(x_1), \dots, f(x_n)\) 以及 \(f'(x_0), \dots, f'(x_n)\)
  • 输出:数 \(Q_{0,0}, Q_{1, 1}, \dots, Q_{2n+1, 2n+1}\),其中

    \[ \begin{align} H(x) = & Q_{0,0} + Q_{1,1}(x - x_0) + Q_{2,2}(x - x_0)^2 + Q_{3,3}(x - x_0)^2(x - x_1) & + Q_{4,4}(x - x_0)^2(x - x_1)^2 + \dots & + Q_{2n+1, 2n+1}(x - x_0)^2(x - x_1)^2 \dots (x - x_{n-1})^2(x - x_n) \end{align} \]
Step 1  for i = 1, 2, ..., n do Step 2 and 3:
    Step 2  Set z[2*i] = x[i];
                z[2*i+1] = x[i];
                Q[2*i][0] = f(x[i]);
                Q[2*i+1][0] = f(x[i]);
                Q[2*i+1][1] = f`(x[i]);
    Step 3  if i != 0 then set Q[2*i][1] = (Q[2*i][0] - Q[2*i-1][0]) / (z[2*i] - z[2*i-1]);
Step 4  for i = 2, 3, ..., 2*n+1:
            for j = 2, 3, ..., i set Q[i][j] = (Q[i][j-1] - Q[i-1][j-1]) / (z[i] - z[i-j]);
Step 5  Output(Q[0][0], Q[1][1], ..., Q[2*n+1][2*n+1]);
        STOP.

Cubic Spline Interpolation⚓︎

例子

考虑关于函数 \(f(x) = \dfrac{1}{1 + x^2}\) 在点 \(x_i = -5 + \dfrac{10}{n}i \in [-5, 5] \ (i = 0, \dots, n)\) 的拉格朗日多项式 \(P_n(x)\)

可以看到,我们无法用多项式(这些彩色曲线)较为准确地近似函数(黑色曲线。正如前面所说,增加多项式的阶数不能保证更好的近似结果,因为高阶多项式具有密切的性质。

一些尝试

思路:在每个子区间 \([x_i, x_{i+1}]\) 上,通过线性多项式近似表示 \(f(x)\),即:

\[ f(x) \approx P_1(x) = \dfrac{x - x_{i+1}}{x_i - x_{i+1}}y_i + \dfrac{x - x_i}{x_{i+1} - x_i} y_{i+1} \text{ for } x \in [x_i, x_{i+1}] \]

\(h = \max |x_{i+1} - x_i|\),那么 \(P_1^h(x) \xrightarrow{\text{uniform}}, h \rightarrow 0\)

缺点:不够平滑

思路:给定 \(x_0, \dots, x_n;\ y_0, \dots, y_n;\ y_0', \dots, y_n'\),在区间 \([x_i, x_{i+1}]\) 的两个端点上构造一个关于 \(y, y'\) 3 阶埃尔米特多项式

缺点:计算导数不太容易

这里介绍一种更好的方法:三次样条插值(cubic spline interpolation)。

定义

给定一个定义在 \([a, b]\) 上的函数 \(f\),以及一组节点 \(a = x_0 < x_1 \dots < x_n = b\),关于 \(f\) 三次样条插值器(cubic spline interpolant) \(S\) 是一个满足下面条件的函数:

  • \(S(x)\) 是一个三次多项式,记作 \(S_i(x)\),在子区间 \([x_i, x_{i+1}]\) 上(\(i = 0, 1, \dots, n - 1\)
  • \(S(x_i) = f(x_i),\ i = 0, 1, \dots, n\)
  • \(S_{i+1}(x_{i+1}) = S_i(x_{i+1}),\ i = 0, 1, \dots, n - 2\)
  • \(S_{i+1}'(x_{i+1}) = S_i'(x_{i+1}),\ i = 0, 1, \dots, n - 2\)
  • \(S_{i+1}''(x_{i+1}) = S_i''(x_{i+1}),\ i = 0, 1, \dots, n - 2\)

Method of Bending Moment⚓︎

\(h_j = x_j - x_{j-1}\) 且 对于 \(x \in [x_{j-1}, x_j],\ S(x) = S_j(x)\),那么 \(S_j''(x)\) 是一个 1 阶多项式,并能通过 \(f\) 上的 2 个节点值确定下来。

  • \(S_j(x)\) 是一个三次多项式

假设 \(S_j''(x_{j-1}) = M_{j-1}, S_j''(x_j) = M_j\)弯矩(bending moment),那么 \(\forall x \in [x_{j-1}, x_j]\)\(S_j''(x) = M_{j-1} \dfrac{x_j - x}{h_j} + M_j \dfrac{x - x_{j-1}}{h_j}\)

\(S_j''\) 积分两次,我们得到了:

  • \(S_j'(x) = -M_{j-1} \dfrac{(x_j - x)^2}{2h_j} + M_{j-1} \dfrac{(x - x_{j-1})^2}{2h_j} + A_j\)
  • \(S_j(x) = M_{j-1} \dfrac{(x_j - x)^3}{6h_j} + M_{j-1} \dfrac{(x - x_{j-1})^3}{6h_j} + A_jx + B_j\)

其中 \(A_j, B_j\) 能通过方程 \(S_j(x_{j-1}) = y_{j-1}, S_j(x_j) = y_j\) 求解。可以得到:

  • \(A_j = \dfrac{y_j - y_{j-1}}{h_j} - \dfrac{M_j - M_{j-1}}{6}h_j\)
  • \(A_j x + B_j = (y_{j-1} - \dfrac{M_{j-1}}{6} h_j^2) \dfrac{x_j - x}{h_j} + (y_j - \dfrac{M_j}{6}h_j^2)\dfrac{x - x_{j-1}}{h_j}\)

现在我们来求解 \(M_j\):因为 \(S'\) \(x_j\) 上是连续的,所以:

  • \([x_{j-1}, x_j]\): \(S_j'(x) = -M_{j-1} \dfrac{(x_j - x)^2}{2h_j} + M_j \dfrac{(x - x_{j-1})^2}{2h_j} + f[x_{j-1}, x_j] - \dfrac{M_j - M_{j-1}}{6}h_j\)
  • \([x_j, x_{j+1}]\): \(S_{j+1}'(x) = -M_j \dfrac{(x_{j+1} - x)^2}{2h_{j+1}} + M_{j+1} \dfrac{(x - x_j)^2}{2h_{j+1}} + f[x_j, x_{j+1}] - \dfrac{M_{j+1} - M_j}{6}h_{j+1}\)

根据 \(S_j'(x_j) = S_{j+1}'(x_j)\),我们可以结合 \(M_{j-1}, M_j, M_{j+1}\) 的系数——定义 \(\lambda_j = \dfrac{h_{j+1}}{h_j + h_{j+1}}, \mu_j = 1 - \lambda_j, g_j = \dfrac{6}{h_j + h_{j+1}} (f[x_j, x_{j+1}] - f[x_{j-1}, x_j])\),可以得到:\(\mu_j M_{j-1} + 2M_j + \lambda_j M_{j+1} = g_j\ (1 \le j \le n - 1)\)。也就是说,我们有 \(n+1\) 个未知数,但只有一个方程;另外我们还需要 2 个额外的边界条件。

  • 固定边界(clamped boundary):\(S'(a) = y_0', S'(b) = y_n'\)

    • \([a, x_1]\): \(S_1'(x) = -M_0 \dfrac{(x_1 - x)^2}{2h_1} + M_1 \dfrac{(x - a)^2}{2h_1} + f[x_0, x_1] - \dfrac{M_1 - M_0}{6}h_1\)
    • \([x_{n-1}, b]\) \(S_n'\) 也是类似的:\(\begin{cases}2M_0 + M_1 = \dfrac{6}{h_1} (f[x_0, x_1] - y_0') = g_0 \\ M_{n-1} + 2M_n = \dfrac{6}{h_n} (y_n' - f[x_{n-1}, x_n]) = g_n\end{cases}\)
    算法:固定三次样条

  • 自由边界(free boundary):\(S''(a) = y_0'' = M_0, S''(b) = y_n'' = M_n\),且 \(M_0 = M_n = 0\)

    • 那么 \(\lambda_0 = 0, g_0 = 2y_0'';\ \mu_n = 0, g_n = 2y_n''\)
    • 此时的样条称为自然样条(natural spline)
    算法:自然三次样条

  • 周期边界(periodic boundary):如果 \(f\) 是周期函数,即 \(y_n = y_0\) \(S'(a^+) = S'(b^-) \Rightarrow M_0 = M_n\)

  • 只要系数矩阵是严格对角占优的,那么三次样条能通过边界被唯一确定
  • 如果 \(f \in C[a, b]\) \(\dfrac{\max h_i}{\min h_i} \le C < \infty\),那么当 \(h_i \rightarrow 0\) 时,\(S(x) \xrightarrow{\text{uniform}} f(x)\)。也就是说,在保证不增加样条阶数的情况下,可通过增加节点个数来提升近似精度

算法概述:三次样条插值法

  1. 计算 \(\mu_j, \lambda_j, g_j\)
  2. 求解 \(M_j\)
  3. 找到包含 \(x\) 的子区间,即找到相应的 \(j\)
  4. 通过 \(S_j(x)\) 得到 \(f(x)\) 的近似值

对应的作业练习📝

评论区

如果大家有什么问题或想法,欢迎在下方留言~