现在我们要回答第 3 章开始时提出的第二个问题和第三个问题, 即函数的逼近问题. 函数逼近问题根据度量逼近函数 $p(x)$ 与被逼近函数 $f(x)$ 的近似程度的标准不同, 可分为一致逼近和平方逼近.

用函数 $f(x)$ 和 $p(x)$ 的最大误差 $\| p - f\|_{\infty} = \max_{a\leqslant x\leqslant b}|p(x) - f(x)|$ 作为度量逼近函数 $p(x)$ 对被逼近函数 $f(x)$ 的逼近程度.若存在一个函数序列 $p_n(x)$ , 满足 $\lim_{n\to \infty}\| p_n(x) - f(x)\|_\infty = 0,$ 则意味着序列 $\{p_n(x)\}$ 在区间 $[a,b]$ 上一致收敛到 $f(x)$ . 序列 $\{p_n(x)\}$ 对 $f(x)$ 的逼近称为一致逼近.

也可以用积分 $\| p - f\| _2 = \left(\int_a^b (p(x) - f(x))^2\mathrm{d}x\right)^{\frac{1}{2}}$ 作为度量逼近函数 $p(x)$ 与被逼近函数 $f(x)$ 的近似程度.若存在一个函数序列 $\{p_n(x)\}$ ，满足 $\lim_{n\to \infty}\| p_n(x) - f(x)\| _2 = 0,$ 则称序列 $\{p_n(x)\}$ 在区间 $[a,b]$ 上平方收敛到 $f(x)$ .序列 $\{p_n(x)\}$ 对 $f(x)$ 的逼近称为平方逼近.

数值上求函数 $f(x)$ 的一致逼近函数和平方逼近函数是困难的, 并且实际问题中人们感兴趣的是最佳一致逼近函数和最佳平方逼近函数. “最佳” 的意思是什么? 将在本章的相应的地方作介绍.

4.1 内积与正交多项式#

4.1.1 权函数和内积#

定义 4.1.1　设 $[a,b]$ 为有限或无限区间, $\rho(x)$ 是定义在 $[a,b]$ 上的非负函数, $\int_{a}^{b} x^{k} \rho(x) \, dx$ 对 $k = 0, 1, \cdots$ 都存在, 对非负的 $f(x) \in C[a,b]$ , 若 $\int_{a}^{b} f(x) \rho(x) \, dx = 0$ , 则 $f(x) \equiv 0$ , 称 $\rho(x)$ 为 $[a,b]$ 上的权函数.

x 处的权函数值 $\rho(x)$ 刻划的是点 x 在 $[a, b]$ 上所处的重要性.

常见的权函数有

注 定义 4.1.2 是 $R^{n}$ 中向量 $\boldsymbol{x} = (x_{1}, \cdots, x_{n})^{\mathrm{T}}$ 与向量 $\boldsymbol{y} = (y_{1}, \cdots, y_{n})^{\mathrm{T}}$ 的数量积 $(\boldsymbol{x}, \boldsymbol{y}) = \sum_{i=1}^{n} x_{i} y_{i}$ 定义的推广.

定义了内积的线性空间就称为内积空间, 连续函数空间 $C[a,b](C[a,b]$ 表示区间 $[a,b]$ 上连续函数的全体) 定义了内积 $(f,g)$ 后就是一个内积空间. 容易验证内积满足以下基本法则:

(1) $(f,g) = (g,f)$ ;

(2) $(c_{1}f + c_{2}g,h) = c_{1}(f,h) + c_{2}(g,h), c_{1},c_{2}$ 是常数;

(3) $(f, f) \geqslant 0$ , 并且当且仅当 $f \equiv 0$ 时, $(f, f) = 0$ .

4.1.2 正交函数系#

例 4.1.5 三角函数系 $1, \sin x, \cos x, \sin 2x, \cos 2x, \cdots$ 在 $[0, 2\pi]$ 上是正交函数系 (权 $\rho(x) \equiv 1$ ).

实际上 $(1,1) = \int_0^{2\pi}\mathrm{d}x = 2\pi$ ，而

下面给出常见的且十分重要的正交多项式.

4.1.3 勒让德多项式#

勒让德 (Legendre) 多项式的表达式为 (4.4), 它是在区间 $[-1, 1]$ 上权函数 $\rho(x) = 1$ 的正交多项式:

$P_{n}(x)$ 的首项 $x^{n}$ 的系数为 $\frac{(2n)!}{2^n(n!)^2}$ , 记

则 $\widetilde{P}_n(x)$ 是首项 $x^n$ 系数为1的勒让德多项式.

勒让德多项式有许多重要性质, 特别有以下三条.

(1) 正交性

只要令 $\varphi(x)=(x^{2}-1)^{n}$ ，则 $\varphi^{(k)}(\pm1)=0, k=0,1,\cdots,n-1$ ，且 $P_{n}(x)=\frac{1}{2^{n}n!}\varphi^{(n)}(x)$ 。设多项式 $Q(x)$ 为次数 $\leqslant n$ 的多项式，用分部积分得

当 $Q(x)$ 为次数不超过 $n - 1$ 的多项式时 $Q^{(n)}(x) = 0,$ 于是有

当 $Q(x) = P_{n}(x)$ 时 $Q^{(n)}(x) = P_n^{(n)}(x) = \frac{(2n)!}{2^n n!}$ , 于是有

这就证明了式 (4.6) 的正确性.

(2) 递推公式

其中， $P_0(x) = 1$ ， $P_{1}(x) = x$

式 (4.7) 可直接利用正交性证明. 由式 (4.7) 可得

图 4-1 给出了 $P_{0}(x)$ , $P_{1}(x)$ , $P_{2}(x)$ , $P_{3}(x)$ 的图形.

图4-1 勒让德多项式的图形， $n = 0,1,2,3$

(3) 奇偶性

4.1.4 切比雪夫多项式#

切比雪夫 (Chebyshev) 多项式定义为 (4.8), 它是在区间 $[-1,1]$ 上权函数 $\rho(x)=\frac{1}{\sqrt{1-x^{2}}}$ 的正交多项式:

若令 $x = \cos \theta$ ，则 $T_{n}(x) = \cos (n\theta), 0 \leqslant \theta \leqslant \pi,$ 这是 $T_{n}(x)$ 的参数表示. 利用三角公式可将 $\cos n\theta$ 展成 $\cos \theta$ 的一个 $n$ 次多项式, 故式(4.8)是 $x$ 的 $n$ 次多项式. 下面给出 $T_{n}(x)$ 的主要性质.

(1) 正交性

只要对积分做变换 $x = \cos\theta$ ，利用三角公式就可得到 (4.9) 的结果.

(2) 递推公式

其中 $T_{0}(x)=1, T_{1}(x)=x.$

由 $x = \cos \theta, T_{n + 1}(x) = \cos (n + 1)\theta$ ，用三角公式

则得式(4.10). 由式(4.10)可推出 $T_{2}(x)$ 到 $T_{8}(x)$ 如下：

图 4-2 给出了 $T_{0}(x)$ , $T_{1}(x)$ , $T_{2}(x)$ , $T_{3}(x)$ , $T_{4}(x)$ 的图形.

图 4-2 切比雪夫多项式的图形, n = 0, 1, 2, 3, 4

(3) 奇偶性

(4) $T_{n}(x)$ 在 $(-1,1)$ 内的 n 个零点为 $x_{k} = \cos \frac{2k-1}{2n} \pi, k = 1, 2, \cdots, n,$ 在 $[-1,1]$ 上有 $n + 1$ 个极值点 $y_{k} = \cos \frac{k}{n} \pi, k = 0, 1, \cdots, n.$

(5) $T_{n}(x)$ 的最高次幂 $x^{n}$ 的系数为 $2^{n - 1}, n \geqslant 1$ .

4.1.5 其他正交多项式#

除上述两个最常用的正交多项式外, 较重要的还有如下两个无穷区间的正交多项式.

(1) 拉盖尔 (Laguerre) 多项式

拉盖尔多项式定义为 (4.11), 它是在区间 $[0, \infty)$ 上权函数 $\rho(x) = \mathrm{e}^{-x}$ 的正交多项式:

它的递推公式为:

其中， $L_{0}(x)=1$ ， $L_{1}(x)=1-x$ 。正交性为

(2) 埃尔米特 (Hermite) 多项式

埃尔米特多项式定义为 (4.13), 它是在区间 $(- \infty, + \infty)$ 上带权函数 $\rho(x) = \mathrm{e}^{-x^2}$ 的正交多项式:

它的递推公式为

其中 $H_0(x) = 1, H_1(x) = 2x$ . 正交性为

4.2 最佳一致逼近与切比雪夫展开#

在数值计算中, 人们感兴趣的是最佳一致逼近. 我们将在本节中讨论最佳的含义, 以及一些相应的方法.

4.2.1 最佳一致逼近多项式#

(1) 一致逼近的定义

(2) 一致逼近多项式的存在性#

对于 $[a, b]$ 上的连续函数 $f(x)$ , 是否存在多项式 $p(x)$ 一致逼近 $f(x)$ 呢? 德国数学家魏尔斯特拉斯 (Weierstrass) 在 1885 年给出下述著名定理.

魏尔斯特拉斯定理从理论上证明了存在多项式 $p(x)$ 可以以任意的精度逼近闭区间上的连续函数, 同时从伯恩斯坦 (Bernstein) 的构造性证明中知道: 若给定的逼近精度高, 则找到的多项式次数就高. 在实际计算中, 人们对下面两类问题感兴趣.

问题 1: 若多项式次数 n 固定, 求一个多项式 $p_{n}^{*}(x)$ , 使 $\max_{a \leqslant x \leqslant b} |p_{n}^{*}(x) - f(x)|$ 最小. 该问题中的多项式即是最佳一致逼近多项式.

问题 2: 若给定逼近精度, 求次数较低的逼近多项式.

(3) 最佳一致逼近多项式

定义 4.2.3　记 $P_{n}$ 为次数不超过 n 的多项式全体, 设 $f(x) \in C[a, b]$ , $p(x) \in P_{n}$ , 记

称为 $p(x)$ 与 $f(x)$ 的偏差. 若存在 $x_0 \in [a, b]$ , 使

则称 $x_0$ 是 $p(x)$ 关于 $f(x)$ 的偏差点. 若 $p(x_0) - f(x_0) = \mu$ , 则称 $x_0$ 为正偏差点, 若 $p(x_0) - f(x_0) = -\mu$ , 则称 $x_0$ 为负偏差点.

可以证明最佳逼近多项式存在且唯一.

证明略. 进一步可以证明下面的重要定理.

定理 4.2.6（切比雪夫定理） $p_n^*(x) \in P_n$ 是 $f(x) \in C[a, b]$ 的最佳逼近多项式的充分必要条件是：在 $[a, b]$ 上至少有 $n + 2$ 个交替为正负的偏差点，即至少有 $n + 2$ 个点 $a \leqslant x_1 < x_2 < \cdots < x_{n+2} \leqslant b$ ，使得

上述点 $\{x_{k}\}_{1}^{n+2}$ 称为切比雪夫交错点组.

定理证明可见参考书目 [11]. 这个定理给出了最佳逼近多项式的基本特性, 它是求最佳逼近多项式的主要依据, 但最佳逼近多项式的计算是困难的, 下面只对 $n = 1$ 的情形进行讨论.

4.2.2 线性最佳逼近多项式的求法#

本节讨论问题：设 $f(x)$ 在 $[a,b]$ 上有二阶导数，且 $f''(x)$ 在 $[a,b]$ 上不变号，求 $f(x)$ 的线性最佳一致逼近多项式 $p_{1}(x)=a_{0}+a_{1}x$ .

由于 $f(x) \in C[a, b]$ 且 $f''(x)$ 在 $[a, b]$ 上不变号, 由定理4.2.6可知, 存在点 $a \leqslant x_1 < x_2 < x_3 \leqslant b$ 使

由于 $f''(x) \neq 0$ , 故 $f'(x)$ 在 $[a, b]$ 上单调, 于是 $f'(x) - p_1'(x) = f'(x) - a_1 = 0$ 在 $[a, b]$ 上只有一个根 $x_2$ , 故 $p_1(x)$ 对 $f(x)$ 的另两个偏差点只能在 $[a, b]$ 的端点, 故有 $x_1 = a$ , $x_3 = b$ . 由此可得 $p_1(a) - f(a) = -[p_1(x_2) - f(x_2)] = p_1(b) - f(b)$ 或

由式 (4.18) 得

代入式 (4.19) 得

其中 $x_{2}$ 由 $f^{\prime}(x_{2}) = a_{1}$ 得到，由此则得 $f(x)$ 在 $[a,b]$ 上的最佳一次逼近多项式 $p_1(x) = a_0 + a_1x$ 。其几何意义如图4-3所示。

图 4-3 最佳一次逼近多项式的几何意义

例 4.2.7 设 $f(x)=\sqrt{1+x^{2}}$ ，在 [0,1] 上求 $f(x)$ 的最佳一次逼近多项式 $p_{1}(x)=a_{0}+a_{1}x$ . 解：由式(4.20)求得 $a_{1}=\frac{\sqrt{2}-1}{1}\approx0.414$ ，由 $f'(x)=\frac{x}{\sqrt{1+x^{2}}}$ ，得到 $f'(x_{2})=\frac{x_{2}}{\sqrt{1+x_{2}^{2}}}=0.414$ ，于是

由式 (4.21) 可求得 $a_0 = 0.955$ , 于是得 $f(x)$ 的最佳一次逼近多项式为 $p_1(x) = 0.955 + 0.414x$ .

4.2.3 切比雪夫展开与近似最佳逼近多项式#

求连续函数的最佳逼近多项式是一件困难的事, 在实际计算中, 人们常用切比雪夫多项式的性质, 求近似最佳逼近多项式. 由切比雪夫多项式 $T_{n}(x)$ 的性质 (5) 知, $T_{n}(x)$ 的最高次幂 $x^{n}$ 的系数为 $2^{n-1} (n \geqslant 1)$ , 于是 $\widetilde{T_{n}}(x) = \frac{1}{2^{n-1}} T_{n}(x)$ 是最高项系数为 1 的 $n$ 次多项式. 对于多项式 $\widetilde{T_{n}}(x)$ 有如下结论.

证：由于 $\widetilde{T_{n}}(x)=\frac{1}{2^{n-1}}T_{n}(x)=x^{n}-p_{n-1}^{*}(x)$ ，当 $x_{k}=\cos\frac{k}{n}\pi(k=0,1,\cdots,n)$ 时，有

这表明 $p_{n-1}^{*}(x)$ 与 $f(x) = x^n$ 有 $n+1$ 个交替为正负的偏差点, 根据定理4.2.6可知, $p_{n-1}^{*}(x) \in P_{n-1}$ 是 $f(x) = x^n$ 的最佳逼近多项式, 即

所以 $\tilde{T}_{n}(x)$ 是 $[-1,1]$ 上与零偏差最小的多项式.

这个定理给出了切比雪夫多项式的一个非常重要的性质, 它表明以 $\widetilde{T}_{n}(x)$ 为余项的误差在整个区间 $[-1,1]$ 上的分布是均匀的, 因此这一性质在求函数的近似最佳逼近多项式中被广泛应用.

例 4.2.9 设 $f(x)=x^{4}$ ，在 $[-1,1]$ 上求 $f(x)$ 在 $P_{3}$ 中的最佳逼近多项式.

解：根据定理 4.2.6, 已知 $\widetilde{T}_{4}(x)=\frac{1}{8}\left(8x^{4}-8x^{2}+1\right)=x^{4}-\left(x^{2}-\frac{1}{8}\right)$ ，于是可知所求最佳逼近多项式为 $p_{3}^{*}(x)=\left(x^{2}-\frac{1}{8}\right)$ .

根据定理 4.2.6, 实际上只要误差 $f(x) - p_{n-1}^{*}(x) = aT_{n}(x)$ ，则 $p_{n-1}^{*}(x) \in P_{n-1}$ 是区间 $[-1,1]$ 上多项式 $f(x)$ 的最佳一致逼近多项式，因此有 $n+1$ 个交替为正负的偏差点。更一般地，若在 $[-1,1]$ 上 $f(x) - p_{n}(x) \approx aT_{n+1}(x)$ ，那么 $p_{n}(x) \in P_{n}$ 可作为 $f(x)$ 在 $P_{n}$ 中的近似最佳逼近多项式。例如，设 $f(x) \in C[-1,1]$ 的拉格朗日插值多项式为 $L_{n}(x)$ ，其余项可表示为

若取 $\omega_{n+1}(x)=(x-x_{0})(x-x_{1})\cdots(x-x_{n})=\widetilde{T}_{n+1}(x)$ . 这时插值节点为 $\widetilde{T}_{n+1}(x)$ 的 $n+1$ 个零点

这样构造的插值多项式 $L_{n}(x)$ 也可作为最佳逼近多项式的近似.

更常用和可行的方法是将 $f(x) \in C[a, b]$ 直接按 $\{T_k(x)\}_0^\infty$ 展开并用它的前 $n + 1$ 项部分和逼近 $f(x)$ , 由于 $\{T_k(x)\}_0^\infty$ 是在区间 $[-1, 1]$ 上带权 $\rho(x) = (1 - x^2)^{-1/2}$ 正交的多项式, 构造 $n$ 次多项式 $c_n^*(x)$ , 其表达式为

其中

可以证明, 如果 $f(x) \in C[a, b]$ , 则 $\lim_{n \to \infty} \|f(x) - c_n^*(x)\|_\infty = 0$ , 由此可得

这表明用 $c_{n}^{*}(x)$ 逼近 $f(x)$ 的误差近似于 $a_{n + 1}^{*}T_{n + 1}(x)$ , 是均匀分布的, 故 $c_{n}^{*}(x)$ 是 $f(x)$ 在 $[-1,1]$ 上近似最佳逼近多项式.

例 4.2.10 设 $f(x)=\mathrm{e}^{x}$ ，在 $[-1,1]$ 上按切比雪夫多项式 $\{T_{k}(x)\}_{0}^{\infty}$ 展开，并计算到 n=3. 解：由式 (4.24) 可得出

此积分可用数值积分计算得表 4-1.

表 4-1 切比雪夫展开多项式系数

于是可得到

直接计算可得

这与 $f(x) = \mathrm{e}^{x}$ 的最佳一致逼近多项式 $p_3^* (x)\in P_3$ 的误差

相差很小, 而计算 $p_3^*(x)$ 通常是相当困难的, 它是根据定理 4.2.6 用逐次迭代求得, 但计算 $c_3^*(x)$ 则容易得多了, 并且效果很好.

4.3 最佳平方逼近#

4.3.1 预备知识#

为讨论最佳平方逼近函数的求法, 我们回忆如下两个线性代数知识.

(1) 设 $\varphi_{0}(x),\varphi_{1}(x),\cdots,\varphi_{n}(x)$ 是 $[a,b]$ 上的线性无关的连续函数, $a_{0},a_{1},\cdots,a_{n}$ 是任意实数, 则 $\Phi=\{s(x)\mid s(x)=a_{0}\varphi_{0}(x)+\cdots+a_{n}\varphi_{n}(x),a_{i}\in\mathbb{R},i=0,1,\cdots,n\}$ , 并称集合 $\Phi$ 是由 $\varphi_{0}(x),\varphi_{1}(x),\cdots,\varphi_{n}(x)$ 所生成的线性空间, $\varphi_{0}(x),\varphi_{1}(x),\cdots,\varphi_{n}(x)$ 是生成空间 $\Phi$ 的一个基底.

(2) 设 $\varphi_{0}(x),\varphi_{1}(x),\cdots,\varphi_{n}(x)\in C[a,b]$ ，则 $\varphi_{0}(x),\varphi_{1}(x),\cdots,\varphi_{n}(x)$ 在 $[a,b]$ 上线性无关的充分必要条件是 $\det G_{n}\neq0$ ，其中

这里 $(\cdot,\cdot)$ 表示内积.

证：由于 $\operatorname{det} G_n \neq 0$ 等价于齐次方程组

只有零解.

当 $\varphi_{0}(x),\varphi_{1}(x),\cdots,\varphi_{n}(x)$ 线性无关时, 有 $\det G_{n}\neq0$ . 否则的话, 存在非零向量 $\boldsymbol{c}=(c_{0},c_{1},\cdots,c_{n})^{\mathrm{T}}\in\mathbb{R}^{n+1}$ 使得 $G_{n}c=0$ , 从而 $(\sum_{i=0}^{n}c_{i}\varphi_{i}(x),\sum_{i=0}^{n}c_{i}\varphi_{i}(x))=\boldsymbol{c}^{\mathrm{T}}G_{n}\boldsymbol{c}=0$ . 因此, $\sum_{i=0}^{n}c_{i}\varphi_{i}(x)=0$ 有非零解, 与假设矛盾.

反之, 若 $\operatorname{det} G_n \neq 0$ , 则由 $\sum_{i=0}^{n} c_i \varphi_i(x) = 0$ , 可得

故 $c_{0}=c_{1}=\cdots=c_{n}=0$ ，由线性无关的定义知 $\varphi_{0}(x),\varphi_{1}(x),\cdots,\varphi_{n}(x)$ 线性无关.

4.3.2 最佳平方逼近#

若 $\Phi = P_{n} = \operatorname{span}\{1, x, \cdots, x^{n}\}$ ，则满足上述定义的 $s^{*}(x)$ 是 $f(x)$ 的 n 次最佳平方逼近多项式.

那么最佳平方逼近函数是否存在呢？下面的定理回答了这个问题.

先构造函数 $s^{*}(x)$ . 由定义知, 求 $s^{*}(x) \in \Phi$ 等价于求多元函数 $I(a_{0}, a_{1}, \cdots, a_{n})$ 的极小值, 其中 $I(a_{0}, a_{1}, \cdots, a_{n}) = \int_{a}^{b} \rho(x) \left( \sum_{j=0}^{n} a_{j} \varphi_{j}(x) - f(x) \right)^{2} \mathrm{d}x$ 是关于 $a_{0}, a_{1}, \cdots, a_{n}$ 的二次函数. 由多元函数取极值的必要条件得

于是有

这是关于 $a_{0}, a_{1}, \cdots, a_{n}$ 的线性方程组, 称为法方程. 由于 $\varphi_{0}(x), \varphi_{1}(x), \cdots, \varphi_{n}(x)$ 线性无关, 由预备知识中的定理知, 方程组 (4.26) 的系数矩阵 $G_{n}$ 非奇异, 故此方程组有唯一解 $a_{k} = a_{k}^{*}, k = 0, 1, \cdots, n$ , 于是有

下面证明由此得到的 $s^{*}(x)$ 满足式 (4.25)，即对任何 $s(x) \in \Phi$ 有

由于 $a_{j}^{*}(j=0,1,\cdots,n)$ 是式 (4.26) 的解, 因此有

于是对任何 $s(x) \in \Phi$ , 有 $(f - s^{*}, s) = 0$ 及 $(f - s^{*}, s^{*} - s) = 0$ , 由此可得

这就证明了式 (4.28). 从而也证明了 $f(x)$ 在 $\Phi$ 中存在唯一的最佳平方逼近函数 $s^*(x)$ . 此时作为特例, 若取 $\varphi_k(x) = x^k$ , $k = 0,1,\dots,n$ , 区间取为 $[0,1]$ , $\rho(x) = 1$ , $f(x) \in C[0,1]$ , 在 $\Phi = P_n = \text{span}\{1,x,\dots,x^n\}$ 上的最佳平方逼近多项式为

此时由于

相应于法方程 (4.26) 的系数矩阵 $G_{n}$ 记为 $H_{n+1}=G(1,x,\cdots,x^{n})$ ，即

其中 $h_{ij} = \frac{1}{i + j - 1}$ . $\pmb{H}_{n + 1}$ 称为希尔伯特 (Hilbert) 矩阵.

若记 $\boldsymbol{a}=(a_{0},a_{1},\cdots,a_{n})^{\mathrm{T}},\boldsymbol{d}=(d_{0},d_{1},\cdots,d_{n})^{\mathrm{T}}$ ，此时法方程为

它的解为 $a_{k}=a_{k}^{*}, k=0,1,\cdots,n,$ 由此则得最佳平方逼近多项式 $p_{n}^{*}(x)$ .

例4.3.3 设 $f(x) = \sqrt{1 + x^2}$ ，求[0,1]上的一次最佳平方逼近多项式 $p_1^*(x) = a_0^* + a_1^* x$ . 解：

法方程为

求得解为 $a_0^* = 0.934, a_1^* = 0.426$ , 于是得

由于 $H_{n+1}$ 是病态矩阵, 在 $n \geqslant 3$ 时直接解法方程 (4.30) 的误差很大, 因此当 $\varphi_{k}(x) = x^{k} (k = 0, 1, \cdots, n)$ 时, 解法方程方法只适合 $n \leqslant 2$ 的情形. 对 $n \geqslant 3$ 的情形, 可用正交多项式作 $\Phi$ 的基的方法来求解最佳平方逼近多项式.

设 $f(x) \in C[a, b]$ , $\Phi = \operatorname{span}\{\varphi_{0}, \varphi_{1}, \cdots, \varphi_{n}\}$ , 若 $\varphi_{0}(x), \varphi_{1}(x), \cdots, \varphi_{n}(x)$ 是正交函数系，则当 $i \neq j$ 时， $(\varphi_{i}, \varphi_{j}) = 0$ ，而 $(\varphi_{j}, \varphi_{j}) > 0$ ，于是法方程 (4.26) 的系数矩阵 $G_{n}$ 为非奇异对角阵，方程的解为

于是 $f(x)$ 在 $\Phi$ 中的最佳平方逼近函数为

称式 (4.32) 为 $f(x)$ 的广义傅里叶 (Fourier) 展开, 相应式 (4.31) 的系数 $a_k^*$ 称为 $f(x)$ 的广义傅里叶系数.

下面考虑特殊情形, 设 $[a,b]=[-1,1]$ , $\rho(x)=1$ , 此时正交多项式为勒让德多项式, 取 $\Phi=\operatorname{span}\left\{P_{0}(x),P_{1}(x),\cdots,P_{n}(x)\right\}$ , 根据式 (4.32) 可得 $f(x)$ 的最佳平方逼近多项式为

其中

这样得到的最佳平方逼近多项式 $s_{n}^{*}(x)$ 与直接由 $\{1, x, \cdots, x^{n}\}$ 为基得到的 $p_{n}^{*}(x)$ 是一致的，但此处不用解病态的法方程 (4.30). 关于病态的概念, 参看 §6.1 节.

例4.3.4 用勒让德正交多项式求 $f(x) = \mathrm{e}^{x}$ 在 $[-1,1]$ 上的最佳平方逼近多项式 (取 $n = 1,3$ ).

解：先计算

由 (4.34) 可算出

于是由 (4.33) 可求得

如果所给的区间不是 $[-1, 1]$ , 而是一般的有限区间 $[a, b]$ , 那么, 可以通过变量替换 $x = \frac{a + b}{2} + \frac{b - a}{2} t$ 将它转化为区间 $[-1, 1]$ 上的情况来处理.

例4.3.5 求函数 $f(x) = \sqrt{x}$ 在[0,1]上的一次最佳平方逼近多项式.

解：先作变换 $x = \frac{1}{2} (1 + t)$ ，则

下面求 $g(t)$ 在区间 $[-1, 1]$ 上的一次最佳平方逼近多项式 $q_{1}(t)$ .

由

可知

再把 $t = 2x - 1$ 代入 $q_{1}(t)$ 就得到 $f(x) = \sqrt{x}$ 在[0,1]上的一次最佳平方逼近多项式为 $s^{*}(x) = \frac{2}{3} +\frac{6}{15} (2x - 1) = \frac{4}{15} +\frac{4}{5} x.$

4.4 曲线拟合的最小二乘法#

4.4.1 最小二乘法#

若 $f(x)$ 是由实验或观测得到的，则其函数通常由函数表 $(x_{i}, f(x_{i})) (i = 0, 1, \cdots, m)$ 给出。由前面叙述的内容，由函数表给出函数关系通常有下列两个方法：方法一，使用多项式插值；方法二，三次样条插值。使用多项式插值会带来两个问题，问题之一，当所给的数值点较多时，多项式次数要高，会出现数值震荡，即所谓的龙格现象；问题之二，由于数值本身带有误差，使用插值条件来确定函数关系不合理。三次样条插值克服了多项式插值的第一个缺点，但求三次样条插值带来了大的计算量。曲线拟合的最小二乘方法可以克服数值震荡，同时不引起大的计算量。那么，什么是曲线拟合的最小二乘方法呢？

对函数表在函数空间 $\Phi$ 中求 $s^{*}(x)$ ，使

就是曲线拟合的最小二乘问题, 其中 $\omega_{i}$ 是点 $x_{i}$ 处的权. 这个问题的实质是 $f(x)$ 为离散情形的最佳平方逼近问题. 这样求连续函数的最佳平方逼近的方法可以用到曲线拟合的最小二乘问题上来. 具体做法简述如下.

求 $s^{*}(x)$ 的问题等价于求多元函数

的极小值, 由取极值的必要条件可得法方程

其中内积 $(\cdot,\cdot)$ 表示和式, 即

由于 $\varphi_{0}(x),\varphi_{1}(x),\cdots,\varphi_{n}(x)$ 线性无关, 可以证明法方程 (4.36) 的系数矩阵非奇异, 于是求解法方程得 $a_{k}=a_{k}^{*}, k=0,1,\cdots,n$ , 从而得到 $s^{*}(x)$ , 它是存在的且唯一的. 称 $s^{*}(x)$ 是 $f(x)$ 在 $\Phi$ 中的最小二乘逼近函数.

在最小二乘逼近中如何选择数学模型是很重要的, 即如何根据给定的 $f(x)$ 来选择 $\Phi$ . 通常要根据物理定义或 $f(x_{i})(i=0,1,\cdots,n)$ 数据分布的大致图形选择相应的数学模型.

例 4.4.1 给定数据 (见表 4-2), 试选择适当模型, 求最小二乘拟合函数 $s^{*}(x)$ .

表 4-2 给定的数据表

解：在坐标平面上描出数据表中的点 $(x_{i},y_{i})(i=0,1,\cdots,9)$ ，根据点的分布情况，所求函数可用 $y=\ln x,y=\cos x$ 和 $y=e^{x}$ 的线性组合表示，因此选拟合函数为

其中 a, b, c 为待定参数, 基底为

利用公式 (4.36) 和 (4.37), 建立法方程组如下:

其解为

所以

为所要求的最小二乘解.

此题若借助 Matlab 求解, 程序如下:

则 $a = Z(1)$ , $b = Z(2)$ , $c = Z(3)$ .

读者运行上述程序后, 将会发现 Matlab 的结果与上述计算的结果一致, 但显然 Matlab 解决问题更加方便、简洁.

作为曲线拟合的一种常见情况, 一般讨论代数多项式拟合, 即取

那么相应的法方程 (4.36) 就是

此时 $s^{*}(x) = \sum_{k=0}^{n} a_{k}^{*} x^{k}$ , 称它为数据拟合多项式, 上述拟合称为多项式拟合.

例 4.4.2 求数据表 (表 4-3) 的二次最小二乘拟合多项式.

表 4-3 二次最小二乘拟合数据

解：设二次拟合多项式为

将数据表代入式 (4.28), 得此问题的法方程为

该方程组的解为

所以, 此数据组的最小二乘拟合二次多项式为

注意：此问题没有给出各点的权函数值 $\omega_{i}$ 。一般情形下，若没有给出各点权函数值 $\omega_{i}$ ，视各点权函数值为 1。

此题借助 Matlab 求解, 有两个方法.

方法一：

>> x = [0 0.25 0.50 0.75 1.00 ]';
>> y = [1.00 1.284 1.6487 2.1170 2.7183]';
>> x1 = ones(5,1);
>> A = [x1 x x.^2];
>> Z = A \ y

此时 $a_{0}=Z(1)$ , $a_{1}=Z(2)$ , $a_{2}=Z(3)$ .

方法二：

此时 $a_{2}=p(1)$ , $a_{1}=p(2)$ , $a_{0}=p(3)$ .

函数 polyfit 是最小二乘多项式拟合函数, 它的用法如下: p=polyfit(x,y,n), 其中 x,y 为要拟合的数据; n 为要拟合的多项式的阶数; p 为一个向量, 各个分量是所求的拟合多项式的幂次从高到低排列的系数. 当要拟合的多项式的阶数等于给定数据点的个数减 1 时, 函数 polyfit 所得的结果即为插值多项式的系数.

实际计算与理论分析表明, 用多项式作最小二乘的基函数, 当 n 较大时, 方程组 (4.38) 的解对初始数据的微小变化非常敏感, 属于 “病态” 问题.

4.4.2 利用正交多项式作最小二乘拟合#

在用多项式作拟合函数时, 为避开求解法方程组 (4.38), 考虑选择正交多项式做基底.

若已给数据 $(x_{i},f_{i})$ 及权 $\omega_{i}(i=0,1,\cdots,m)$ ，可以构造出带权 $\{\omega_{i}\}_{i=0}^{m}$ 的正交多项式，用递推关系可表示为

其中 $\alpha_{k},\beta_{k - 1}$ 为

当取关于点集 $\{x_{i}\}_{i=0}^{m}$ 带权 $\{\omega_{i}\}_{i=0}^{m}$ 的正交的多项式族 $p_{0}(x), p_{1}(x), \cdots, p_{n}(x)$ 作为拟合多项式的基底时, 法方程组 (4.38) 便化简为

求解可得

则 $n$ 次多项式

为拟合给定的数据点 $\{(x_i, y_i)\}_{i=0}^m$ 的最小二乘解.

利用正交多项式作最小二乘拟合时, 可以将构造正交多项式族 $\{p_{k}(x)\}_{k=0}^{n}$ 与解法方程组求 $a_{k}$ 以及形成拟合多项式 $g_{n}(x)$ 穿插进行, 见下面的例题.

例 4.4.4 给定数据点 $(x_{i}, y_{i})$ 及 $\omega_{i}$ 如表 4-4. 试用最小二乘法求拟合这组数据的多项式.

表 4-4 多项式拟合数据表

解: 为确定拟合多项式的次数, 首先描点, 如图 4-4. 根据数据点的分布情况, 用二次多项式拟合这组数据. 函数类的基底取正交多项式 $p_0(x), p_1(x), p_2(x)$ , 所以拟合函数为 $g(x) = a_0 p_0(x) + a_1 p_1(x) + a_2 p_2(x)$ . 以下计算一律取 6 位有效数字.

图 4-4 例 4.4.4 数据点的分布情况

取 $p_{0}(x)=1$ ，则 $p_{0}(x_{i})=1, i=0,1,\cdots,6.$

则

从而拟合多项式为

用 Matlab 工具解决此问题的方法如下:

$$

\mathrm{x} = [ 0 \quad 0.5 \quad 0.6 \quad 0.7 \quad 0.8 \quad 0.9 \quad 1.0 ]; $$

则运行该程序的结果为

显然借助 Matlab 工具, 求数据的最小二乘拟合多项式非常简单.

4.4.3 非线性最小二乘问题#

有时需要用指数函数 $s(x)=ae^{bx}$ 、幂函数 $g(x)=ax^{b}$ 或三角函数 $h(x)=a\sin bx$ 等非多项式函数拟合给定的一组数据, 这时拟合函数是关于待定参数的非线性函数, 按最小二乘准则, 用极值原理建立的法方程组将是关于待定参数的非线性方程组, 称这类数据拟合问题为非线性最小二乘问题, 其中某些特殊的情形可以转化为线性最小二乘问题求解.

例 4.4.5 设给定数据 $(x_{i}, f_{i})(i = 0, 1, 2, 3, 4)$ (表 4-5)，试用最小二乘法求拟合这组数据的函数.

表 4-5 未知函数数据表

解：所给数据接近一个指数曲线, 因而可选择数学模型为指数函数 $y = a\mathrm{e}^{bx}$ ( $a, b$ 为待定常数). 这是一个关于 $a, b$ 的非线性模型, 为此对 $y = a\mathrm{e}^{bx}$ 两边取对数得 $\ln y = \ln a + bx$ , 令 $u = \ln y$ , $A = \ln a$ , 于是有 $u = A + bx$ , 这是一个线性模型, 可用最小二乘求解.

取 $\varphi_0 = 1, \varphi_1 = x,$ 要求 $u = A + bx$ 与 $(x_i, u_i)(i = 0, 1, 2, 3, 4)$ 做最小二乘拟合，由式(4.38)得法方程

解得 A = 1.122, b = 0.5056, 从而 $a = e^{A} = 3.071$ , 于是得最小二乘拟合曲线方程为

对某些问题, 同时含有线性理论和非线性理论的因素, 这时不能将非线性问题转化为线性最小二乘问题, 而只能按最小二乘原则, 用极值原理建立法方程组. 这里得到的法方程组将是关于待定参数的非线性方程组. 用合适的求解非线性方程组的方法求解即可得非线性最小二乘问题的解. 下面是非线性最小二乘问题的简单例子.

例4.4.6 用函数 $y = a\sin bx$ 拟合数据 (表4-6).

表 4-6 函数 $y = a\sin{bx}$ 拟合数据

解：这是一个非线性最小二乘问题, 按照最小二乘原理, 应选取参数 a, b 使得表达式 $I = \sum_{i=1}^{8}(a \sin bx_i - y_i)^2$ 达到极小值. 通过对 I 关于 a 和 b 求偏导数, 并置这些偏导数等于 0 得

[\left\{ \begin{array} { l } { \frac { \partial I } { \partial a } =  \sum _ { i = 1 } ^ { 8 } 2 ( a \sin b x _ { i } - y _ { i } ) \sin b x _ { i } \qquad = 0 , } \\ { \frac { \partial I } { \partial b } =  \sum _ { i = 1 } ^ { 8 } 2 ( a \sin b x _ { i } - y _ { i } ) a ( \cos b x _ { i } ) x _ { i } = 0 . } \end{array} \right.]

从上述两个方程中分别解出 a 并置这两个值相等, 得方程

再从这个方程中用非线性方程求根的数值方法 (如弦截法) 解出参数 $b$ . 最后可以计算这个方程的任一边作为 $a$ 的值.

不借助计算机, 上述问题得不到具体的 $a$ 和 $b$ 的值. 用 Matlab 解决此问题时, 需要借助 fminsearch 函数 (求多元函数的极小值) 来求 $I = \sum_{i=1}^{8}(a\sin bx_i - y_i)^2$ 的极小值. 具体做法分两步: 第一步创建函数 fitfun, fitfun 实际上返回 $I$ 的值; 第二步用 fminsearch 求 fitfun 的极小值. fminsearch 函数的使用方法如下:

>> [xmin, value, flag, output] = fminsearch(f, x0)

其中 f 是向量参数 x 的标量函数, x0 是搜索开始的向量. 输出参数有 4 个: 最小值出现的点 xmin, 在最小值点的函数值 value, 一个表明运行成功的标志符 flag, 以及一个算法统计结构 output.

函数 fitfun 以向量 $(a,b)^{\mathrm{T}}$ 作为自变量, 以计算的 x 和 y 值为参数值, 返回误差.

nction [err,a,b] = nlfit(x,y)
    if nargin<2,
    x = [1:8]'/10;
    y = [0.6 1.1 1.6 1.8 2.0 1.9 1.7 1.3]';
end
c = fminsearch(@fitfun,[0;0],optimset,x,y);
fprintf('The nonlinear least square fitting y=a*sin(b*x) for data\n\n');
fprintf('%6.1f',x);

fprintf('\n');
fprintf('%6.1f',y);
fprintf('\n\n is\n\t y = %7.4f * sin( %7.4f * x )\n\n',c(1),c(2));
z = linspace(x(1),x(end),100);
plot(x,y,'r+',z,c(1)*sin(c(2)*z),'b-.')

function err = fitfun(c,x,y)
    a = c(1); % coefficients
    b = c(2);
    err = y - a * sin(b*x);
    err = err'*err;

运算结果如下:

>> nlfit
The nonlinear least square fitting y=a*sin(b*x) for data
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8
0.6 1.1 1.6 1.8 2.0 1.9 1.7 1.3

is
y = 1.9751 * sin(3.0250 * x)

图 4-5 就是利用 fminsearch 进行非线性最小二乘拟合的结果, '+' 为原始数据点, 虚线为拟合得到的曲线 $y = a \sin bx = 1.9751 \sin 3.0250x$ , 由图 4-5 可以看出, 它的拟合结果是相当好的.

图 4-5 非线性最小二乘拟合

也可以用 Matlab 函数 nlinfit(非线性拟合函数), 具体程序如下:

function [err,a,b] = nlfitb(x,y)
    if nargin<2,
    x = [1:8]'/10;
    y = [0.6 1.1 1.6 1.8 2.0 1.9 1.7 1.3]';

end
beta0 = [1 1]';
beta = nlinfit(x,y,@mymodel,beta0);
fprintf('The nonlinear least square fitting y=a*sin(b*x) for data\n\n');
fprintf('%6.1f',x);
fprintf('\n');
fprintf('%6.1f',y);
fprintf('\n is\n\t y = %7.4f * sin( %7.4f * x )\n\n',beta);
z = linspace(x(1),x(end),100);
plot(x,y,'ro',z,beta(1)*sin(beta(2)*z),'b-.')

function yb = mymodel(beta, xb)
    yb = beta(1) * sin(beta(2) * xb);

4.4.4 矛盾方程组#

在实际问题中, 人们会碰到矛盾方程组, 即形式如下的方程组

如果将给定的 n 元组 $(x_{1}, x_{2}, \cdots, x_{n})$ 代入式 (4.44) 的左边, 第 k 个方程的两边之差称为第 k 个残差. 理想中所有的残差都应该是 0. 如果不可能取到 $(x_{1}, x_{2}, \cdots, x_{n})$ 使所有的残差都是 0, 则方程组 (4.44) 为矛盾的或不相容的. 在这种情况下, 替而代之的是要求残差的平方和极小. 于是, 我们要做的就是适当选取 $(x_{1}, x_{2}, \cdots, x_{n})$ , 使表达式

取到极小值, 对 $x_{i}$ 取偏导数, 并令其等于 0, 就得到法方程

这是 n 个未知数 $x_{1}, x_{2}, \cdots, x_{n}$ 的 n 个方程的线性方程组. 可以证明, 如果原始系数矩阵中的列向量是线性无关的, 则这个方程组是相容的, 从而可以解出方程组 (4.45).

4.5 周期函数逼近与快速傅里叶变换#

4.5.1 周期函数的最佳平方逼近#

当 f 为周期函数时, 用三角多项式逼近比用代数多项式更合适. 现假定 $f(x) \in C(-\infty, +\infty)$ , 并且 $f(x + 2\pi) = f(x)$ , 在 $\Phi = \operatorname{span}\{1, \cos x, \sin x, \cdots, \cos nx, \sin nx\}$ 上求最佳平方逼近函数

由于函数族 $\{1,\cos x,\sin x,\cdots,\cos nx,\sin nx\}$ 在 $[0,2\pi]$ 上是正交函数族 (见例 4.1.5)，因此由 (4.31) 可得， $f(x)$ 在 $[0,2\pi]$ 上的最佳平方逼近函数 (4.46) 中的系数为

$a_{j}, b_{j}$ 称为 $f(x)$ 的傅里叶系数. 借助最佳平方逼近的性质可得下列贝塞尔 (Bessel) 不等式

由于右边不依赖 $n$ ，故正项级数 $\frac{1}{2} a_0^2 +\sum_{k = 1}^{\infty}(a_k^2 +b_k^2)$ 收敛，并有 $\lim_{k\to \infty}a_k = \lim_{k\to \infty}b_k = 0.$ 显然三角多项式(4.46)是 $f(x)$ 的傅里叶级数

的部分和. 当 $f(x)$ 连续并以 $2\pi$ 为周期, 则由系数 (4.47) 定义的 $s_n(x)$ 平方收敛于 $f(x)$ . 若还假定 $f(x)$ 在 $x$ 处可微, 则 $s_n(x)$ 在 $x$ 处收敛到 $f(x)$ .

实际问题中, 有时 $f(x)$ 仅在离散点集 $\left\{x_{k}=\frac{2\pi}{N}k\right\}_{0}^{N-1}$ 上给出函数值 $f\left(\frac{2\pi}{N}k\right), k=0,1,\cdots,N-1$ . 可以证明, 当 $2n+1\leqslant N$ 时, 三角函数族 $\{1,\cos x,\sin x,\cdots,\cos nx,\sin nx\}$ 为离散点集 $\left\{x_{k}=\frac{2\pi}{N}k\right\}_{0}^{N-1}$ 的正交函数族, 即对任何 $k,l=0,1,\cdots,n$ 有

于是由离散点集 $\left\{x_{k} = \frac{2\pi}{N} k\right\}_{0}^{N - 1}$ 给出的 $f(x)$ 在三角函数族

中的最小二乘解, 仍可用 (4.46) 中的 $s_{n}(x)$ 表示, 其中系数为

这里 $2n + 1 \leqslant N$ . 当 $2n + 1 = N$ 时, 则有

此时 $s_n(x)$ 就是三角插值多项式.

更一般的情形, 假定 $f(x)$ 是以 $2\pi$ 为周期的复值函数, 在 N 个节点 $\left\{x_{k}=\frac{2\pi}{N}k\right\}_{0}^{N-1}$ 上的函数值 $f\left(\frac{2\pi}{N}k\right)$ 已知, 令 $\psi_{k}(x)=\mathrm{e}^{ikx}=\cos kx+i\sin kx, i=\sqrt{-1}$ 是虚根, $k=0,1,\cdots,N-1$ , 则 $\{\psi_{k}\}_{0}^{N-1}$ 关于节点集 $\{x_{j}\}_{0}^{N-1}$ 正交, 即

因此 $f(x)$ 在点集 $\left\{x_{k} = \frac{2\pi}{N} k\right\}_{0}^{N - 1}$ 上的最小二乘解为

其中，

如果取 n = N - 1，则 $s_{n}(x)$ 为 $f(x)$ 在点 $x_{j}$ ( $j = 0, 1, \cdots, N - 1$ ) 上的插值函数，即有 $s_{n}(x_{j}) = f(x_{j}), j = 0, 1, \cdots, N - 1$ ，利用 (4.52) 有

(4.53) 是由 $\{f(x_j)\}_{0}^{N - 1}$ 求 $\{C_k\}_{0}^{N - 1}$ 的过程, 称为 $f$ 的离散傅里叶变换 (Discrete Fourier Transformation), 简称 DFT, 而 (4.54) 是由 $\{C_k\}_{0}^{N - 1}$ 求 $\{f(x_j)\}_{0}^{N - 1}$ 的过程, 称为 DFT 的逆变换. 它们是用计算机进行频谱分析的主要方法, 在数字信号处理、全息技术、光谱和声谱分析等领域都有广泛的应用.

4.5.2 快速傅里叶变换 (FFT)#

由 $\{f(x_j)\}_{0}^{N - 1}$ 求 $\{C_k\}_{0}^{N - 1}$ 的(4.53)式可以当作傅里叶变换

的离散近似, 其逆变换是 (4.54), 这两个表达式的计算在傅里叶分析中非常重要, 它们都可以归结为计算

其中 $\{B_{k}\}_{0}^{N-1}$ 是给定的复数序列, $W = e^{-i\frac{2\pi}{N}}$ (正变换) 或 $W = e^{i\frac{2\pi}{N}}$ (逆变换). 由 (4.55) 看到, 直接计算一个 $C_{j}$ 需要 N 次复数乘法和 N-1 次复数加法, 计算全部 $\{C_{j}\}_{0}^{N-1}$ 需要 $N^{2}$ 次复数乘法和 $N(N-1)$ 次复数加法. 当 N 很大时, 计算 $N^{2}$ 个乘法很费时, 很多问题用高速计算机也无法完成. 20 世纪 60 年代中期产生了 FFT 算法, 减少了计算量, 才使 DFT 的计算得以实现. FFT 的思想是尽量减少 (4.55) 中乘法的次数. 由于 W 是 N 等分复平面单位圆上的一点, 且 $W^{N} = 1$ , 所以 $\{W^{jk}\}_{j,k=0}^{N-1}$ 实际上仍是单位圆上的 N 个点, 用 N 去除 jk, 可得 $jk = qN + r (0 \leqslant r \leqslant N - 1)$ , 故 $W^{jk} = W^{r}$ 只有 N 个不同的值 $W^{0}, W^{1}, \cdots, W^{N-1}$ , 特别当 $N = 2^{m}$ 时, 只有 $\frac{N}{2}$ 个不同值, 因此可把同一个 $W^{r}$ 对应的 $B_{k}$ 相加后再乘 $W^{r}$ , 这就能大量减少乘法的次数.

下面仅介绍 $N = 2^{m}$ 时的算法. 把 $k, j$ 分别用二进制表示为

其中 $k_{r}, j_{r} (r = 0, 1, \cdots, m - 1)$ 只能取 0 或 1, 相应地令

于是 (4.55) 可分解为 $m$ 层求和, 即

上式的 $m$ 层求和由里往外, 分别引入记号

由此看到

为简化每个和式的计算, 利用 $W^{j_{0}2^{m-1}} = W^{j_{0}\frac{N}{2}} = (-1)^{j_{0}}$ , 并将二进制 $(0k_{m-2}\cdots k_{0})_{2} = k$ 表示为 $k = k_{m-2}2^{m-2} + \cdots + k_{0}2^{0}$ , 即为十进制数, 于是

由于 $j_{0}=0$ 或 1, 并将上式中的二进制数表示为十进制数可得, 对于 $k=0,1,\cdots,2^{m-1}-1$ 有

同理可推出, 对于 j=0,1 以及 $k=0,1,\cdots,2^{m-2}-1$ 有

一般情况可得, 对于 $l=1,2,\cdots,m, j=0,1,\cdots,2^{l-1}-1$ 以及 $k=0,1,\cdots,2^{m-l}-1$ , 有

公式 (4.59) 就是计算 DFT 的快速算法公式, 它比原始的 DFT 算法 (相应于公式 (4.53)) 有所改进, 方法的详细推导见 [14]. 此方法计算全部 $\{C_j\}_0^{N-1}$ 共用 $\frac{N}{2}(m-1)$ 次复数乘法和 $Nm$ 次加法运算, 当 $N=2^{10}$ 时, 它是非快速算法运算量的 $\frac{1}{230}$ .

在 Matlab 软件中, FFT(f) 是一个内部函数, 其值是向量 f 的离散傅里叶系数的 N 倍, 在数字信号处理、全息技术、光谱和声谱分析等领域称该值为向量 f 的离散谱 C. IFFT(C) 也是一个内部函数, 它是 DFT 的逆变换.

例 4.5.1 给出一张记录 $\{f_{j}\}=(2,1,0,1)$ ，用公式 (4.55) 和 Matlab 软件计算 $\{f_{j}\}$ 的离散谱 $\{C_{j}\}$ .

解：此时公式 (4.55) 中的 $N = 4, W = \mathrm{e}^{-i\frac{2\pi}{4}} = -i$ ，故

用 Matlab 计算该题的命令为

$$

f = [ 2101 ] $$

$$

\mathrm{C} = \mathrm{fft(f)} $$

评 注#

数值逼近是数值计算中的基本问题, 有着广泛的应用背景. 本章除介绍最佳一致逼近和最佳平方逼近的基本思想及算法外, 还介绍了非线性最小二乘问题和矛盾方程组的求解. 我们认为对此问题的介绍有助于更好地理解线性最小二乘问题. 另外, 考虑到傅里叶变换在工程中的广泛应用, 本章对周期函数的最佳平方逼近多项式与周期函数的傅里叶级数之间的关系、周期的离散数据的最小二乘解、快速傅里叶变换的算法思想作了简要的介绍.

习题四#

1. 求下列函数在区间 $[-1,1]$ 上的线性最佳一致逼近多项式： (1) $f(x)=x^{2}+3x-5;$ (2) $f(x)=\mathrm{e}^{x}.$

2. 令 $M$ 和 $m$ 分别代表连续函数 $f(x)$ 在区间 $[a, b]$ 上的最大值和最小值. 证明 $f(x)$ 在区间 $[a, b]$ 上的零次最佳一致逼近多项式是 $p(x) = \frac{1}{2}(M + m)$ .

3. 试分别求函数 $f(x) = \sqrt{1 + x}$ 在区间[0,1]上的一次最佳一致逼近多项式和一次最佳平方逼近多项式.

4. 求函数 $f(x) = \sin (\pi x)$ 在区间[0,1]上的二次最佳平方逼近多项式.

5. 用最小二乘法, 找出形如 $y = ax^2 + b$ 的抛物线方程, 使之最佳地代表表4-7中的数据.

表 4-7 习题 5 中给定的拟合数据

6. 证明如果用最小二乘法使一条直线拟合数据表 $(x_{i},y_{i})$ ，那么这条直线必通过点 $(x^{*},y^{*})$ 这里 $x^{*}$ 和 $y^{*}$ 分别是 $x_{i}$ 和 $y_{i}$ 的平均值.

7. 已知液体的粘性度 $V$ 随温度按照 $V = a + bT + cT^2$ 变化, 利用表4-8的数值, 找出 $a, b$ 和 $c$ 的最佳值.

表 4-8 习题 7 中给定的粘度-温度数据

8. 找出形如 $a \sin \pi x + b \cos \pi x$ 的函数, 使之在最小二乘意义下拟合表 4-9 中的数据点:

表 4-9 习题 8 中给定的拟合数据

数值实验四#

1. 拟合形如 $f(x) \approx \frac{a + bx}{1 + cx}$ 的函数的一种快速方法是将最小二乘法用于下列问题： $f(x)(1 + cx) \approx a + bx$ ，试用这一方法拟合表4-10给出的中国人口数据.

表 4-10 中国人口数据

2. 已知液体的表面张力 $s$ 是温度 $T$ 的线性函数 $s = aT + b$ , 对某种液体有表 4-11 的实验数据. 试用最小二乘法确定系数 $a, b$ .

表 4-11 表面张力-温度数据

3. 给出一张记录 $\{f_{k}\}=(4,3,2,1,0,1,2,3)$ ，用 FFT 算法求出 $\{f_{k}\}$ 离散谱 $\{C_{k}\}$ .

第5章 数值积分与数值微分#