我们提出三个有关函数表达式的问题, 来说明本章及第 4 章 (函数逼近) 中论述的主题.

第一个问题, 假定已有一个函数的数值表如表 3-1, 要问: 是否能找到一个简单且便于计算的公式, 利用它可以算出给定区间上任意点的值.

表 3-1 一般插值数据表

第二个问题, 设给定一个函数 $f$ , $f$ 的表达式非常复杂, 计算 $f$ 的值很不经济. 在这种情况下, 就要寻找另一个函数 $p$ , 它既易于求值且又是对 $f$ 的一个合理的逼近.

第三个问题与第一个问题类似, 但假定表中给出的数值带有误差. 比如当这些值来自于物理实验时, 就可能出现这种情况. 现在要寻找一个公式, 使得它可以近似地表示这些数据.

在这三个问题中, 都能得到一个简单的函数 $p$ 来表示或逼近给定的数值表或函数 $f$ . 简单的函数 $p$ 的类型可以是多项式、样条函数、三角函数及有理函数等.

对第一个问题的解答, 就是插值方法. 本章重点讨论函数 p 为多项式和样条函数的插值方法.

对第二个问题的解答, 就是连续函数逼近的方法. 根据度量逼近的标准不同, 逼近又分为一致逼近和平方逼近, 将在第 4 章讨论. 第三个问题实际上是离散函数的逼近问题, 也将在第 4 章讨论.

3.1 多项式插值#

3.1.1 多项式插值问题的定义#

按条件 (3.1) 求函数 $f(x)$ 的近似表达式 $p(x)$ 的方法称为插值法, 称条件 (3.1) 为插值条件, $x_{i}(i=0,1,2,\cdots,n)$ 为插值节点, 包含插值节点的区间 [a,b] 为插值区间.

多项式插值有明确的几何意义, 即通过平面上给定的 $n+1$ 个互异点 $(x_{i},y_{i})(i=0,1,2,\cdots,n)$ , 作一条 n 次多项式曲线 $y=p(x)$ 近似地表示曲线 $y=f(x)$ , 见图 3-1. 多项式插值问题也称为代数插值问题.

图 3-1 多项式插值的几何意义

证：设 n 次多项式 $p(x)=a_{0}+a_{1}x+\cdots+a_{n}x^{n}$ 是函数 $y=f(x)$ 在 $[a,b]$ 上的 $n+1$ 个互异的节点 $x_{i}(i=0,1,2,\cdots,n)$ 上的插值多项式，则求 $p(x)$ 的问题就可归结为求以它的系数 $a_{i}(i=0,1,2,\cdots,n)$ 为未知元的 $n+1$ 阶线性方程组

其系数行列式是范德蒙德 (Vandermonde) 行列式:

因为 $x_{i}$ 互不相同, 所以式 (3.3) 不为零, 根据解线性方程组的克拉默 (Cramer) 法则, 方程组的解 $a_{i}$ 存在且唯一, 从而 $p(x)$ 被唯一确定, 这就证明了 $n$ 次代数插值问题的解是存在且唯一的.

上述证明实质上给出了求代数插值多项式的一个方法. 但此方法不适用于计算机求解: 范德蒙德行列式的计算需要的工作量太大. 唯一性保证不论用什么方法求 $p(x)$ , 其结果总是一样. 本章主要介绍拉格朗日 (Lagrange) 插值和牛顿 (Newton) 插值.

3.1.3 插值基函数#

所要求的插值多项式 $p(x)$ 是线性空间 $P_{n}(x)$ (次数小于等于 $n$ 的代数多项式的全体) 中的一个点, 根据线性代数中线性空间的理论, 全体次数小于等于 $n$ 的代数多项式构成的 $n + 1$ 维线性空间的基底是不唯一的, 因此, $n$ 次代数插值多项式 $p(x)$ 可以写成多种形式. 由线性空间的不同基底出发, 构造满足插值条件的多项式的方法称为基函数法. 基函数法求插值多项式分两个步骤. 首先, 定义 $n+1$ 个线性无关的特殊代数多项式, 它们在插值理论中称为插值基函数, 用 $\varphi_{0}(x), \cdots, \varphi_{n}(x)$ 表示. 其次, 利用插值条件, 确定插值基函数的线性组合表示的 n 次插值多项式

的系数 $a_{0},\cdots,a_{n}$ .

3.2 拉格朗日插值#

3.2.1 拉格朗日插值基函数#

$l_{k}(x)$ 很容易找到. 因为 n 个节点 $x_{i}(i=0,1,\cdots,k-1,k+1,\cdots,n)$ 是 $l_{k}(x)$ 的零点, 因此可设

其中 $A_{k}$ 为待定系数. 由条件 $l_{k}(x_{k})=1$ 可定 $A_{k}$ ，于是

取 $k=0,1,2,\cdots,n,$ 就得到 $n+1$ 个拉格朗日插值基函数.

3.2.2 拉格朗日插值多项式#

以 $l_{0}(x), l_{1}(x), \cdots, l_{n}(x)$ 的线性组合表示的 n 次多项式为

其中， $a_{0},a_{1},\cdots,a_{n}$ 为待定系数，令

利用 $l_{j}(x)$ 的性质得 $a_{i}=y_{i}\ (i=0,1,\cdots,n)$ ，将其代入式 (3.7) 所得的多项式称为拉格朗日插值多项式，记作 $L_{n}(x)$ 。即

当 $n = 1$ 时, 拉格朗日插值多项式 (3.8) 为 $L_{1}(x) = l_{0}(x)y_{0} + l_{1}(x)y_{1}$ , 即

用 $L_{1}(x)$ 近似代替 $f(x)$ 称为线性插值, 公式 (3.9) 称为线性插值多项式或一次插值多项式. 当 $n = 2$ 时, 拉格朗日插值多项式 (3.8) 为 $L_{2}(x) = l_{0}(x)y_{0} + l_{1}(x)y_{1} + l_{2}(x)y_{2}$ , 即

用 $L_{2}(x)$ 近似代替 $f(x)$ 称为二次插值或抛物线插值, 称式 (3.10) 为二次插值多项式.

例3.2.2 已知函数 $y = \ln x$ 的几个函数值如表3-2所示. 试分别用线性插值和抛物线插值求 $\ln 11.75$ 的近似值.

表 3-2 函数 $y = \ln x$ 数据表

解：在插值计算中, 为了减少截断误差, 在选择插值节点时, 应尽量选取与插值点 $x$ 距离较近的一些节点. 本题中 $x = 11.75$ 介于点 11 与 12 之间, 故作线性插值时应取节点 $x_0 = 11$ , $x_1 = 12$ , 用线性插值公式 (3.9) 有

将 x = 11.75 代入, 即得

类似地, 在抛物线插值时取节点 $x_{0}=11, x_{1}=12, x_{2}=13$ , 所得 $\ln 11.75$ 的近似值为

查对数表得, $\ln 11.75 \approx 2.46385$ .

由此可见，抛物线插值的精度较线性插值好.

下面我们介绍拉格朗日插值法的 MatLab 实现:

function yh = lagrange(x,y,xh)
    n = length(x);
    m = length(xh);
    yh = zeros(1,m);
    c1 = ones(n-1,1);
    c2 = ones(1,m);
    for i = 1:n

xp = x([1:i-1 i+1:n]);
yh = yh + y(i) * prod((c1*xh-xp'*c2)./(x(i)-xp'*c2));

上面的程序中, prod 命令对输入的一个向量返回其所有分量的乘积. 例题的 Matlab 程序如下.

线性插值:

>> x = [11 12];
>> y = [2.3979 2.4849];
>> xh = [11.75];
>> yh = lagrange(x, y, xh)
抛物线插值：
>> x = [11 12 13];
>> y = [2.3979 2.4849 2.5649];
>> xh = [11.75];
>> yh = lagrange(x, y, xh)

与准确值相比较, 发现线性插值的误差比抛物线插值的误差大. 这里提出两个问题. (1) 通常 $f(x)$ 的准确值是不知道的, 怎样估计用 $L_{n}(x)$ 近似代替 $f(x)$ 时所产生的误差? (2) 是不是插值多项式的次数越高, 其计算结果就越精确?

3.2.3 插值余项#

我们首先来解决上述的第一个问题, 即怎样估计用 $L_{n}(x)$ 近似代替 $f(x)$ 时所产生的误差? 在区间 $[a,b]$ 上用插值多项式 $L_{n}(x)$ 近似 $f(x)$ 时, 由插值条件, $f(x)$ 和 $L_{n}(x)$ 在插值节点 $x_{0}, x_{1}, \cdots, x_{n}$ 上的差应为 0, 设 $R_{n}(x) = f(x) - L_{n}(x)$ , 则在节点 $x_{i}$ 上有

即 $R_{n}(x)$ 在 $[a,b]$ 上至少有 $n + 1$ 个零点, 故设

其中, $K(x)$ 为待定系数, $\Pi(x)=(x-x_{0})(x-x_{1})\cdots(x-x_{n})$ .

引进辅助函数 $\varphi(t)=f(t)-L_{n}(t)-K(x)\Pi(t)$ ，视 x 为 $(a,b)$ 上的一个固定点，且 $x\neq x_{i}(i=0,1,\cdots,n)$ 。则 $\varphi(t)$ 在 $n+2$ 个点 $x,x_{0},x_{1},\cdots,x_{n}$ 上取值为 0。根据罗尔 (Rolle) 定理，在 $\varphi(t)$ 的任何两个零点之间必存在一点 $\eta$ ，使 $\varphi'(\eta)=0$ ，于是 $\varphi'(x)$ 在 $(a,b)$ 上至少有 $n+1$ 个零点。对 $\varphi'(t),\varphi''(t),\cdots,\varphi^{(n)}(t)$ 反复运用罗尔定理，最后推出，在 $(a,b)$ 上至少存在一个点 $\xi$ ，使 $\varphi^{(n+1)}(\xi)=0$ ，而

其中 $L_{n}(x)$ 是次数不超过 n 的多项式, $\Pi(x)$ 是 $n+1$ 次多项式, 所以 $L_{n}^{(n+1)}(\xi)=0$ , $\Pi^{(n+1)}(\xi)=(n+1)!$ , 于是

代入式 (3.11) 得

其中 $\xi \in (a,b)$ 且依赖于 $x$ .

称 $R_{n}(x)$ 为 n 次插值多项式的余项.

综上所述, 得到如下定理.

$R_{n}(x)$ 的表达式 (3.12) 中的 $f^{(n + 1)}(\xi)$ 难于确定, 但 $\left|f^{(n + 1)}(x)\right|$ 的上界往往可以估计, 记 $M_{n + 1} = \max_{a\leqslant x\leqslant b}\left|f^{(n + 1)}(x)\right|$ , 则

这是一个实用的插值余项估计式. 如要求估计 $|R_n(x)|$ 在 $[a, b]$ 上的最大值, 还需计算 $\max_{a \leqslant x \leqslant b} |\Pi(x)|$ , 并有

余项表达式 (3.13) 表示 $|R_n(x)|$ 的大小与 $M_{n+1}$ 及 $|\Pi(x)| = |\prod_{i=0}^n (x - x_i)|$ 都有关. 对于指定的插值点 $x$ , 当节点的个数 $m$ 大于 $n + 1(n$ 为插值多项式的次数)时, 应当选取靠近 $x$ 的节点做插值多项式, 以使 $|\Pi(x)|$ 中的诸 $|x - x_i|$ 都较小, 从而余项 $|R_n(x)|$ 也较小.

例 3.2.4 估计上例中用线性插值和抛物线插值计算 ln11.75 的误差.

解：由于 $R_{1}(x) = \frac{f^{\prime\prime}(\xi)}{2!}(x - x_{0})(x - x_{1}),\xi \in [x_{0},x_{1}],f(x) = \ln (x),f^{\prime\prime}(x) = -\frac{1}{x^{2}},\xi$ 介于11与12之间，所以 $|f''(\xi)| = \frac{1}{\xi^2} < \frac{1}{11^2}$ .于是，

由于 $R_{2}(x) = \frac{f^{(3)}(\xi)}{3!} (x - x_{0})(x - x_{1})(x - x_{2}),\xi \in [x_{0},x_{2}],f^{(3)}(x) = \frac{2}{x^{3}},|f^{(3)}(\xi)| = \frac{2}{\xi^{3}} <$ $\frac{2}{11^3}$ ，于是

可见, $L_{2}(11.75)$ 比 $L_{1}(11.75)$ 的误差小.

那么, 是否插值多项式的次数越高, 插值多项式近似原来函数的精度越高? 这个问题我们将在随后的内容中回答.

3.3 牛顿插值#

3.3.1 差商#

若将插值基函数取为

用它们组合成如下形式的 $n$ 次多项式

例如

等等, 其中 $a_{0}, a_{1}, \cdots, a_{n}$ 为待定系数. 为使式 (3.15) 表示的 $N_{n}(x)$ 成为 $f(x)$ 的插值多项式, 需要按插值条件

确定参数 $a_{i}(i=0,1,\cdots,n)$ . 令

可以求得，

等等. 为得到参数 $a_{i}$ 的一般表达式, 我们来给出差商的定义.

差商有以下三条性质.

(1) 差商可以表示为函数值的线性组合. 用数学归纳法可以证明

(2) 差商关于所含节点是对称的. 例如: 在 $k$ 阶差商 $f[x_0, x_1, \cdots, x_{k-1}, x_k]$ 中, 任意调换节点 $x_i$ 与 $x_j$ 的次序, 其值不变, 即

(3) 设 $f(x)$ 在含有 $x_{0}, x_{1}, \cdots, x_{n}$ 的区间 $[a, b]$ 上具有 n 阶导数，则在这一区间上至少有一点 $\xi$ ，使

这一性质揭示了差商与导数之间的关系, 我们将在后面给出它的证明.

差商可以列表计算, 以 $n = 4$ 为例, 差商表如表3-3所示. 表3-3的计算可以按照自上而下, 从低阶到高阶的顺序完成.

表 3-3 差商表

用 Matlab 计算差商时, 用二维数组 $P(n, n + 1)$ 保存上述差商表. 其中, $n$ 是函数表中节点的个数. 数组 $P$ 的第 1 列为节点值 $x$ , 第二列为函数值 $y$ , 第 $j$ 列 $(3 \leqslant j \leqslant n + 1)$ 为 $j - 2$ 阶差商. 一维数组 $q(n)$ 为上述差商表中带下划线的差商 (规定 $f(x_i) = y_i$ 为 $x_i$ 处的零阶差商). Matlab 计算差商的自定义函数如下.

function [p,q]=chashang(x,y)
    n = length(x);
    p(:,1) = x;
    p(:,2) = y;
    for j = 3:n+1
    p(1:n+2-j,j) = diff(p(1:n+3-j,j-1)) ./(x(j-1:n)-x(1:n+2-j));
    end
    q = p(1,2:n+1);

Matlab 中的命令 diff 给出一个向量所有相邻分量的差 (后一个分量减去前一个分量), 因此计算结果是一个比输入向量长度少 1 的向量, 例如

>> diff([1 3 6 10 15])
ans =

3.3.2 牛顿插值公式及其余项#

有了差商的定义后, 可以证明式 (3.15) 中的系数 $a_{0}, a_{1}, \cdots, a_{n}$ 分别为

从而获得数据 $\{(x_{i},f(x_{i})\}_{i = 0}^{n}$ 上的 $n$ 次牛顿插值多项式

现在转向解决牛顿插值余项问题, 即估计用 $N_{n}(x)$ 近似代替 $f(x)$ 时所产生的误差.

设 $x \in [a, b]$ , $x \neq x_{i} (i = 0, 1, \cdots, n)$ , 视 x 为一个节点. 设 $N_{n}(t)$ 是数据 $\{(x_{i}, f(x_{i}))\}_{i=0}^{n}$ 上的 n 次牛顿插值多项式, $N_{n+1}(t)$ 是数据 $\{(x_{0}, f(x_{0})), (x_{1}, f(x_{1})), \cdots, (x_{n}, f(x_{n})), (x, f(x))\}$ 上的 $n+1$ 次牛顿插值多项式. 则由式 (3.16), $n+1$ 次多项式 $N_{n+1}(t)$ 与 n 次多项式 $N_{n}(t)$ 有下列关系

由于 $N_{n+1}(t)$ 是数据 $\{(x_{0},f(x_{0})),(x_{1},f(x_{1})),\cdots,(x_{n},f(x_{n})),(x,f(x))\}$ 上的 $n+1$ 次牛顿插值多项式, 则在点 x 处, 一定满足插值条件

在式 (3.17) 的两边, 自变量 $t$ 用 $x$ 来代替, 得

即

从而牛顿插值多项式的余项

根据插值多项式的唯一性, 对于同一组数据 $\{(x_{i}, f(x_{i})\}_{i=0}^{n}$ 上的 n 次插值多项式 $L_{n}(x)$ 和 $N_{n}(x)$ , 应有 $N_{n}(x) = L_{n}(x)$ . 因此, 若 $f(x)$ 在 $[a, b]$ 上的 $n + 1$ 阶导数存在, 则余项

且 $\xi$ 依赖于 $x$ . 同时得到差商与导数有关系

因这里视 x 为一个节点, 故可以推广到一般情形, 对于 $f[x_{0}, x_{1}, \cdots, x_{k}], k = 1, 2, \cdots, n,$ 必存在位于这些节点之间的 $\xi$ , 使得

式 (3.18) 即为差商的性质 3.

例 3.3.2 用牛顿线性插值和抛物线插值计算 $\ln 11.75$ ，已知函数 $y = \ln x$ 的函数表如前面的表 3-2 所示.

解: 由于 11.75 位于 11 与 12 之间, 可以取节点 $x_0 = 11$ , $x_1 = 12$ 作线性插值, 取节点 $x_0 = 11$ , $x_1 = 12$ , $x_2 = 13$ 作二次插值, 由给定数据作差商表 (表 3-4). 因此可得

线性插值计算 $\ln 11.75\approx N_1(11.75)\approx 2.4632.$

二次插值计算 $\ln 11.75 \approx N_2(11.75) \approx 2.4638$ .

表 3-4 函数 $y = \ln x$ 的差商表

从牛顿插值公式 (3.16) 看出, 当增加一个节点时, 公式只在原来基础上增加一项, 前面计算的结果依然可以用, 因而与拉格朗日插值相比较, 牛顿插值公式具有灵活增加节点的优点.

3.3.3 差分与等距节点的插值公式#

牛顿插值公式可以应用于非等距节点的情形, 但实际应用时, 常采用等距节点. 这时, 利用差分可以进一步简化牛顿插值公式, 导出计算上更为有效的等距节点的插值公式.

设已知函数 $f(x)$ 在等距节点 $x_{i}=x_{0}+ih\ (i=0,1,2,\cdots,n)$ 上的函数值为 $f(x_{i})=f_{i}$ ，式中 h>0 称为步长.

定义3.3.3 称 $\Delta f_{i} = f(x_{i} + h) - f(x_{i}) = f_{i + 1} - f_{i}$ 为函数 $f(x)$ 在点 $x_{i}$ 处步长为 $h$ 的一阶向前差分.

称 $\Delta^2 f_i = \Delta (\Delta f_i) = \Delta f_{i + 1} - \Delta f_i$ 为函数 $f(x)$ 在点 $x_{i}$ 处的二阶向前差分.

一般地, 设 $n - 1$ 阶差分已定义, 则称

为函数 $f(x)$ 在点 $x_{i}$ 处的 n 阶向前差分, 并规定 $\Delta^{0}f_{i}=f_{i}$ 为 $f(x)$ 在点 $x_{i}$ 处的零阶差分.

类似地, 还可以定义向后差分和中心差分.

由差分的定义, 也可列出类似于差商表的差分表, 并可导出差分与差商的关系. 例如, k 阶差商与 k 阶差分之间有关系

利用差分与差商的关系式 (3.19), 在牛顿插值公式 (3.16) 中, 将差商替换为差分, 并令 $x = x_0 + th (0 \leqslant t \leqslant 1)$ , 可得

这个用各阶向前差分表示的插值公式 (3.20) 称为牛顿向前插值公式.

3.4 埃尔米特插值#

前面介绍的代数插值, 只要求插值多项式 $p(x)$ 满足条件 (3.1), 但是, 这种插值多项式往往还不能全面反映被插值函数 $f(x)$ 的性态, 许多实际问题不但要求插值函数在某些节点或全部节点上与 $f(x)$ 的导数值也相等, 甚至要求高阶导数值也相等. 这种插值问题为埃尔米特插值问题. 本节不去讨论一般的埃尔米特插值问题, 而只讨论下面这种最简单也是最常用的情形.

3.4.1 两点三次埃尔米特插值#

这个问题的几何意义是十分明显的, 即不仅要求代数曲线 $y = H(x)$ 与函数曲线 $y = f(x)$ 在 $n + 1$ 个互异的点 $(x_{i}, y_{i})$ 处完全重合, 而且还要求它们有公切线.

已知函数 $f(x)$ 在节点 $x_0, x_1$ 上的函数值以及一阶导数值如表 3-5. 求一个三次埃尔米特插值多项式 $H_3(x)$ ，使其满足

表 3-5 两点埃尔米特插值数据

设 $\alpha_0(x),\alpha_1(x),\beta_0(x),\beta_1(x)$ 都是三次多项式, 并且满足

它们的线性组合记为 $H_{3}(x)$ ，即

其中 $a_0, a_1, b_0, b_1$ 为待定参数. 令 $H_3(x)$ 满足插值条件 (3.21), 根据 $\alpha_j(x)$ 和 $\beta_j(x)$ 的性质 (3.22)~(3.25), 立即得到 $a_j = y_j$ ( $j = 0, 1$ ), $b_j = y_j'$ ( $j = 0, 1$ ). 于是

由 $\alpha_{j}(x)$ 的性质 (3.22) 和 (3.23) 可以看出, 节点 $x_{1}, x_{0}$ 依次是 $\alpha_{0}(x)$ 和 $\alpha_{1}(x)$ 的二重零点, 并且 $\alpha_{0}(x_{0}) = \alpha_{1}(x_{1}) = 1$ . 因此, $\alpha_{j}(x)$ 可以用一次拉格朗日插值基函数表示. 又因 $\alpha_{j}(x)$ 是三次多项式, 故可设

其中, $a$ 和 $b$ 为待定参数. 由 $\alpha_{j}(x_{j}) = 1$ 及 $\alpha_{j}^{\prime}(x_{j}) = 0$ , 解出

代入式 (3.27) 经整理得

类似地, 设

利用条件 (3.24) 和 (3.25) 得

代入式 (3.30) 得

将式 (3.29) 和 (3.31) 代入式 (3.26), 得两个节点的三次埃尔米特插值多项式

3.4.2 埃尔米特插值多项式的余项#

设 $R_{3}(x) = f(x) - H_{3}(x)$ ，由插值条件(3.21)知

因此, 节点 $x_{0}, x_{1}$ 都是 $R_{3}(x)$ 的二重零点, 故设

其中 $K(x)$ 为待定函数. 类似于拉格朗日插值余项的推导, 引进辅助函数

则 $\varphi(t)$ 至少有 5 个零点 (二重零点算 2 个), 利用罗尔定理, 至少存在一点 $\xi \in [x_{0}, x_{1}]$ , 使

于是 $K(x) = \frac{f^{(4)}(\xi)}{4!}$ , 代入式 (3.33) 得

综上所述得到如下定理.

定理3.4.2 设 $f^{(4)}(x)$ 在 $[x_0, x_1]$ 上存在且连续，则对于任意的 $x \in [x_0, x_1]$ ，都存在 $\xi \in [x_0, x_1]$ ，使式(3.34)成立.

3.4.3 $n + 1$ 个点 $2n + 1$ 次埃尔米特插值多项式 $H_{2n + 1}(x)$ 及其余项 $R_{2n + 1}(x)$#

类似两点三次埃尔米特插值多项式 $H_{3}(x)$ 和余项 $R_{3}(x)$ 的推导, 我们可以得到 $n + 1$ 个点 $2n + 1$ 次埃尔米特插值多项式 $H_{2n + 1}(x)$ 的表达式:

其中 $l_{k}(x)$ 为节点 $\{x_{i}\}_{i=0}^{n}$ 上的拉格朗日基函数.

$2n + 1$ 次埃尔米特插值多项式 $H_{2n + 1}(x)$ 的余项 $R_{2n + 1}(x)$ 的表达式为

其中, $\Pi(x) = \prod_{j=0}^{n}(x - x_j)$ , $\xi \in (a,b)$ , 且 $\xi$ 与 $x$ 有关. (这里假定 $f(x)$ 在插值区间 $[a,b]$ 上有 $2n+2$ 阶导数.)

3.5 三次样条插值#

3.5.1 样条插值概念的产生#

1. 插值多项式的次数问题#

在例 3.2.2 中, 线性插值的误差比抛物线插值的误差大, 由此人们可能会产生这样的认识: 插值多项式的次数越高, 精度越好. 但是实际情况并非如此.

20 世纪初, 龙格 (Runge) 对函数 $f(x)=\frac{1}{1+x^{2}}(-5\leqslant x\leqslant5)$ , 取等距的插值节点 $x_{k}=-5+kh\left(h=\frac{10}{n}, k=0,1,\cdots,n\right)$ 做拉格朗日插值多项式

当 n = 10, 观察 $L_{10}(x)$ 与 $f(x)$ 之间的差别可以用如下程序来实现.

$$

x = [ - 5: 1: 5 ]; $$

$$

\mathrm{y} = 1. / (1 + \mathrm{x}. ^{\wedge} 2); $$

$$

\mathrm{xh} = [ - 5: 0.25: 5 ]; $$

$$

\mathrm{yh} = \text{lagrange} (\mathrm{x}, \mathrm{y}, \mathrm{xh}); $$

$$

\mathrm{x} 1 = [ - 5: 0.25: 5 ]; $$

$$

\mathrm{y1} = \text{interp1} (\mathrm{x}, \mathrm{y}, \mathrm{x1}); $$

$$

\mathrm{x} = \mathrm{x} 1; $$

$$

\mathrm{y} = 1. / (1 + \mathrm{x}. ^{\wedge} 2); $$

$$

\text{plot} (x h, y h, ^{\prime} r -- ^{\prime}, x 1, y 1, ^{\prime} g o^{\prime}, x, y) $$

上述程序中, lagrange 函数是在 §3.2 中用 Matlab 生成的自定义函数. 运行该程序后, 所得图形见图 3-2.

图 3-2 龙格函数及其 10 阶拉格朗日插值多项式与分段线性多项式

从图 3-2 看出：在接近区间两端点附近， $L_{10}(x)$ 与 $f(x)$ 的偏离很大。比如 $L_{10}(\pm4.8)=1.804\ 38$ ，而 $f(\pm4.8)=0.416\ 0$ 。

龙格还进一步证明了, 在节点等距的条件下, 当 $n \to \infty$ 时, 由式 (3.8) 表示的插值多项式 $L_{n}(x)$ 只在 $|x| \leqslant 3.63$ 内收敛, 在此范围以外, 有

一般地, 把插值多项式不收敛的现象称作龙格现象. 龙格现象说明, 并非插值多项式的次数越高, 其精度就越高.

在处理实际问题时, 人们总希望将所得到的数据点用得越多越好. 解决这一问题的办法是改用直线将节点处的函数值连接起来. 图 3-2 中由空心圆点连接起来的折线近似 $f(x)$ 的程度比 $L_{10}(x)$ 近似 $f(x)$ 的效果要好. 这个解决问题的思想即是分段低次插值.

2. 分段低次插值#

基本思想：用分段多项式来代替单个多项式作插值.

具体作法：插值时，先把整个区间分成若干个小区间；然后在每个小区间上分别作低次插值多项式，例如可用线性插值、抛物线插值、三次插值；再将每个小区间上的插值函数拼接在一起作为整个插值区间上的插值函数.

虽然, 分段低次插值具有公式简单、运算量节省、稳定性好、收敛性有保证等优点, 只要每一个小区间的长度取得足够小, 分段低次插值总能满足所要求的精度. 但是分段低次插值的明显的缺点是节点处的导数值不连续, 即插值多项式不光滑. 在实际问题中, 船体、飞机的外型曲线出于美观和减少阻力的要求必须保持光滑, 即提出这样一个问题: 如何使用分段低次插值多项式, 同时使插值函数具有一定的光滑性. 解决这个问题的办法是使用样条插值函数.

3. 样条插值函数#

样条 (spline) 在英语中指的是富有弹性的细长木条. 样条曲线是指工程师制图时, 用压铁将样条固定在样点上, 其他地方让它自由弯曲, 然后画下的长条曲线. 样条函数的数学实质是由一些按照某种光滑性条件分段拼接起来的多项式组成的函数. 最常用的样条函数是三次样条. 在这个情形, 我们将三次多项式拼接在一起使得所得的样条函数处处有二阶连续导数.

3.5.2 三次样条函数#

下面转向用三次样条插值一个给定的函数值表, 即给定函数 $f(x)$ 在 $[a,b]$ 上的节点 $a = x_{0} < x_{1} < \cdots < x_{n} = b$ 及节点上的函数值

求 $f(x)$ 的三次样条插值函数 $s(x)$ ，使其满足

在给出三次样条插值函数 $s(x)$ 的具体求法以前, 先讨论待定参数的个数与已知的条件个数是否相等, 即该问题的适定性问题.

根据三次样条插值函数的定义, $s(x)$ 是 $[a, b]$ 上的分段三次插值多项式, 即

其中 $s_k(x)$ 应是子区间 $[x_{k-1}, x_k]$ 上的两点三次插值多项式, 故在每个子区间上待定参数的个数为 4, 因共有 $n$ 个子区间, 所以待定参数的个数为 $4n$ .

由于 $s(x)$ 是插值函数, 故满足插值条件

条件 (3.38) 共有 $n + 1$ 个, 又因为 $s(x)$ 是三次样条函数, $s(x) \in C^2[a, b]$ , 故有

其中, $p$ 表示导数的阶数. 这就是说, 在 $(n - 1)$ 个内节点上有条件 (3.39) 成立, 即有 $3(n - 1)$ 个条件. 综合上述两条, 已知的条件个数为 $4n - 2$ .

已知的条件个数 $4n - 2$ 小于待定参数的个数 $4n$ , 该问题在数学上是不适定的. 为使该问题适定, 须补充两个条件. 在实际问题中, 常在边界上补充条件, 故补充的条件常被称为边界条件. 通常使用的边界条件有如下三类.

第一类边界条件是

其中 $f_{0}^{\prime}, f_{n}^{\prime}$ 为给定值.

第二类边界条件是

其中 $f_{0}^{\prime\prime}, f_{n}^{\prime\prime}$ 为给定值. 当 $s^{\prime\prime}(x_{0}) = s^{\prime\prime}(x_{n}) = 0$ 时, 样条函数在两端点不受力, 呈自然状态, 故称之为自然边界条件.

第三类边界条件是周期性条件. 设 $f(x)$ 是周期函数, 不妨设以 $x_{n} - x_{0}$ 为一个周期, 这时 $s(x)$ 也应是以 $x_{n} - x_{0}$ 为周期的周期函数, 于是 $s(x)$ 在端点处满足条件

系统求解三次样条插值函数一般有两种方法：三转角方法和三弯矩方法。三转角方法是从样条函数的一阶导数出发而得到三次样条插值函数的方法，某点的一阶导数在力学上的意义为该点对应截面处的转角，故称该方法为三转角方法。三弯矩方法是从样条函数的二阶导数出发而得到三次样条插值函数的方法，某点的二阶导数在力学上的意义为该点对应截面处的弯矩，故称该方法为三弯矩方法。本章此处给出三弯矩方法的具体过程，三转角方法可参考 §5.7.2 节基于样条的求导方法。

记 $h_j = x_j - x_{j-1}$ , 并设 $s(x)$ 在节点 $x_j$ 处的二阶导数值为 $M_j$ , 即

由于 $s(x)$ 在每个子区间 $[x_{j - 1},x_j]$ 上是一个三次多项式 $s_j(x)$ ，故 $s_j''(x)$ 在区间 $[x_{j - 1},x_j]$ 上是 $x$ 的线性函数，且有 $s_j''(x_{j - 1}) = M_{j - 1}$ ， $s_j''(x_j) = M_j$ ，用线性插值可知其表达式为

将式 (3.43) 积分两次, 得

利用插值条件

确定积分常数 $c_{1}$ 和 $c_{2}$ , 然后代入式 (3.44) 中并整理, 得到 $s(x)$ 在 $[x_{j-1}, x_j]$ 上的表达式为

其中, $j = 1, 2, \cdots, n$ . 由此可知, 只要确定 $M_{0}, M_{1}, \cdots, M_{n}$ 这 $n + 1$ 个待定参数, 所求的三次样条插值函数 $s(x)$ 在各个子区间上的表达式就由式 (3.45) 确定, 从而就求得了在整个 [a, b] 上的三次样条插值函数 $s(x)$ . 并由以上构造过程知, 这时的 $s(x)$ 保证了逐段三次插值在 [a, b] 上的连续性, 并保证了 $s''(x)$ 在内节点的连续性.

为了确定 $M_{0}, M_{1}, \cdots, M_{n}$ ，可利用函数 $s(x)$ 在插值区间 $[a, b]$ 上各内节点 $x_{j} (j = 1, 2, \cdots, n - 1)$ 处有一阶导数连续的条件 $s_{j}^{\prime}(x_{j} - 0) = s_{j+1}^{\prime}(x_{j} + 0)$ ，在区间 $[x_{j-1}, x_{j}]$ 上对 $s_{j}(x)$ 求导，由式 (3.45) 得

从而在右端点 $x_{j}$ 上有

在左端点 $x_{j - 1}$ 上有

将式 (3.47) 中的 $j - 1$ 改为 $j$ , 即得

利用 $s'(x)$ 在内节点连续的条件 $s_{j}'(x_{j}-0)=s_{j+1}'(x_{j}+0)$ ，就得到在 n-1 个内节点 $x_{j}(j=1,2,\cdots,n-1)$ 处关于 $n+1$ 个参数 $M_{j}(j=0,1,2,\cdots,n)$ 的方程

其中

因式 (3.49) 是关于 $n + 1$ 个参数 $M_{j}$ 的 $n - 1$ 个方程, 所以有无穷多组解, 要唯一确定 $n + 1$ 个未知数 $M_{j}$ , 尚须附加边界条件, 补充两个方程.

例如, 对第一种边界条件, 已知插值区间两端的一阶导数值:

利用式 (3.46) 和式 (3.48), 可得到以下两个方程:

将方程 (3.49)、方程 (3.51)、方程 (3.52) 合在一起, 构成了确定 $n+1$ 个参数 $M_{0}, M_{1}, \cdots, M_{n}$ 的 $n+1$ 阶线性方程组:

它的系数矩阵是严格对角占优的, 故有唯一解, 并可用追赶法求解.

对第二种边界条件, $s''(x_{0}) = y''(x_{0})$ , $s''(x_{n}) = y''(x_{n})$ , 即已知 $M_{0} = y_{0}''$ , $M_{n} = y_{n}''$ , 所以, 方程组 (3.49) 只含有 n - 1 个未知元 $M_{1}, M_{2}, \cdots, M_{n-1}$ , 可以表示成

特别当 $M_{0}=M_{n}=0$ 为自然边界条件时, 方程组 (3.54) 的右端成为 $(d_{1},d_{2},\cdots,d_{n-1})^{\mathrm{T}}$ , 其形式特别简单.

对第三种边界条件， $s''(x_0 + 0) = s''(x_n - 0)$ ， $s'(x_0 + 0) = s'(x_n - 0)$ ，可以导出方程 $M_0 = M_n$ 和方程 $\lambda_n M_1 + \mu_n M_{n-1} + 2M_n = d_n$ ，其中

与式 (3.49) 合在一起, 构成以 $M_{1}, M_{2}, \cdots, M_{n}$ 为未知元的线性代数方程组:

方程组 (3.55) 的系数矩阵也是非奇异的, 因而此方程组也有唯一解.

综上所述, 对于给定的函数表 $(x_{i}, f(x_{i}))$ $(i = 0, 1, 2, \cdots, n)$ , 满足第一 (第二或第三) 种边界条件的三次样条插值函数 $s(x)$ 是存在且唯一的, 且其计算过程可归纳如下.

(1) 根据给定的数据 $(x_{i}, f(x_{i}))$ $(i = 0, 1, 2, \cdots, n)$ 及相应的边界条件建立方程组 (3.53) [(3.54) 或 (3.55)].

(2) 解上述线性方程组, 求出 $M_{0}, M_{1}, \cdots, M_{n}$ .

(3) 把求出的 $M_{0}, M_{1}, \cdots, M_{n}$ 代入 $s_{j}(x)$ 的表达式 (3.45)，即得 $s(x)$ 在每一个小区间 $[x_{j-1}, x_{j}]$ 上的分段表达式.

(4) 整个区间 $[a, b]$ 上的三次样条插值函数可表示为

例 3.5.2 设 $f(x)$ 为定义在区间 [0,3] 上的函数, 剖分节点为 $x_{i}=i (i=0,1,2,3)$ , 并给出 $f(x_{0})=0$ , $f(x_{1})=0.5$ , $f(x_{2})=2.0$ , $f(x_{3})=1.5$ , $f'(x_{0})=0.2$ , $f'(x_{3})=-1$ . 试求区间 [0,3] 上满足上述条件的三次样条插值函数 $s(x)$ .

解：利用三弯矩方程组 (3.53) 进行求解, 易知 $h_i = 1, i = 0, 1, 2$ .

于是三弯矩方程组为

求解此方程组, 得 $M_0 = -0.36$ , $M_1 = 2.52$ , $M_2 = -3.72$ , $M_3 = 0.36$ , 代入 (3.45), 经简化得到

用 MatLab 做此题, 程序为

>> x = [0 1 2 3];
>> y = [0.2 0 0.5 2.0 1.5 -1];
>> pp = csape(x,y,'complete')
>> [breaks,coefs,npolys,ncoefs,dim] = unmkpp(pp)
运行结果如下:
pp =
    form: 'pp'
    breaks: [0 1 2 3]
    coefs: [3x4 double]
    pieces: 3
    order: 4
    dim: 1
breaks =
coefs =
    -1.0400   1.2600   1.2800    0.5000
    0.6800   -1.8600   0.6800    2.0000
npolys =
    3
ncoefs =
    4
dim =
    1

其中, coefs 即为 $s(x)$ 的三个分段三次多项式的系数.

上述程序中, 函数 csape(x,y,'complete') 是指定边界条件的样条插值函数, csape 返回一个包含三次样条插值的 pp 形, 或者说是分段多项式形的结构. 这个结构包含了我们希望的任何插值点数值的三次样条值需要的所有信息. 字符串 'complete' 表示所给边界条件是第一边界条件: 若将字符串 'complete' 换作 'second', 表示所给边界条件是第二边界条件; 若将字符串 'complete' 换作 'periodic', 表示所给边界条件是第三边界条件. 在第一边界条件和第二边界条件中, 边界条件值放在 y 的第一个分量和最后一个分量的位置上, 即 y 的分量个数比 x 的分量个数多 2. 比如例题 3.5.2 中, x 与 y 的取法为

>> x = [0 1 2 3];
>> y = [0.2 0 0.5 2.0 1.5 -1];

周期边界条件中, 无须指定边界条件值, x 和 y 即为要插值的函数的节点值及其对应的函数值. 在计算一个三次样条表达式的时候, 必须将 pp 形中的不同域提取出来进行计算, 这个过程可以由函数 unmkpp 很方便地完成. 函数 unmkpp 的使用方法为:

其中输入变量 pp 是样条插值函数 csape 的输出变量, pp 是分段多项式形的结构; unmkpp 的输出变量有 5 个: breaks, coefs, nploys, ncoefs, dim. breaks 包含了插值节点, coefs 是一个矩阵, 其第 i 行是第 i 个三次多项式的系数, npolys 是多项式的个数, ncoefs 是每有多项式的系数的个数, dim 是样条的维数. 注意: coefs 中的第 i 个三次多项式的系数是指形如

中的系数 $a_{j}(j=0,\cdots,3)$ .

用 Matlab 进行三次样条插值的更详细内容请参看 Matlab 使用手册.

例 3.5.3 已知函数表 (表 3-6), 在区间 [1,5] 上求满足表 3-6 所给出的插值条件的三次自然样条插值函数.

表 3-6 给定数据表

解：因为这是在第二种边界条件下的插值问题，故确定 $M_{0}, M_{1}, M_{2}, M_{3}$ 的方程组形如式 (3.54)，由已知自然边界条件，有 $M_{0} = M_{3} = 0$ ，则得求解 $M_{1}, M_{2}$ 的方程组为

为确定方程组的右端项 $d_{1}, d_{2}$ ，先作差商表 (表 3-7) 得

再由 $x_0 = 1, x_1 = 2, x_2 = 4, x_3 = 5$ 得 $h_1 = x_1 - x_0 = 1, h_2 = 2, h_3 = 1$ . 所以， $\lambda_1 = \frac{h_2}{h_1 + h_2} = \frac{2}{3}, \mu_2 = \frac{h_2}{h_2 + h_3} = \frac{2}{3}$ ，代入方程组 (3.54)，得

表 3-7 例 3.5.3 的差商表

用 MatLab 完成此题的程序如下:

>> x = [1 2 4 5];
>> y = [0 1 3 4 2 0];
>> pp = csape(x,y,'second')
>> [breaks,coefs,npolys,ncoeffs,dim] = unmkpp(pp)
pp =
    form: 'pp'
    breaks: [1 2 4 5]
    coefs: [3x4 double]
    pieces: 3
    order: 4
    dim: 1
breaks =
coefs =
    -0.1250    0    2.1250    1.0000
    -0.1250   -0.3750    1.7500    3.0000
    0.3750   -1.1250   -1.2500    4.0000
npolys =
    3
ncoeffs =
    4
dim =
    1
其中 coefficients 的第 i 行为式 (3.56) 中 s(x) 的第 i 个多项式的系数

评 注#

插值是数值计算中最基础的部分, 本身又有广泛的应用. 本章除了介绍基本的算法理论外, 还重点介绍了算法产生的思想及算法的 Matlab 实现. 由于本书的主要读者对象为工科研究生, 故没有将 B- 样条作为本章的内容. 在商业软件越来越流行的今天, 数值数学教学的目的应从单纯的算法理论的教学转为注重算法思想及算法的软件实现.

习题三#

1. 设有数据表如表 3-8, 用线性插值找出 $\sin 0.705$ 和 $\cos 0.702$ 的近似值.

表 3-8 习题 1 的数据

2. 给定函数 $y = \sin x$ 的函数表如表 3-9, 试分别用线性插值与二次插值求 $\sin 0.57891$ 的近似值.

表 3-9 习题 2 的数据

3. 设给定数值表如表 3-10.

表 3-10 习题 3 的数据

(1) 构造出差商表;

(2) 用 4 次牛顿插值多项式求出 $f(4.2)$ 的近似值.

4. 多项式 $p(x)=x^{4}-x^{3}+x^{2}-x+1$ 有数值表如表 3-11.

表 3-11 习题 4 中 $p\left( x\right)$ 的数据

试找一次数不大于 5 的多项式 $q(x)$ , 它取数值表如表 3-12.

表 3-12 习题 4 中 $q(x)$ 的数据

5. 如何判断下面数值表 3-13 是否来自一个三次多项式？试作解释.

表 3-13 习题 5 中给定的数据

6. 如果用一个 20 次的插值多项式在 $[0,2]$ 上逼近 $\mathrm{e}^{-x}$ , 那么精确性如何?

7. 设用在区间 [1,2] 上 10 个均匀分布节点的 9 次插值多项式逼近函数 $f(x) = \ln x$ ，误差界是多少？

8. 设 $l_{0}(x)$ , $l_{1}(x)$ , $\cdots$ , $l_{n}(x)$ 是以 $x_{0}$ , $x_{1}$ , $\cdots$ , $x_{n}$ 为节点的 n 次拉格朗日插值基函数, 证明

9. 假设对函数 $f(x)$ 在步长为 $h$ 的等距点上造表, 且 $|f''(x)| \leqslant M$ , 证明: 在表中任意相邻两点间做线性插值, 误差不超过 $\frac{1}{8} M h^2$ . 设 $f(x) = \sin x$ , 问 $h$ 应取多大才能保证线性插值的误差不大于 $\frac{1}{2} \times 10^{-6}$ .

10. 已知 $f(x)=x^{8}+x^{5}-32$ ，求 $f[3^{0},3^{1},\cdots,3^{8}]$ ， $f[3^{0},3^{1},\cdots,3^{9}]$ .

11. 设 $l_{0}(x)$ , $l_{1}(x)$ , $\cdots$ , $l_{n}(x)$ 是以 $x_{0}$ , $x_{1}$ , $\cdots$ , $x_{n}$ 为节点的 n 次拉格朗日插值基函数, 对 n = 1 直接验证

然后对任意 n 值建立上述等式.

12. 确定 $a, b, c, d$ , 使得 $f$ 是一个三次样条函数, 且 $\int_0^2 (f''(x))^2 \mathrm{d}x$ 最小:

13. 用笔算找出表 3-14 的自然三次样条插值.

表 3-14 习题 13 中给定的三次样条数据

14. 找出节点 -1, 0, 1 上的三次样条 $s(x)$ , 使得 $s''(-1) = s''(1) = 0$ , $s(-1) = s(1) = 0$ 和 $s(0) = 1$ .

数值实验三#

1. 找出函数 $\arctan x$ 在区间 [1,6] 的 11 个等距点上插值的 10 次多项式, 打印出这个多项式的牛顿形式中的系数. 计算并打印出这个多项式与 $\arctan x$ 之差在区间 [0,8] 的 33 个等距点上的值. 由此能得出什么结论?

2. 检验正文中给出的决定差商表的程序. 例如, 计算出表 3-15 的差商, 并由此找出表 3-15 的插值多项式.

表 3-15 实验题 2 给定的数据

3. 使用区间 $[-5, 5]$ 上的 21 个等距节点, 找出函数 $f(x) = (x^2 + 1)^{-1}$ 的 20 阶插值多项式 $p(x)$ . 打印出 $f(x)$ 和 $p(x)$ 的图形, 观察 $f(x)$ 和 $p(x)$ 的最大偏差.

4. 在计算机上, 对上一题使用切比雪夫节点 $x_{i} = 5\cos (i\pi /20)$ , $0 \leqslant i \leqslant 20$ , 找出函数 $f(x) = (x^{2} + 1)^{-1}$ 的 20 阶插值多项式 $q(x)$ . 打印出 $f(x)$ 和 $q(x)$ 的图形. 由上一题和本题, 你能得出什么结论?

5. 找出函数 $f(x) = (x^2 + 1)^{-1}$ 在区间 $[-5, 5]$ 上的 41 个等距节点的三次样条函数 $s(x)$ , 打印出 $f(x)$ 和 $s(x)$ 的图形.

第4章 函数逼近#