数值计算方法笔记

一、误差

1. 误差的产生

从纯数学的角度，我们求解问题的答案是完全精确的但在实际工程问题中，我们是很难真正得到所谓的精确解，常常只需要近似解就可以了。

误差的来源一般有：

模型误差：从实际问题中抽象出的数学模型，一般只研究主要因素，所以往往是实际问题的近似描述；（由于这类误差难以定量分析，所以本课程假定研究的数学模型是合理的）
观测误差：通过测量得到模型中参数的值，由于测量工具和测量手段的限制，也会有误差；（本课程也不对此类误差作过多研究）
截断误差（方法误差）：求近似解时产生的误差；
舍入误差：由于机器字长有限导致的误差；

例：近似计算 $\int_{0}^{1} e^{- x^{2}} d x$ .

解：原函数不存在，考虑将 $e^{- x^{2}}$ 作Taylor展开后再积分， $\int_{0}^{1} e^{- x^{2}} d x = \int_{0}^{1} (1 - x^{2} + \frac{x ^{4}}{2 !} - \frac{x ^{6}}{3 !} + \frac{x ^{8}}{4 !} - \dots) d x = 1 - \frac{1}{3} + \frac{1}{2 !} \times \frac{1}{5} - \frac{1}{3 !} \times \frac{1}{7} + \frac{1}{4 !} \times \frac{1}{9} - \dots$ ，若将前四项记作 $S_{4}$ ，取 $\int_{0}^{1} e^{- x^{2}} d x \approx S_{4} = 1 - \frac{1}{3} + \frac{1}{10} - \frac{1}{42} \approx 1 - 0.333 + 0.1 - 0.024 = 0.743$ ，则计算分数时的四舍五入导致了舍入误差，此时 $∣ 舍入误差 ∣ < 0.0005 \times 2 = 0.001$ ，而 $R_{4} = \frac{1}{4 !} \times \frac{1}{9} - \frac{1}{5 !} \times \frac{1}{11} + \dots$ 称为截断误差，此时 $∣ R_{4} ∣ < \frac{1}{4 !} \times \frac{1}{9} < 0.005$ ，因此 $∣ 计算 \int_{0}^{1} e^{- x^{2}} d x 的总体误差 ∣ < 0.005 + 0.001 = 0.006$ .

2. 误差的有效数字

定义：我们把 $e (x) = x - x^{*}$ 称为绝对误差，其中 $x$ 为精确值， $x^{*}$ 为 $x$ 的近似值。由于精确值往往无法得知，我们记 $∣ e (x) ∣$ 的上限为 $ε (x)$ ，并称其为绝对误差限，工程上常记为 $x = x^{*} \pm ε$ ，例如 $x = 200 \pm 0.25$ 。

$e (x)$ 理论上是唯一确定的，可能取正也可能取负。但 $ε (x) > 0$ 并不唯一，当然 $ε (x)$ 越小越具有参考价值。

定义：我们把 $e_{r} (x) = \frac{e ( x )}{x}$ 称为相对误差， $x$ 的相对误差上限定义为 $ε_{r} (x) = \frac{ε ( x )}{∣ x ^{*} ∣}$ 。

有时会把 $ε (x)$ 简写为 $ε$ ，把 $ε_{r} (x)$ 简写为 $ε_{r}$ 。

定义：若近似值 $x^{*}$ 的绝对误差限是某一位的半个单位，就称 $x^{*}$ 精确到这一位，且若从该位直到 $x^{*}$ 的第一位非零数字共有 $n$ 位，则称近似值 $x^{*}$ 有 $n$ 位有效数字。即，假设 $x^{*} = a_{1} a_{2} \dots a_{m} . a_{m + 1} \dots a_{n}$ （其中 $a_{1} \neq = 0$ ），则称 $x^{*}$ 有 $n$ 位有效数字，可用科学计数法表示为 $x^{*} = \pm 0. a_{1} a_{2} \dots a_{n} \times 1 0^{m}$ ，则 $∣ x - x^{*} ∣ \leq 0.5 \times 1 0^{m - n}$ 。

例： $π = 3.1415926535897932 \dots$ ，若取 $π^{*} = 3.1415$ ，问 $π^{*}$ 有几位有效数字？请证明。

解： $∵ π^{*} = 0.31415 \times 1 0^{1}$ ，且 $∣ π^{*} - π ∣ = 0.0000926 \dots < 0.5 \times 1 0^{- 3} = 0.5 \times 1 0^{1 - 4}$ ， $∴ π^{*}$ 有4位有效数字，精确到小数点后3位。

注意： $0.2300$ 有4位有效数字，而 $00023$ 只有2位有效数字。12300如果写成 $0.123 \times 1 0^{5}$ ，则表示只有3位有效数字。（数字末尾的0不可随意省去）

有效数字与相对误差的关系1：如果 $x^{*}$ 有 $n$ 位有效数字，那么它的相对误差限为 $ε_{r} \leq \frac{1}{2 a _{1}} \times 1 0^{- n + 1}$ 。

证明：根据定义， $ε_{r} = ∣ \frac{ε}{x ^{*}} ∣ \leq \frac{0.5 \times 1 0 ^{m - n}}{0. a _{1} a _{2} \dots a _{n} \times 1 0 ^{m}} = \frac{1 0 ^{- n}}{2 \times 0. a _{1} \dots} \leq \frac{1}{2 a _{1}} \times 1 0^{- n + 1}$ .

有效数字与相对误差的关系2：已知 $x^{*}$ 的相对误差限可写为 $ε_{r} = \frac{1}{2 ( a _{1} + 1 )} \times 1 0^{- n + 1}$ ，则 $x^{*}$ 至少有 $n$ 位有效数字。

证明： $∣ x - x^{*} ∣ = ε_{r} \cdot ∣ x^{*} ∣ = \frac{1 0 ^{- n + 1}}{2 ( a _{1} + 1 )} \times 0. a_{1} a_{2} \dots \times 1 0^{m} < \frac{1 0 ^{- n + 1}}{2 ( a _{1} + 1 )} \cdot (a_{1} + 1) \times 1 0^{m - 1} = 0.5 \times 1 0^{m - n}$ ，因此，根据有效数字的定义， $x^{*}$ 至少有 $n$ 位有效数字。

例：为使 $π^{*}$ 的相对误差小于 $0.001%$ ，至少应取几位有效数字？

解：假设 $π^{*}$ 取到 $n$ 为有效数字，则其相对误差上限为 $ε_{r} \leq \frac{1}{2 a _{1}} \times 1 0^{- n + 1}$ ，要保证其相对误差小于 $0.001%$ ，只要保证其上限满足 $ε_{r} \leq \frac{1}{2 a _{1}} \times 1 0^{- n + 1} < 0.001%$ ，已知 $a_{1} = 3$ ，则从以上不等式可解得 $n > 6 - lo g 6$ ，即 $n \geq 6$ ，应取 $π^{*} = 3.14159$ 。

3. 误差的积累与传播

误差的积累

例：计算 $I_{n} = \int_{0}^{1} x^{n} e^{x - 1} d x$ ， $n = 0, 1, 2, \dots$ .

解：估算 $I_{n}$ 的范围： $∵ e^{- 1} < e^{x - 1} < e^{1 - 1}$ （ $x \in (0, 1)$ ）， $∴ \int_{0}^{1} x^{n} \cdot e^{- 1} d x < I_{n} < \int_{0}^{1} x^{n} \cdot e^{1 - 1} d x$ ，即 $\frac{1}{e ( n + 1 )} < I_{n} < \frac{1}{n + 1}$ 。

公式一：用分部积分法得 $I_{n} = \int_{0}^{1} x^{n} e^{x - 1} d x = x^{n} e^{x - 1} ∣_{0}^{1} - \int_{0}^{1} n x^{n - 1} e^{x - 1} d x = 1 - n I_{n - 1}$ ， $I_{0} = \int_{0}^{1} e^{x - 1} d x = 1 - \frac{1}{e} \approx 0.63212056 = 记为 I_{0}^{*}$ ，则初始误差 $∣ E_{0} ∣ = ∣ I_{0} - I_{0}^{*} ∣ < 0.5 \times 1 0^{- 8}$ ，依次计算 $I_{1}^{*} \dots I_{15}^{*}$ ，得 $I_{1}^{*} = 1 - 1 \cdot I_{0}^{*} = 0.36787944 \dots I_{10}^{*} = 1 - 10 \cdot I_{9}^{*} = 0.08812800 I_{11}^{*} = 1 - 11 \cdot I_{10}^{*} = 0.03059200 I_{12}^{*} = 1 - 12 \cdot I_{11}^{*} = 0.63289600 I_{13}^{*} = 1 - 13 \cdot I_{12}^{*} = - 7.2276480 I_{14}^{*} = 1 - 14 \cdot I_{13}^{*} = 94.959424 I_{15}^{*} = 1 - 15 \cdot I_{14}^{*} = - 1423.3814$ ，可见从 $I_{12}^{*}$ 开始就已经不在之前估算出来的积分范围内了，考察第 $n$ 步的误差 $E_{n}$ ， $∣ E_{n} ∣ = ∣ I_{n} - I_{n}^{*} ∣ = ∣ (1 - n I_{n - 1}) - (1 - n I_{n - 1}^{*}) ∣ = n ∣ E_{n - 1} ∣ = \dots = n! ∣ E_{0} ∣$ ，可见初始的小扰动 $∣ E_{0} ∣ < 0.5 \times 1 0^{- 8}$ 迅速积累，误差呈递增走势。

公式二： $I_{n} = 1 - n I_{n - 1} \Rightarrow I_{n - 1} = \frac{1}{n} (1 - I_{n})$ ，此公式理论上与公式一等价。要使用此公式，首先要估算一个 $I_{m}$ ，再反推要求的 $I_{n}$ （ $n ≪ m$ ）， $∵ \frac{1}{e ( m + 1 )} < I_{m} < \frac{1}{m + 1}$ ，可取区间的中点作为 $I_{m}$ 的估计值，即 $I_{m}^{*} = \frac{1}{2} [\frac{1}{e ( m + 1 )} + \frac{1}{m + 1}] \approx I_{m}$ ，估算出 $I_{15}^{*} = \frac{1}{2} [\frac{1}{e \cdot 16} + \frac{1}{16}] \approx 0.042746233$ ，则 $I_{14}^{*} = \frac{1}{15} (1 - I_{15}^{*}) \approx 0.063816918 I_{13}^{*} = \frac{1}{14} (1 - I_{14}^{*}) \approx 0.066870220 I_{12}^{*} = \frac{1}{13} (1 - I_{13}^{*}) \approx 0.0071779214 \dots I_{1}^{*} = \frac{1}{2} (1 - I_{2}^{*}) \approx 0.36787944 I_{0}^{*} = \frac{1}{1} (1 - I_{1}^{*}) \approx 0.063212056$ ，考察反推一步的误差 $∣ E_{m - 1} ∣ = ∣ \frac{1}{m} (1 - I_{m}) - \frac{1}{m} (1 - I_{m}^{*}) ∣ = \frac{1}{m} ∣ E_{m} ∣$ ，依此类推，对 $n < m$ 有 $∣ E_{n} ∣ = \frac{1}{m ( m - 1 ) \dots ( n + 1 )} ∣ E_{m} ∣$ ，误差逐步递减，这样的算法称为稳定的算法。

四则运算的误差传播

两个近似数进行加、减、乘、除运算后的绝对误差：

$ε (x_{1} \pm x_{2}) \approx ε (x_{1}) + ε (x_{2})$ ；
$ε (x_{1} x_{2}) \approx ∣ x_{2} ∣ ε (x_{1}) + ∣ x_{1} ∣ ε (x_{2})$ ；
$ε (\frac{x _{1}}{x _{2}}) \approx \frac{ε ( x _{1} x _{2} )}{∣ x _{2} ∣ ^{2}} \approx \frac{∣ x _{2} ∣ ε ( x _{1} ) + ∣ x _{1} ∣ ε ( x _{2} )}{∣ x _{2} ∣ ^{2}}$ ，（ $x_{2} \neq = 0$ ）。

两个近似数进行加、减、乘、除运算后的相对误差：

$ε_{r} (x_{1} \pm x_{2}) \approx \frac{ε ( x _{1} \pm x _{2} )}{∣ x _{1} \pm x _{2} ∣} \approx \frac{ε ( x _{1} ) + ε ( x _{2} )}{∣ x _{1} \pm x _{2} ∣} \approx \frac{x _{1} ε _{r} ( x _{1} ) + x _{2} ε _{r} ( x _{2} )}{∣ x _{1} \pm x _{2} ∣}$ ；
$ε_{r} (x_{1} x_{2}) \approx ε_{r} (x_{1}) + ε_{r} (x_{2})$ ；
$ε_{r} (\frac{x _{1}}{x _{2}}) \approx ε_{r} (x_{1}) + ε_{r} (x_{2})$ 。

函数的误差传播

对于 $y = f (x)$ ，若 $x$ 存在误差，对 $y$ 产生的影响：

绝对误差： $∣ e (y) ∣ \approx ∣ f^{'} (x^{*}) ∣ \cdot ∣ e (x) ∣$ ；

推导： $e (y) = f (x^{*}) - f (x) = 中值定理 f^{'} (ξ) (x^{*} - x)$ ，由于 $x^{*}$ 与 $x$ 非常接近，可认为 $f^{'} (ξ) \approx f^{'} (x^{*})$ ，再由于 $e (x) = x^{*} - x$ ，故 $∣ e (y) ∣ \approx ∣ f^{'} (x^{*}) ∣ \cdot ∣ e (x) ∣$ ，即 $x^{*}$ 产生的误差经过 $f$ 作用后被放大或缩小了 $∣ f^{'} (x^{*}) ∣$ 倍，故称 $∣ f^{'} (x^{*}) ∣$ 为放大因子或绝对条件数。
相对误差： $∣ e_{r} (y) ∣ \approx ∣ \frac{x ^{*} \cdot f ^{'} ( x ^{*} )}{f ( x ^{*} )} ∣ \cdot ∣ e_{r} (x) ∣$ 。

推导：由于 $∣ e_{r} (x) ∣ = ∣ \frac{e ( x )}{x ^{*}} ∣$ ，故 $∣ e_{r} (y) ∣ = ∣ \frac{e ( y )}{f ( x ^{*} )} ∣ = ∣ \frac{f ( x ^{*} ) - f ( x )}{x ^{*} - x} \cdot \frac{x ^{*}}{f ( x ^{*} )} \cdot \frac{x ^{*} - x}{x ^{*}} ∣ \approx ∣ \frac{x ^{*} \cdot f ^{'} ( x ^{*} )}{f ( x ^{*} )} ∣ \cdot ∣ e_{r} (x) ∣$ ，我们称 $∣ \frac{x ^{*} \cdot f ^{'} ( x ^{*} )}{f ( x ^{*} )} ∣$ 为相对误差条件数。

$f$ 的条件数在 $x$ 处小，则称 $f$ 在该点是好条件（ $f$ 的条件数在 $x$ 处大，则称 $f$ 在该点是坏条件）。

4. 误差的控制（算法的稳定性）

为把误差控制在一定范围内，构造算法时要遵循以下准则：

避免大数“吃”小数。避免方法：求和时从小到大相加，可使和的误差减小；

因为如果有的数相差数量级，而计算机字长有限，有时小数会被“吃掉”。

例如：用单精度的计算器按从小到大、以及从大到小的顺序分别计算 $1 + 2 + 3 + \dots + 40 + 1 0^{9}$ ，

如果按从小到大算，先计算 $1 + \dots + 40 = 820$ ，由于 $820$ 与 $1 0^{9}$ 相差数量级，先将 $1 0^{9}$ 写成 $1 0^{9} = 0.1 \times 1 0^{10}$ ，再把 $820$ 写成 $820 = 0.0000 000820 \times 1 0^{10}$ ，由于单精度只能保留八位有效数字， $20$ 会被四舍五入舍弃，因此 $820 + 1 0^{9} \approx 0.1000 0008 \times 1 0^{10}$ ；

如果从大到小的算，先计算 $40 + 1 0^{9} = 0.0000 00004 \times 1 0^{10} + 0.1 \times 1 0^{10} \approx 0.1000 0000 \times 1 0^{10}$ ，因此 $40$ 会被 $1 0^{9}$ 吃掉，其他数同理，最后的计算结果为 $1 + 2 + 3 + \dots + 40 + 1 0^{9} \approx 1 0^{9}$ ，只有一位有效数字。
避免两个相近的数相减。避免方法：用等价的公式替换，如 $x + ε - x = \frac{ε}{x + ε + x}$ 、 $ln (x + ε) - ln x = ln (\frac{x + ε}{x})$ 、当 $∣ x ∣ ≪ 1$ 时 $1 - cos x = 2 sin^{2} \frac{x}{2}$ 、当 $∣ x ∣ ≪ 1$ 时 $e^{x} - 1 = x (1 + \frac{1}{2} x + \frac{1}{6} x^{2} + \dots)$ ；

两个相近的数相减会导致有效数字的严重损失，导致误差增大。

例如： $a_{1} = 0.12345$ ， $a_{2} = 0.12346$ ，各有 $5$ 位有效数字，而 $a_{2} - a_{1} = 0.00001$ ，只剩下一位有效数字，如果精度不够，很容易丢失。
避免小分母。避免方法：用等价的公式替换，如当 $x$ 接近零时将 $\frac{1 - c o s x}{s i n x}$ 转化为等值式 $\frac{s i n x}{1 + c o s x}$ ；

分母小会造成浮点溢出。
先化简再计算，减少步骤，避免误差积累；

例如：直接计算 $x^{127}$ 需要做 $126$ 次乘法，但若写成 $x^{127} = x \cdot x^{2} \cdot x^{4} \cdot x^{8} \cdot x^{16} \cdot x^{32} \cdot x^{64}$ ，只需要12次乘法运算。

再如：计算多项式 $P_{n} (x) = a_{0} + a_{1} x + a_{2} x^{2} + \dots + a_{n} x^{n}$ ，若直接计算需 $\frac{n ( n + 1 )}{2}$ 次乘法和 $n$ 次加法，但若改写为 $P_{n} (x) = a_{0} + x (a_{1} + x (a_{2} + \dots + x (a_{n - 1} + x a_{n})) \dots)$ ，则只需做 $n$ 次乘法和 $n$ 次加法。
选用稳定的算法。

例如前面误差的积累中的题目，如果用 $I_{n} = 1 - n I_{n - 1}$ 就不是好算法，但若反过来，用 $I_{n - 1} = \frac{1}{n} (1 - I_{n})$ 计算，则误差就会越来越小，所以这是一个好算法，是稳定的。

二、非线性方程的数值解法

1. 方程求根的基本思想

对于方程 $f (x) = 0$ ，当 $f (x)$ 不是 $x$ 的线性函数时，称对应的函数方程为非线性方程，如果 $f (x)$ 是多项式函数，则称为代数方程，否则称为超越方程（如三角方程，指数、对数方程等）。

非线性方程根的定义：非线性方程 $f (x) = 0$ 的根，亦称为函数 $f (x)$ 的零点。如果 $f (x)$ 可以分解为 $f (x) = (x - x^{*})^{m} g (x)$ ，其中 $m$ 为正整数且 $g (x^{*}) \neq = 0$ ，则称 $x^{*}$ 是 $f (x)$ 的m重零点，或称方程 $f (x) = 0$ 的m重根，当 $m = 1$ 时称 $x^{*}$ 为单根。若 $f (x)$ 存在 $m$ 阶导数，则是方程 $f (x)$ 的 $m$ 重根（ $m > 1$ ），当且仅当 $f (x^{*}) = f^{'} (x^{*}) = \dots = f^{(m - 1)} (x^{*}) = 0$ ， $f^{(m)} (x^{*}) \neq = 0$ 。

求根的算法思路：

（1）根的存在性：方程有没有根？如果有，有几个根？

（2）这些根大致在哪里？如何把根隔离开来？

（3）根的精确化。

根的存在性定理（零点定理）：设函数 $f (x)$ 在区间 $[a, b]$ 上连续，如果 $f (a) \cdot f (b) < 0$ 。则方程 $f (x) = 0$ 在 $[a, b]$ 内至少有一实根 $x^{*}$ 。

根据该定理找根的方法一：画出 $f (x)$ 的略图，从而看出曲线与 $x$ 轴交点的位置。方法二：从左端点 $x = a$ 出发，按某个预先选定的步长 $h$ 一步一步地向右跨，每跨一步都检验每步起点 $x_{0}$ 与终点 $x_{0} + h$ 的函数值，若 $f (x_{0}) \cdot f (x_{0} + h) \leq 0$ ，那么所求的根 $x^{*}$ 必在 $x_{0}$ 与 $x_{0} + h$ 之间，这里可取 $x_{0}$ 或 $x_{0} + h$ 作为根的初始近似。

2. 非线性方程求根-二分法

对于方程 $f (x) = 0$ ，设函数 $f (x)$ 在闭区间 $[a, b]$ 上连续，且 $f (a) f (b) < 0$ ，根据连续函数的性质可知， $f (x) = 00$ 在 $(a, b)$ 内必有实根，称区间 $[a, b]$ 为有根区间。这里为了方便讨论，假定方程 $f (x) = 0$ 在区间 $[a, b]$ 内有唯一实根 $x^{*}$ 。

二分法的思想：先确定有根区间，将区间二等分，根据零值定理，通过判断 $f (x)$ 的符号，逐步将有根区间缩小，直至有根区间足够小，便可求出满足精度要求的近似根。

圈定根所在的范围的过程叫做圈定根或根的隔离，圈定好根后，采取适当的数值方法确定有一定精度要求的初值。

对于代数方程，其根的个数（实根或复根）与其次数相同；对于超越方程，其根可能是一个、几个或无解，并没有固定的圈根方法。

求方程根的问题，从几何上讲，就是求曲线 $y = f (x)$ 与 $x$ 轴交点的横坐标，方法有画图法、逐步搜索法等：

画图法：比如 $x lo g x - 1 = 0$ 可改写为 $lo g x = \frac{1}{x}$ ，分别画出 $lo g x$ 和 $\frac{1}{x}$ 的图像可知它们交点的横坐标在 $[2, 3]$ 内；

逐步搜索法：对于给定的 $f (x)$ ，设有根区间为 $[A, B]$ ，从 $x_{0} = A$ 出发，以步长 $h = \frac{( B - A )}{n}$ （ $n$ 是正整数），在 $[A, B]$ 内取定节点 $x_{i} = x_{0} + ih$ （ $i = 0, 1, 2, \dots n$ ），从左至右检查 $f (x_{i})$ 的符号，如发现 $x_{i}$ 与端点 $x_{0}$ 的函数值异号，则得到一个缩小的有根子区间 $[x_{i - 1}, x_{i}]$ 。

例：确定方程 $f (x) = x^{3} - x - 1 = 0$ 的有根区间.

解：不难发现 $f (0) < 0$ ， $f (2) > 0$ ，故在区间 $(0, 2)$ 内至少有一个实根，设从 $x = 0$ 出发，取 $h = 0.5$ 为步长，进行根的搜索，列表如下：

$x$ $0$ $0.5$ $1.0$ $1.5$ $2$

$f (x)$ $-$ $-$ $-$ $+$ $+$

可以看出，在 $[1.0, 1.5]$ 内必有一根。

$x$	$0$	$0.5$	$1.0$	$1.5$	$2$
$f (x)$	$-$	$-$	$-$	$+$	$+$

二分法是对逐步搜索法的改进，二分法的具体步骤如下：

（1）取有根区间 $[a, b]$ 的中点，将它分为两半，分点为 $x_{0} = \frac{a + b}{2}$ ，判断 $a$ 、 $b$ 、 $x_{0}$ 的符号，根据零点定理确定新的有根区间 $[a_{1}, b_{1}]$ ，这样就缩小了有根区间。

（2）对压缩了的有根区间 $[a_{1}, b_{1}]$ ，施行同样的手法，即再取中点 $x_{1} = \frac{a _{1} + b _{1}}{2}$ ，将区间 $[a_{1}, b_{1}]$ 再分为两半，然后再确定有根区间 $[a_{2}, b_{2}]$ ，其长度是 $[a_{1}, b_{1}]$ 的二分之一。

（3）依此类推，若不出现 $f (x_{k}) = 0$ ，即可得到一系列有根区间序列 $[a, b] \supset [a_{1}, b_{1}] \supset [a_{2}, b_{2}] \supset \dots \supset [a_{k}, b_{k}]$ ，上述的每个区间都是前一个区间的一半，因此 $[a_{k}, b_{k}]$ 的长度 $b_{k} - a_{k} = \frac{1}{2} (b_{k - 1} - a_{k - 1}) = \dots = \frac{1}{2 ^{k}} (b - a)$ ，当 $k \to \infty$ 时趋于零，这些区间最终会收敛于一点 $x^{*}$ ，即为所求的根。

（4）终止条件：每次二分后，将 $[a_{k}, b_{k}]$ 的中点 $x_{k} = \frac{1}{2} (a_{k} + b_{k})$ 作为根的近似值，可以得到一个近似根的序列 $x_{0}, x_{1}, \dots, x_{k}, \dots$ ，该序列以根 $x^{*}$ 为极限，只要二分次数足够多（即 $k$ 足够大），便有 $∣ x^{*} - x_{k} ∣ < ε$ （ $ε$ 为给定的精度），由于 $x^{*} \in [a_{k}, b_{k}]$ ，故 $∣ x^{*} - x_{k} ∣ \leq \frac{b _{k} - a _{k}}{2} = \frac{b - a}{2 ^{k + 1}}$ ，若给定的精度 $ε > 0$ ，要使 $∣ x^{*} - x_{k} ∣ < ε$ 成立，只要取 $k$ 满足 $\frac{1}{2 ^{k + 1}} (b - a) < ϵ$ 即可，也就是说当 $k \geq \frac{l g ( b - a ) - l g ε}{l g 2} - 1$ 时，做到第 $k + 1$ 次二分，计算得到的 $x_{k}$ 就是满足精度要求的近似根。（程序中常用相邻的 $x_{k}$ 与 $x_{k - 1}$ 的差的绝对值或 $a_{k}$ 与 $b_{k}$ 的差的绝对值是否小于 $ε$ 来决定二分的次数）

例：证明方程 $f (x) = x^{3} - 2 x - 5 = 0$ 在区间 $[2, 3]$ 内有一个根，使用二分法求误差不超过 $0.5 \times 1 0^{- 3}$ 的根要二分几次？

解：因为 $f (2) = - 1 < 0$ 、 $f (3) = 16 > 0$ ，且 $f (x)$ 在 $[2, 3]$ 上连续，故 $f (x) = 0$ 在 $[2, 3]$ 内至少有一个根；又 $f^{'} (x) = 3 x^{2} - 2$ 当 $x \in [2, 3]$ 时， $f^{'} (x) > 0$ ，故 $f (x)$ 在 $[2, 3]$ 内单调递增，从而 $f (x)$ 在 $[2, 3]$ 内有且仅有一根。

给定误差 $ε = 0.5 \times 1 0^{- 3}$ ，使用二分法时，只要取 $k$ 满足 $\frac{1}{2 ^{k + 1}} (b - 1) \leq \frac{1}{2} \times 1 0^{- 3} \Rightarrow 2^{k} \geq 1 0^{3} \Rightarrow k \geq \frac{3 l g 10}{l g 2} = 9.97$ ，所以需二分10次便可达到要求。

二分法的优点是不管有根区间 $[a, b]$ 多大，总能求出满足精度要求的根，且对 $f (x)$ 要求不高，只要连续即可，计算也简单。它的局限性是只能用于求函数的实根，不能用于求复根及重根，它的收敛速度与比值为 $\frac{1}{2}$ 的等比级数相同。

3. 非线性方程求根-迭代法

迭代法是一种逐次逼近的方法，用某个固定公式反复校正根的近似值，使之逐步精确化，最后得到满足精度要求的结果。

迭代法的思想：为求解非线性方程 $f (x) = 0$ 的根，先将其写成便于迭代的等价方程 $x = φ (x)$ ，其中 $φ (x)$ 为 $x$ 的连续函数，即，如果数 $x^{*}$ 使 $f (x^{*}) = 0$ ，则也有 $x^{*} = φ (x^{*})$ ，反之，若 $x^{*} = φ (x^{*})$ ，则也有 $f (x^{*}) = 0$ ，那么称 $φ (x)$ 为迭代函数。任取一个初值 $x_{0}$ ，代入式 $x = φ (x)$ 的右端，得到 $x_{1} = φ (x_{0})$ ，再将 $x_{1}$ 代入式 $x = φ (x)$ 的右端，得到 $x_{2} = φ (x_{1})$ ，依此类推，得到一个数列 $x_{1}, x_{2}, x_{3}, \dots$ ，其中 $x_{k + 1} = φ (x_{k})$ （ $k = 0, 1, 2, \dots$ ），此式称为求解非线性方程的简单迭代法，如果由迭代格式 $x_{k + 1} = φ (x_{k})$ 产生的序列 ${x_{k}}$ 收敛，即 $n \to \infty lim x_{n} = x^{*}$ ，则称迭代法收敛（实际计算中不可能也没必要无穷多步地算下去，对预先给定的精度 $ε$ ，只要满足 $∣ x_{k} - x_{k - 1} ∣ < ε$ 即可结束计算，并取 $x^{*} \approx x_{k}$ ）。

例：用迭代法求方程 $x^{3} - x - 1 = 0$ 在 $x = 1.5$ 附近的一个根.

解：将方程改写成如下两种等价形式： $x = φ_{1} (x) = 3 x + 1 x = φ_{2} (x) = x^{3} - 1$ ，相应地可得到两个迭代公式 $x_{k + 1} = φ_{1} (x_{k}) = 3 x_{k} + 1 x_{k + 1} = φ_{2} (x_{k}) = x_{k}^{3} - 1$ ，如果取初始值 $x_{0} = 1.5$ ，用上述的两个迭代公式分别迭代，计算得到：

（1）对于 $x_{0} = 1.5$ ， $x_{k + 1} = 3 x_{k} + 1$ （ $k = 0, 1, 2, \dots$ ），则 $x_{1} = 1.35721$ 、 $x_{2} = 1.33086$ 、 $x_{3} = 1.32588$ 、 $x_{4} = 1.32494$ 、 $x_{5} = 1.32476$ 、 $x_{6} = 1.32473$ 、 $x_{7} = 1.32472$ 、 $\dots$

（2）对于 $x_{0} = 1.5$ ， $x_{k + 1} = x_{k}^{3} - 1$ （ $k = 0, 1, 2, \dots$ ），则 $x_{1} = 2.375$ 、 $x_{2} = 12.39$ 、 $\dots$

可见第一种迭代方式是逐步收敛的，而第二种会随着迭代次数增加差异越来越大，于是可以知道，有的迭代公式可以求出根，有的公式则不能，因此判断迭代公式是否收敛尤为重要。

迭代法的几何意义：通常将方程 $f (x) = 0$ 化为与它同解的方程 $x = φ (x)$ 的方法不止一种，有的收敛，有的不收敛，这取决于 $φ (x)$ 的形态，类似画图法，方程 $x = φ (x)$ 的求根问题在几何上就是确定曲线 $y = φ (x)$ 与直线 $y = x$ 的交点 $P^{*}$ 的横坐标，如下图所示：

定理（迭代法的收敛性）：设函数在 $[a, b]$ 上有连续的一阶导数，且满足①对所有的 $x \in [a, b]$ 有 $φ (x) \in [a, b]$ ；②存在 $0 < L < 1$ ，使所有的 $x \in [a, b]$ 有 $∣ φ^{'} (x) ∣ \leq L$ ；则方程 $x = φ (x)$ 在 $[a, b]$ 上的解 $x^{*}$ 存在且唯一，对任意的 $x_{0} \in [a, b]$ ，迭代过程 $x_{k + 1} = φ (x_{k})$ 均收敛于 $x^{*}$ ，并有如下的误差估计公式：① $∣ x^{*} - x ∣ \leq \frac{L}{1 - L} ∣ x_{k} - x_{k - 1} ∣$ ；② $∣ x^{*} - x_{k} ∣ \leq \frac{L ^{k}}{1 - L} ∣ x_{1} - x_{0} ∣$ 。

例：对方程 $x^{5} - 4 x - 2 = 0$ ，构造收敛的迭代格式.

解：容易判断 $[1, 2]$ 是方程的有根区间，且在此区间内 $f^{'} (x) = 5 x^{4} - 4 > 0$ ，所以此方程在 $[1, 2]$ 内有且仅有一根，将原方程改写为以下两种等价形式：

（1）对于 $x = \frac{x ^{5} - 2}{4}$ ，即 $φ (x) = \frac{x ^{5} - 2}{4}$ ，对于 $x \in [1, 2]$ ， $∣ φ^{'} (x) ∣ = \frac{5 x ^{4}}{4} > 1$ ，不满足收敛条件；

（2）对于 $x = 5 4 x + 2$ ，即 $φ (x) = 5 4 x + 2$ ，对于 $x \in [1, 2]$ ， $∣ φ^{'} (x) ∣ = \frac{1}{5 5 ( 4 x + 2 ) ^{4}} < \frac{1}{5 5 ( 4 + 2 ) ^{4}} \approx 0.2 < 1$ ，此时迭代公式满足迭代收敛条件。

定理：设 $φ (x)$ 在 $x = φ (x)$ 的根 $x^{*}$ 的邻域中有连续的一阶导数，且 $∣ φ^{'} (x^{*}) ∣ < 1$ ，则迭代过程 $x_{k + 1} = φ (x_{k})$ 具有局部收敛性。

当迭代函数较复杂时，通常只能设法使迭代过程在根的邻域（局部）收敛。

例：已知方程 $x = φ (x)$ 在 $[a, b]$ 内有根 $x^{*}$ ，且在 $[a, b]$ 上满足 $∣ φ^{'} (x) - 3∣ < 1$ ，利用 $φ (x)$ 构造一个迭代函数 $g (x)$ ，使 $x_{k + 1} = g (x_{k})$ （ $k = 0, 1, 2, \dots$ ）局部收敛于 $x^{*}$ .

解：由 $x = φ (x)$ 得 $x - 3 x = φ (x) - 3 x$ ，令 $x = - \frac{1}{2} (φ (x) - 3 x) = g (x)$ ，则对于 $x \in [a, b]$ 有 $∣ g^{'} (x) ∣ = ∣ - \frac{1}{2} (φ^{'} (x) - 3) ∣ < \frac{1}{2} < 1$ ，故 $∣ g^{'} (x^{*}) ∣ < 1$ ，迭代公式 $x_{k + 1} = g (x_{k}) = - \frac{1}{2} (φ (x_{k}) - 3 x_{k})$ 局部收敛。

定义（收敛速度）：设迭代过程 $x_{k + 1} = φ (x_{k})$ 收敛于 $x = φ (x_{k})$ 的根 $x^{*}$ ，记迭代误差 $e_{k} = x^{*} - x_{k}$ ，若存在常数 $p$ （ $p \geq 1$ ）和 $c$ （ $c \neq = 0$ ）使得 $k \to \infty lim \frac{∣ e _{k + 1} ∣}{∣ e _{k} ∣ ^{p}} = c$ ，则称序列 ${x_{k}}$ 是 $p$ 阶收敛的， $c$ 称渐进误差常数，特别地， $p = 1$ 时称为线性收敛， $p = 2$ 时称为平方收敛， $1 < p < 2$ 时称为超线性收敛。

$p$ 的大小反映了迭代法收敛速度的快慢， $p$ 越大则收敛速度越快，故迭代法的收敛阶是对迭代法收敛速度的一种度量。

定理：设迭代过程 $x_{k + 1} = φ (x_{k})$ ，若 $φ^{(p)} (x)$ 在所求根 $x^{*}$ 的邻域连续且 $φ^{'} (x^{*}) = φ^{''} (x) = \dots = φ^{(p - 1)} (x^{*}) = 0$ ， $φ^{(p)} (x^{*}) \neq = 0$ ，则迭代过程在 $x^{*}$ 邻域是 $p$ 阶收敛的。

4. 牛顿迭代法

牛顿迭代法的思想：方程求根的另一种思路是用近似方程代替原方程去求根，牛顿迭代法就是将非线性方程线性化，用线性方程的解逼近非线性方程的解。取 $x_{0} \approx x^{*}$ ，将 $f (x)$ 在 $x_{0}$ 作一阶Taylor展开，得 $f (x) = f (x_{0}) + f^{'} (x_{0}) (x - x_{0}) + \frac{f ^{''} ( ξ )}{2 !} (x - x_{0})^{2}$ ，其中 $ξ$ 在 $x_{0}$ 和 $x$ 之间，令 $0 = f (x^{*}) = f (x_{0}) + f^{'} (x_{0}) (x^{*} - x_{0}) + \frac{f ^{''} ( ξ )}{2} (x^{*} - x_{0})^{2}$ ，如果将 $(x^{*} - x_{0})^{2}$ 看成高阶小量，则有 $0 = f (x^{*}) \approx f (x_{0}) + f^{'} (x_{0}) (x^{*} - x_{0})$ ，这样就变成了求线性方程的解，解得 $x^{*} \approx x_{0} - \frac{f ( x _{0} )}{f ^{'} ( x _{0} )}$ 。

牛顿迭代法的几何意义：如下图，曲线是 $f (x)$ ，它与 $x$ 轴交点为 $x^{*}$ ，而 $x^{*} \approx x_{0} - \frac{f ( x _{0} )}{f ^{'} ( x _{0} )}$ 是曲线 $y = f (x)$ 在 $(x_{0}, f (x_{0}))$ 点处的切线与 $x$ 轴的交点的横坐标 $x_{1}$ （故牛顿迭代法又被称为切线法），依此类推可得到迭代公式 $x_{k + 1} = x_{k} - \frac{f ( x _{k} )}{f ^{'} ( x _{k} )}$ 。

牛顿迭代公式： $x_{k + 1} = x_{k} - \frac{f ( x _{k} )}{f ^{'} ( x _{k} )}$ ，只要 $f \in C^{1}$ （即 $f (x)$ 一阶连续），每一步迭代都有 $f^{'} (x_{k}) \neq = 0$ ，而且 $k \to \infty lim x_{k} = x^{*}$ ，则 $x^{*}$ 就是 $f$ 的根。

定理（收敛的充分条件）：设 $f \in C^{2} [a, b]$ ，若① $f (a) f (b) < 0$ ；②在整个 $[a, b]$ 上 $f^{''}$ 不变号且 $f^{'} (x) \neq = 0$ ；③选取 $x_{0} \in [a, b]$ 使得 $f (x_{0}) f^{''} (x_{0}) > 0$ ；则牛顿迭代法产生的序列 ${x_{k}}$ 收敛到 $f (x)$ 在 $[a, b]$ 上的唯一根。

条件一保证了根的存在性，条件二保证了根的唯一性（二阶导不变号则一阶导单调），条件三保证了产生的序列单调有界，保证收敛。

定理（局部收敛性）设 $f \in C^{2} [a, b]$ ，若 $x^{*}$ 为 $f (x)$ 在 $[a, b]$ 上的根，且 $f^{'} (x^{*}) \neq = 0$ ，则存在 $x^{*}$ 的邻域 $B_{δ} (x^{*}) : ∣ x - x^{*} ∣ < δ$ 使得任取初始值 $x_{0} \in B_{δ} (x^{*})$ ，牛顿迭代法产生的序列 ${x_{k}}$ 收敛到 $x^{*}$ ，且满足 $k \to \infty lim \frac{x ^{*} - x _{k + 1}}{( x ^{*} - x _{k} ) ^{2}} = - \frac{f ^{''} ( x ^{*} )}{2 f ^{'} ( x ^{*} )}$ 。

前面两个定理都说明了，牛顿迭代法的收敛性依赖于 $x_{0}$ 的选取， $x_{0}$ 越接近 $x^{*}$ 就越有可能收敛，如下图：

例：用牛顿迭代法求方程 $f (x) = x^{3} - 7.7 x^{2} + 19.2 x - 15.3$ 在 $x_{0} = 1$ 附近的根.

解： $x_{k + 1} = x_{k} - \frac{x _{k}^{3} - 7.7 x _{k}^{2} + 19.2 x _{k} - 15.3}{3 x _{k}^{2} - 15.4 x _{k} + 19.2}$ ，列表如下：

$k$ $x_{k}$ $f (x)$

0 1.00 -2.8

1 1.41176 -0.727071

2 1.62424 -0.145493

3 1.6923 -0;.0131682

4 1.69991 -0.0001515

5 1.7 0

$k$	$x_{k}$	$f (x)$
0	1.00	-2.8
1	1.41176	-0.727071
2	1.62424	-0.145493
3	1.6923	-0;.0131682
4	1.69991	-0.0001515
5	1.7	0

例：用牛顿迭代法计算 $2$ .

解：要求 $x = 2$ 相当于求 $f (x) = x^{2} - 2 = 0$ 的正根，由牛顿迭代公式得 $x_{k + 1} = x_{k} - \frac{x _{k}^{2} - 2}{2 x _{k}} = \frac{1}{2} (x_{k} + \frac{2}{x _{k}})$ ，其中 $k = 0, 1, 2, \dots$ ，取初始值 $x_{0} = 1$ ，代入公式进行迭代得 $x_{1} = 1.5$ 、 $x_{2} = 1.416666667$ 、 $x_{3} = 1.414215686$ 、 $x_{4} = 1.414213562$ 、 $x_{5} = 1.414213562$ ，故 $2 = 1.414213562$ 。

5. 牛顿迭代法的改进

重根-加速收敛法

引入：前面的两个收敛定理（充分性定理和局部收敛定理）都要求在 $x^{*}$ 附近 $f^{'} (x) \neq = 0$ ，那么若 $f^{'} (x^{*}) = 0$ ，牛顿法是否仍收敛？

$f^{'} (x) = 0$ 意味着会出现重根，设 $x^{*}$ 是 $f$ 的 $n$ 重根，则可以把 $f (x)$ 写成 $f (x) = (x - x^{*})^{n} q (x)$ 的形式，且 $q (x^{*}) \neq = 0$ ，因为牛顿迭代法实际上是一种特殊的不动点迭代， $g (x) = x - \frac{f ( x )}{f ^{'} ( x )}$ ，则 $∣ g^{'} (x^{*}) ∣ = ∣1 - \frac{f ^{'} ( x ^{*} ) ^{2} - f ( x ^{*} ) f ^{''} ( x ^{*} )}{f ^{'} ( x ^{*} ) ^{2}} ∣ = 1 - \frac{1}{n} < 1$ ，则牛顿迭代法仍有局部收敛性，但重数 $n$ 越高，收敛越慢。

加速收敛法的思路：根据前面的分析，对于重根的情况，我们希望加速收敛速度，既然单根收敛速度较快，所以考虑将要求的 $f$ 的重根转化为另一函数的单根，因为 $x^{*}$ 是 $f (x)$ 的 $n$ 重根，则 $x^{*}$ 是 $f^{'} (x)$ 的 $n - 1$ 重根，令 $μ (x) = \frac{f ( x )}{f ^{'} ( x )}$ ，则 $f$ 的重根等于 $μ$ 的单根，此时对 $μ (x)$ 使用牛顿迭代法求出的单根就是 $f (x)$ 的重根。

牛顿下山法

例：牛顿迭代法求方程 $x^{3} - x - 1 = 0$ 在 $x_{0} = 1.5$ 附近一个根.

解：建立牛顿迭代公式 $x_{k + 1} = x_{k} - \frac{x _{k}^{3} - x _{k} - 1}{3 x _{k}^{2} - 1}$ ，算出 $x_{0} = 1.5$ 、 $x_{1} = 1.34783$ 、 $x_{2} = 1.32520$ 、 $x_{3} = 1.32472$ ，此时是收敛的。

前面的两个收敛定理可知，牛顿迭代法要收敛则要求初始值 $x_{0}$ 选在 $x^{*}$ 附近，这题如果选取初始值 $x_{0} = 0.6$ 的话，则 $x_{0} = 0.6$ 、 $x_{1} = 17.9$ ，算法发散。为防止迭代发散，在迭代的过程中附加一项要求 $∣ f (x_{k + 1}) ∣ < ∣ f (x_{k}) ∣$ ，以保证函数值单调下降，满足这项要求的方法称为下山法。

牛顿下山法：将牛顿法和下山法结合起来使用，在下山法保证函数值稳定下降的前提下，用牛顿法加快收敛速度，就叫做牛顿下山法。具体步骤为：若由 $x_{k}$ 得到的 $x_{k + 1}$ 不能使 $∣ f ∣$ 减小，则在 $x_{k}$ 和 $x_{k + 1}$ 之间找一个更好的点 $\overline{x_{k + 1}}$ ，使得 $∣ f (\overline{x_{k + 1}}) ∣ < ∣ f (x_{k}) ∣$ ，其中 $\overline{x_{k + 1}} = λ x_{k + 1} + (1 - λ) x_{k} = λ (x_{k} - \frac{f ( x _{k} )}{f ^{'} ( x _{k} )}) + (1 - λ) x_{k} = x_{k} - λ \frac{f ( x _{k} )}{f ^{'} ( x _{k} )}$ ， $λ \in [0, 1]$ ，这里的 $λ$ 称为下山因子。

下山因子 $λ$ 的选择是一个逐步探索的过程，一般先取 $λ = 1$ ，此时就是牛顿迭代公式，若满足下山条件 $∣ f (\overline{x_{k + 1}}) ∣ < ∣ f (x_{k}) ∣$ ，则取 $x_{k + 1} = \overline{x_{k + 1}}$ ，再进行下一次迭代；若不满足下山条件，则将 $λ$ 减半计算，再判断是否满足下山条件，依此类推，直到满足为止；若 $λ$ 取得很小时，即 $λ < ε$ 时，下山条件仍不满足，则称”下山失败“，这时需另选初始值重新计算。

接前例，若取 $x_{0} = 0.6$ ，则按牛顿迭代公式求得的迭代值 $\overline{x_{1}} = 17.9$ ，通过反复调整下山因子 $λ$ ，得知当 $λ = \frac{1}{32}$ 时满足下山条件，代入到 $x_{k + 1} = λ \overline{x_{k + 1}} + (1 - λ) x_{k}$ 得 $x_{1} = 1.140625$ 。

正割法

牛顿迭代法的每一步都要计算 $f$ 和 $f^{'}$ ，相当于计算两个函数值，比较费时，如果用 $f$ 的值近似求 $f^{'}$ ，可少算一个函数。具体近似方法是用割线斜率近似切线斜率，即 $f^{'} (x_{k}) \approx \frac{f ( x _{k} ) - f ( x _{k - 1} )}{x _{k} - x _{k - 1}}$ ，如下图：

将近似公式代入牛顿迭代公式中得 $x_{k + 1} = x_{k} - \frac{f ( x _{k} ) ( x _{k} - x _{k - 1} )}{f ( x _{k} ) - f ( x _{k - 1} )}$ ，这个方法就叫做正割法。

正割法不需要求导，但它需要两个初值 $x_{0}$ 和 $x_{1}$ ，当导数难以求得或计算量较大时比较具有优势。

正割法收敛速度比牛顿迭代法慢，且对初值要求同样较高。

三、线性方程组的解法

引入：对于有 $n$ 个未知量 $x_{1}, x_{2}, \dots, x_{n}$ 的线性方程组 $⎩ ⎨ ⎧ a_{11} x_{1} + a_{12} x_{2} + \dots + a_{1 n} x_{n} = b_{1} a_{21} x_{1} + a_{22} x_{2} + \dots + a_{2 n} x_{n} = b_{2} \dots a_{n 1} x_{1} + a_{n 2} x_{2} + \dots + a_{nn} x_{n} = b_{n}$ ，其中 $a_{ij}$ 、 $b_{i}$ 为方程组的系数，可以写成矩阵的形式 $A x = b$ ，当 $D = det A \neq = 0$ 时，根据线性代数中所学的克莱姆法则，方程 $A x = b$ 的解存在且唯一，其中 $x_{i} = \frac{D _{i}}{D}$ ， $D_{i}$ 是将 $A$ 的第 $i$ 列元素以 $b$ 代替的矩阵行列式的值。由于当方程组较多时，克莱姆法则的计算量较大，所以一般用直接法求解。

直接法：如果所有计算都是精确进行的，那么经过有限步算术运算就可以求出方程组精确解的方法。

在实际计算过程中由于舍入误差的存在与影响，直接法一般只能求得方程组的近似解。对于中等规模的线性方程组（n<200），由于直接法的准确性和可靠性高，一般都用直接法求解。直接法主要有高斯消元法和三角分解法。

对于系数为大型稀疏矩阵的方程组，还可以使用迭代法进行求解。

1. 高斯消元法及其改进

高斯消元法

先通过例题了解以前我们如何求解线性方程组：

例： ${x_{1} - x_{2} = 1 3 x_{1} + 2 x_{2} = 8$ .

解：方程1乘以 $- 3$ 加到方程2中，得 ${x_{1} - x_{2} = 1 0 + 5 x_{2} = 5$ ，解出 $x_{2} = 1$ 代入方程1中得 $x_{1} = 2$ 。

我们把上述步骤写成矩阵的形式 $(A, b) = [13 - 1 2 18] \to [10 - 1 5 15] = (A, b)$ ，可见 $A$ 的部分是一个上三角矩阵。于是我们得到了高斯消元法的思路：先将 $A$ 化为上三角矩阵，再回代求解。

其中我们在变换方程组的时候可以作如下两种操作：（1）对换某两个方程的次序；（2）把某一个方程两边同乘一个不为零的常数后加到另一个方程的两边。这两种操作分别对应了 $(A, b)$ 的两种初等行变换。

高斯消元法步骤：为表示方便，我们将线性方程组初始的系数矩阵记为 $A^{(1)} = A = (a_{ij}^{(1)})_{n \times n}$ ，初始的右端项记为 $\overline{b^{(1)}} = \overline{b} = b_{1}^{(1)} ⋮ b_{n}^{(1)}$ ：

（1）设 $a_{11}^{(1)} \neq = 0$ ，计算因子 $m_{i 1} = \frac{a _{i 1}^{(1)}}{a _{11}^{(1)}}$ （ $i = 2, \dots, n$ ），将增广阵第 $i$ 行 $- m_{i 1} \times$ 第 $1$ 行，得到 $(a_{11}^{(1)} 0 a_{12}^{(1)} \dots A^{(2)} a_{1 n}^{(1)} b_{1}^{(1)} \overline{b^{(2)}})$ ，其中 ${a_{ij}^{(2)} = a_{ij}^{(1)} - m_{i 1} a_{1 j}^{(1)} b_{i}^{(2)} = b_{i}^{(1)} - m_{i 1} b_{1}^{(1)}$ （ $i, j = 2, \dots, n$ ）；

（2）依此类推，到第 $k$ 步：设 $a_{kk}^{(k)} \neq = 0$ ，计算因子 $m_{ik} = \frac{a _{ik}^{(k)}}{a _{kk}^{(k)}}$ （ $i = k + 1, \dots, n$ ），则 ${a_{ij}^{(k + 1)} = a_{ij}^{(k)} - m_{ik} a_{kj}^{(k)} b_{i}^{(k + 1)} = b_{i}^{(k)} - m_{ik} b_{k}^{(k)}$ ，经过 $n - 1$ 步就可以将 $A$ 的部分化为上三角；

（3）回代：对于 $a_{11}^{(1)} a_{12}^{(1)} a_{22}^{(2)} \dots \dots \dots a_{1 n}^{(1)} a_{2 n}^{(2)} ⋮ a_{nn}^{(n)} x_{1} x_{2} ⋮ x_{n} = b_{1}^{(1)} b_{2}^{(2)} ⋮ b_{n}^{(n)}$ ，可以看出 $x_{n} = \frac{b _{n}^{(n)}}{a _{nn}^{(n)}}$ 、 $x_{n - 1} = \frac{b _{n - 1}^{(n - 1)} - a _{n - 1, n}^{(n - 1)} \cdot x _{n}}{a _{n - 1, n - 1}^{(n - 1)}}$ ，类推得 $x_{i} = \frac{b _{i}^{(i)} - \sum _{j = i + 1}^{n} a _{ij}^{i} x _{j}}{a _{ii}^{(i)}}$ （ $i = n - 1, \dots, 1$ ）。

定理：若 $A$ 的所有顺序主子式均不为零，则高斯消元法无需换行即可进行到底，且得到唯一解。（顺序主子式：前 $k$ 行前 $k$ 列构成的矩阵）

例：用高斯消元法解 $⎩ ⎨ ⎧ 2 x_{1} + x_{2} - 0.1 x_{3} + x_{4} = 2.7 0.4 x_{1} + 0.5 x_{2} + 4 x_{3} - 8.5 x_{4} = 21.9 0.3 x_{1} - x_{2} + x_{3} + 5.2 x_{4} = - 3.9 x_{1} + 0.2 x_{2} + 2.5 x_{3} - x_{4} = 9.9$ .

解： $(A, b) = 2 0.4 0.3 1 1 0.5 - 1 0.2 - 0.1 41 2.5 1 - 8.5 5.2 - 1 2.7 21.9 - 3.9 9.9 \to 2000 1 0.3 - 1.15 - 0.3 - 0.1 4.02 1.015 2.55 1 - 8.7 5.05 - 1.5 2.7 21.36 - 4.305 8.55 \to 2000 1 0.3 00 - 0.1 4.02 16.425 6.57 1 - 8.7 - 28.3 - 10.2 2.7 21.36 77.576 29.91$

$\to 2000 1 0.3 00 - 0.1 4.02 16.425 0 1 - 8.7 - 28.3 - 1.12 2.7 21.36 77.576 - 1.12$ ，回代可得 $x = x_{1} x_{2} x_{3} x_{4} = 123 - 1$ .（PPT中计算错误，过程有错）

注意：当 $∣ a_{kk} ∣$ 如果比同一列的元素小得多的话，会导致很大的舍入误差。因此在使用高斯消元法时，除了要避免主元素 $a_{kk}^{(k)} = 0$ 外，还要避免主元素绝对值很小。

例：用单精度解方程组 ${1 0^{- 9} x_{1} + x_{2} = 1 x_{1} + x_{2} = 2$ .

解： $m_{21} = \frac{a _{21}}{a _{11}} = 1 0^{9}$ ， $a_{22} = 1 - m_{21} \times 1 = 0.0000 00001 \times 1 0^{9} - 1 0^{9} = - 1 0^{9}$ ， $b_{2} = 2 - m_{21} \times 1 = - 1 0^{9}$ ，即 $[1 0^{- 9} 1 1112] \to [1 0^{- 9} 0 1 - 1 0^{9} 1 - 1 0^{9}] \Rightarrow x = [x_{1} x_{2}] = [01]$ ，此结果中 $x_{1}$ 显然与正确结果偏差较大。

全主元消元法

全主元消元法是对高斯消元法的改进：

（1）第一次消元时，在系数矩阵 $A$ 中寻找绝对值最大的元素作为主元素来代替 $a_{11}^{(1)}$ ，设 $∣ a_{i_{1} j_{1}} ∣ = 1 \leq i \leq n, 1 \leq j \leq n max ∣ a_{ij} ∣$ ，再交换 $i_{1}$ 行与第 $1$ 行，交换 $j_{1}$ 列与第 $1$ 列，这样就将 $a_{i_{1} j_{1}}$ 换到了左上角的 $a_{11}^{(1)}$ 位置，再进行消元；

（2）依此类推，到第 $k$ 次消元时为，主元素 $a_{i_{k} j_{k}}$ 从系数矩阵右下角 $(n - k + 1)$ 阶子矩阵中选取 $∣ a_{i_{k} j_{k}} ∣ = k \leq i \leq n, k \leq j \leq n max ∣ a_{ij} ∣$ ，交换至主元素位置并消元。

换行不影响方程组求解，但是换列交换了未知数的系数，所以必须把次序记录下来，解完方程组后再把次序换回来。

算法的每一步都选取绝对值最大的元素作为主元素，以保证 $∣ m_{ik} ∣ < 1$ ，这样精度比较高，算法比较稳定。但是该算法需要花大量时间去记录列的次序。

列主元消元法

列主元消元法是对全主元消元法的改进，它省去了换列的步骤，每次选一列中绝对值最大的元素作为主元素： $∣ a_{i_{k}, k} ∣ = k \leq i \leq n max ∣ a_{ik} ∣ \neq = 0$ ，以 $a_{i_{k}, k}$ 代替 $a_{kk}^{k}$ ，可以避免列的交换，而仅交换行。

例：用列主元消元法解方程组 $⎩ ⎨ ⎧ 4 x_{1} - x_{2} + x_{3} = 5 - 18 x_{1} + 3 x_{2} - x_{3} = - 15 x_{1} + x_{2} + x_{3} = 6$ .

解： $(A, b) = 4 - 18 1 - 1 31 1 - 1 1 5 - 15 6 \to - 18 41 3 - 1 1 - 1 11 - 15 56 \to - 18 00 3 - \frac{1}{3} \frac{7}{6} - 1 \frac{7}{9} \frac{17}{18} - 15 \frac{5}{3} \frac{31}{6}$

$\to - 18 00 3 \frac{7}{6} - \frac{1}{3} - 1 \frac{17}{18} \frac{7}{9} - 15 \frac{31}{6} \frac{5}{3} \to - 18 00 3 \frac{7}{6} 0 - 1 \frac{17}{18} \frac{22}{21} - 15 \frac{31}{6} \frac{22}{7}$ ，回代得 $x = x_{1} x_{2} x_{3} = 123$ 。

列主元消元法仍然有精度的问题，它没有全主元消元法稳定。

比如（同前面高斯消元法的例题）： $[1 0^{- 9} 1 1112] \to [1 1 0^{- 9} 1121] \to [101121]$ ，解出 $x = [x_{1} x_{2}] = [11]$ ，显然这次离精确值非常靠近了。

但是如果将第一行乘以 $1 0^{9}$ 得到完全等价的方程组 $[11 1 0^{9} 1 1 0^{9} 2] \to [10 1 0^{9} 1 0^{- 9} 1 0^{9} 1 0^{- 9}]$ ，解出 $x = [x_{1} x_{2}] = [01]$ ，这时又出现了精度问题了。

标度化列主元消元法

标度化列主元消元法对列主元消元法的改进：对每一行计算 $s_{i} = 1 \leq j \leq n max ∣ a_{ij} ∣$ ，为节省时间， $s_{i}$ 只在初始时计算一次，以后每次考虑子列 $a_{kk} ⋮ a_{nk}$ 中 $∣ \frac{a _{ik}}{s _{i}} ∣$ 最大的 $a_{ik}$ 为主元。

标度化列主元消元法稳定性介于列主元消元法与全主元消元法之间。

例：用标度化列主元消元法解 ${x_{1} + 1 0^{9} x_{2} = 1 0^{9} x_{1} + x_{2} = 2$ .

解：写出 $(A, b) = [11 1 0^{9} 1 1 0^{9} 2]$ ，找出每一行系数矩阵 $A$ 部分中的最大值 $s_{1} = 1 0^{9}$ 、 $s_{2} = 1$ ，对每一行都作标度化处理，即每一行所有元素都除以对应的 $s_{i}$ ，得 $[1 0^{- 9} 1 1112]$ ，此时再用列主元消元法进行计算即可，比如 $∣ \frac{a _{11}}{s _{1}} ∣ = 1 0^{- 9}$ 与 $∣ \frac{a _{21}}{s _{2}} ∣ = 1$ 相比，显然应取 $a_{21}$ 为主元素，于是交换1、2行， $\dots$ （参考前面例题），解出 $x = [x_{1} x_{2}] = [11]$ ，所得解很接近精确解了。

高斯-若尔当消元法

高斯-若尔当消元法与之前的高斯消元法不同，它是将 $a_{kk}$ 变为 $1$ ，与此同时 $a_{kj}$ 变为 $\frac{a _{kj}}{a _{kk}}$ （ $j = k, k + 1, \dots, n + 1$ ），再将除 $a_{kk}$ 外的第 $k$ 列的所有元素都变为 $0$ ，即 $a_{ij}$ 等于 $a_{ij} - a_{ik} a_{kj}$ （ $i \neq = k, j = k + 1, \dots, n + 1$ ），经过 $n$ 步计算后得到 $100 ⋮ 0 010 ⋮ 0 \dots \dots \dots ⋱ \dots 000 ⋮ 1 a_{1, n + 1} a_{2, n + 1} a_{3, n + 1} ⋮ a_{n, n + 1}$ ，这样就不需要回代即可得出解。

其实这个过程就相当于 $A x = b \Rightarrow E x = A^{- 1} b$ ，其中 $E$ 为单位阵。

例：用高斯-若尔当消元法与全主元消元法相结合解方程 $⎩ ⎨ ⎧ - 3 x_{1} + 8 x_{2} + 5 x_{3} = 6 2 x_{1} - 7 x_{2} + 4 x_{3} = 9 x_{1} + 9 x_{2} - 6 x_{3} = 1$ .

解：写出 $- 3 21 8 - 7 9 54 - 6 691$ ，按照全主元消元法找出绝对值最大的元素为 $9$ ，按照高斯-若尔当消元法将第三行元素都除以 $9$ 得 $- 3 2 \frac{1}{9} 8 - 7 1 54 - \frac{6}{9} 69 \frac{1}{9}$ ，将 $1$ 所在列变为零得 $- \frac{35}{9} \frac{25}{9} \frac{1}{9} 001 \frac{31}{3} - \frac{2}{3} - \frac{6}{9} \frac{46}{9} \frac{88}{9} \frac{1}{9}$ ；再从除 $1$ 所在列外的系数矩阵 $A$ 部分 $[- \frac{35}{9} \frac{25}{9} \frac{31}{3} - \frac{2}{3}]$ 中找出最大的元素是 $\frac{31}{3}$ ，重复前述步骤得 $- \frac{35}{93} \frac{25}{9} \frac{1}{9} 001 1 - \frac{2}{3} - \frac{6}{9} \frac{46}{93} \frac{88}{9} \frac{1}{9} \to - \frac{35}{93} \frac{235}{93} - \frac{13}{93} 001100 \frac{46}{93} \frac{940}{93} \frac{41}{93} \to - \frac{35}{93} 1 - \frac{13}{93} 001100 \frac{46}{93} \frac{940}{235} \frac{41}{93} \to 010001100241$ ，故解为 $x_{3} = 2$ 、 $x_{1} = 4$ 、 $x_{2} = 1$ 。

可将，将高斯-若尔当消元法与全主元消元法相结合，可以得到不用换行、不用换列，也不用回代，并且精度很高的算法。

例：用高斯-若尔当消元法与；列主元消元法结合求矩阵 $A$ 的逆矩阵， $A = 123245356$ .

解：原理是 $A^{- 1} (A ∣ E) = (A^{- 1} A ∣ A^{- 1} E) = (E ∣ A^{- 1})$ ，即线性代数中的用增广矩阵求逆的方法。

写出 $123245356100010001$ ，找出第一列中最大的元素是 $3$ ，交换第1行和第3行得 $321542653001010100$ ，第一行元素除以 $3$ （标准化）再将第一列其他位置化为 $0$ （消元）得 $100 \frac{5}{3} \frac{2}{3} \frac{1}{3} 211001010 \frac{1}{3} - \frac{2}{3} - \frac{1}{3}$ ，重复前述步骤，最后得到 $100010001 1 - 3 2 - 3 3 - 1 2 - 1 0$ .

2. 三角分解法

道立特(Doolittle)分解法

引入：高斯消元法实际上是对系数矩阵和常向量进行初等变换，而初等变换可以看作是左乘初等矩阵，即 $L_{n - 1} L_{n - 2} \dots L_{1} (A, b) = a_{11}^{(1)} a_{12}^{(1)} a_{22}^{(2)} \dots \dots \dots a_{1 n}^{(1)} a_{2 n}^{(2)} ⋮ a_{nn}^{(n)} b_{1}^{(1)} b_{2}^{(2)} ⋮ b_{n}^{(n)}$ ，或 $L_{n - 1} L_{n - 2} \dots L_{1} A x = L_{n - 1} L_{n - 2} \dots L_{1} b$ ，令 $L_{n - 1} L_{n - 2} \dots L_{1} A = U$ ，则 $A$ 可以写成 $A = L_{1}^{- 1} L_{2}^{- 1} \dots L_{n - 1}^{- 1} U = LU$ ，其中 $L_{k}^{- 1} = 100000 ⋱ 0000 1 m_{k + 1, k} ⋮ m_{n, k} ⋱ ⋱ 0 ⋱ 0 1$ ，而 $L_{1}^{- 1} L_{2}^{- 1} \dots L_{n - 1}^{- 1} = L$ 也是一个下三角阵且对角线元素都是 $1$ ，称为单位下三角阵，而 $U = a_{11}^{(1)} a_{12}^{(1)} a_{22}^{(2)} \dots \dots \dots a_{1 n}^{(1)} a_{2 n}^{(2)} ⋮ a_{nn}^{(n)}$ 。

$LU$ 分解又叫三角分解，当 $L$ 为单位下三角而 $U$ 为非奇异上三角的分解称为道立特（Doolittle）分解；当 $L$ 为一般下三角而 $U$ 为单位上三角的分解称为克路特（Crout）分解。

定理：若 $A$ 的所有顺序主子式均不为零，则 $A$ 的 $LU$ 分解唯一（其中 $L$ 为单位下三角阵）。

证明：由前面可知 $LU$ 分解存在，只需证明其唯一性。假设不唯一，设 $A = L_{1} U_{1} = L_{2} U_{2}$ ，则 $U_{1} = L_{1}^{- 1} L_{1} U_{1} = L_{1}^{- 1} L_{2} U_{2}$ ，两端同时右乘 $U_{2}^{- 1}$ 得 $U_{1} U_{2}^{- 1} = L_{1}^{- 1} L_{2}$ ，观察可知， $U_{1}$ 和 $U_{2}^{- 1}$ 是上三角，相乘仍为上三角，而 $L_{1}^{- 1}$ 和 $L_{2}$ 是下三角，相乘仍为下三角，显然只有当等式两端都是单位阵时才能成立，即 $U_{1} = U_{2}$ 、 $L_{1} = L_{2}$ 。

Doolittle分解法的步骤：

（1）三角分解：通过比较法直接导出 $L$ 和 $U$ 的计算公式，观察 $a_{11} ⋮ ⋮ a_{n 1} \dots \dots a_{1 n} ⋮ ⋮ a_{nn} = 1 l_{21} ⋮ l_{n 1} 1 \dots \dots ⋱ \dots 1 u_{11} \dots \dots u_{1 n} ⋮ ⋮ u_{nn}$ ，可看出 $a_{ij} = \sum_{k = 1}^{m i n (i, j)} l_{ik} u_{kj}$ ，若 $L$ 的 $r - 1$ 列与 $U$ 的前 $r - 1$ 行已经算好，则 $u_{r, j} = a_{r, j} - \sum_{k = 1}^{r - 1} l_{r k} \cdot u_{kj}$ （ $j = r, \dots, n$ ）， $l_{i r} = \frac{1}{u _{rr}} (a_{i r} - \sum_{k = 1}^{r - 1} l_{ik} u_{k r})$ （ $i = r + 1, \dots, n$ ）；

比如 $1 \cdot u_{1 j} = a_{1 j} \Rightarrow u_{1 j} = a_{1 j}$ （ $j = 1, 2, \dots, n$ ）， $l_{i 1} \cdot u_{11} = a_{i 1} \Rightarrow l_{i 1} = \frac{a _{i 1}}{u _{11}}$ （ $i = 2, \dots, n$ ）， $l_{21} \cdot u_{1 j} + u_{2 j} = a_{2 j} \Rightarrow u_{2 j} = a_{2 j} - l_{i 1} u_{1 j}$ （ $j = 2, \dots, n$ ）， $l_{i 1} \cdot u_{12} + l_{i 2} \cdot u_{22} = a_{i 2} \Rightarrow l_{i 2} = \frac{a _{i 2} - l _{i 1} u _{12}}{u _{22}}$ （ $i = 3, \dots, n$ ）；可以看出 $l_{r 1} \cdot u_{1 j} + l_{r 2} \cdot u_{2 j} + \dots + l_{r, r - 1} \cdot u_{r - 1, j} + u_{r, j} = a_{r, j}$ ，可推出 $u_{r, j} = a_{r, j} - \sum_{k = 1}^{r - 1} l_{r k} \cdot u_{kj}$ （ $j = r, \dots, n$ ）；同理，由 $l_{i 1} \cdot u_{1 r} + l_{i 2} \cdot u_{2 r} + \dots + l_{i, r - 1} \cdot u_{r - 1, r} + l_{i r} \cdot u_{rr} = a_{i r}$ ，可得 $l_{i r} = \frac{a _{i r} - \sum _{k = 1}^{r - 1} l _{ik} u _{k r}}{u _{rr}}$ （ $i = r + 1, \dots, n$ ）。

（2）求出 $A$ 的 $LU$ 分解后， $A x = b \to LUx = b$ ，设 $Ux = y$ ，得 $L y = b$ ，先解出 $y$ 再解 $Ux = y$ 即可。

对于 $L y = b$ ，即 $1 l_{21} ⋮ l_{n 1} 1 ⋮ \dots ⋱ \dots 1 y_{1} y_{2} ⋮ y_{n} = b_{1} b_{2} ⋮ b_{n}$ ，则 $y_{i} = b_{i} - \sum_{k = 1}^{i - 1} l_{ik} y_{k}$ （ $i = 1, 2, \dots, n$ ）；

对于 $Ux = y$ ，即 $u_{11} u_{12} u_{22} \dots \dots ⋱ u_{1 n} u_{2 n} ⋮ u_{nn} x_{1} x_{2} ⋮ x_{n} = y_{1} y_{2} ⋮ y_{n}$ ，则 $x_{i} = \frac{1}{u _{ii}} (y_{i} - \sum_{k = i + 1}^{n} u_{ik} \cdot x_{k})$ （ $i = n, n - 1, \dots, 1$ ）。

例：用Doolittle分解法求解方程组 $⎩ ⎨ ⎧ 2 x_{1} + x_{2} + x_{3} = 4 x_{1} + 3 x_{2} + 2 x_{3} = 6 x_{1} + 2 x_{2} + 2 x_{3} = 5$ .

解：对于 $211132122 = 1 l_{21} l_{31} 1 l_{32} 1 u_{11} u_{12} u_{22} u_{13} u_{23} u_{33}$ ，根据前面的分析知 $u_{1 j} = a_{1 j}$ ，故 $u_{11} = 2$ 、 $u_{12} = 1$ 、 $u_{13} = 1$ ，同时 $l_{i 1} = \frac{a _{i 1}}{u _{11}}$ ，故 $l_{21} = 0.5$ 、 $l_{31} = 0.5$ ，同理 $u_{22} = a_{22} - l_{21} u_{12} = 3 - 0.5 = 2.5$ ， $u_{23} = a_{23} - l_{21} u_{13} = 2 - 0.5 = 1.5$ ，同时 $l_{32} = \frac{a _{32} - l _{31} u _{12}}{u _{22}} = \frac{2 - 0.5}{2.5} = 0.6$ ，同理 $u_{33} = a_{33} - l_{31} u_{13} - l_{32} u_{23} = 0.6$ ；

$L y = b$ ，即 $1 0.5 0.5 1 0.6 1 y_{1} y_{2} y_{3} = 465 \Rightarrow y_{1} y_{2} y_{3} = 44 0.6$ ， $Ux = y$ ，即 $200 1 2.5 0 1 1.5 0.6 x_{1} x_{2} x_{3} = 44 0.6 \Rightarrow x_{1} x_{2} x_{3} = 111$ .

平方根法-对称正定系数矩阵

定义：对于矩阵 $A = (a_{ij})_{n \times n}$ ，如果 $a_{ij} = a_{ji}$ ，则称 $A$ 为对称阵。

定义：对于矩阵 $A = (a_{ij})_{n \times n}$ ，如果 $x^{T} A x > 0$ 对任意非零向量 $x$ 都成立，则称 $A$ 为正定阵。

回顾对称正定阵的几个性质：（1） $A^{- 1}$ 也正定，且 $a_{ii} > 0$ ；（2） $A$ 的顺序主子阵 $A_{k}$ 也对称正定；（3） $A$ 的特征值 $λ_{i} > 0$ ；（4） $A$ 的全部顺序主子式 $det (A_{k}) > 0$ ；

定理：设矩阵 $A$ 对称正定，则存在非奇异下三角阵 $L \in R^{n \times n}$ 使得 $A = L L^{T}$ ，若限定 $L$ 对角元为正，则分解唯一。

证明：先对 $A$ 作 $LU$ 分解，即 $A = LU$ ；

先证明 $u_{ii}$ 都大于零，将 $A = LU$ 写成分块运算的形式 $A_{k} = L_{k} U_{k}$ （ $k = 1, 2, \dots, n$ ），根据正定的性质， $A$ 的各阶主子式的行列式值都大于零，所以 $det A_{k} = det L_{k} \cdot det U_{k} \Rightarrow \prod_{i = 1}^{k} u_{ii} > 0$ ，所以 $u_{ii} > 0$ ；

进一步将 $U$ 分解为 $U = D U^{'} = u_{11} u_{22} ⋱ u_{nn} 1 \frac{u _{12}}{u _{11}} 1 \dots \dots ⋱ \frac{u _{1 n}}{u _{11}} \frac{u _{2 n}}{u _{22}} ⋮ 1$ ，则 $A = L D U^{'}$ ，由于 $A$ 对称所以 $A^{T} = A$ ，故 $L D U^{'} = U^{' T} D^{T} L^{T}$ ，由于 $D$ 是对角阵 $D = D^{T}$ ，观察可知 $L = U^{' T}$ ，即 $A = L D^{T} L^{T}$ ，由于前面证明了对称正定矩阵的 $u_{ii} > 0$ ，于是可以将 $D$ 进一步分解为 $D = D^{'} D^{'}$ ，其中 $D^{'} = u_{11} ⋱ u_{nn}$ ，令 $L^{*} = L D^{'}$ ，则 $A = L D^{'} D^{' T} L^{T} = L^{*} L^{*}^{T}$ 。

平方根法的步骤（要求 $A$ 对称正定）：

（1）利用 $A = L L^{T}$ 分解，对于 $l_{11} ⋮ l_{n 1} ⋱ \dots l_{nn} l_{11} \dots ⋱ l_{n 1} ⋮ l_{nn} = a_{11} ⋮ a_{n 1} \dots \dots a_{1 n} ⋮ a_{nn}$ ，当 $i > j$ 时 $l_{ij} = \frac{a _{ij} - \sum _{k = 1}^{j - 1} l _{ik} \cdot l _{jk}}{l _{jj}}$ ，当 $i = j$ 时 $l_{ii} = (a_{ii} - \sum_{k = 1}^{i - 1} l_{ik}^{2})$ ；

比如 $l_{11} \cdot l_{11} = a_{11} \Rightarrow l_{11} = a_{11}$ ，类似地第一列有 $l_{i 1} \cdot l_{11} = a_{i 1} \Rightarrow l_{i 1} = \frac{a _{i 1}}{l _{11}}$ （ $i = 2, \dots, n$ ），依此类推有 $l_{i 1} \cdot l_{j 1} + l_{i 2} \cdot l_{j 2} + \dots + l_{i, j} \cdot l_{j, j} = a_{ij} \Rightarrow l_{ij} = \frac{a _{ij} - \sum _{k = 1}^{j - 1} l _{ik} \cdot l _{jk}}{l _{jj}}$ （ $i > j$ ），同理 $l_{i 1} \cdot l_{i 1} + l_{i 2} \cdot l_{i 2} + \dots + l_{ii} \cdot l_{ii} = a_{ii} \Rightarrow l_{ii} = (a_{ii} - \sum_{k = 1}^{i - 1} l_{ik}^{2})$ （ $i = j$ ）。

同时可以发现 $l_{ii}^{2} = a_{ii} - \sum_{k = 1}^{i - 1} l_{ik}^{2} > 0 \Rightarrow a_{ii} > \sum_{k = 1}^{i - 1} l_{ik}^{2} \Rightarrow ∣ l_{ik} ∣ < a_{ii} < 1 \leq i \leq n max a_{ii}$ ，因此 $L$ 的元素有界且 $l_{ii} > 0$ ，故平方根法是比较稳定的方法，且精度比较好。

（2）求出 $L$ 后， $A x = b \to L L^{T} x = b$ ，令 $L^{T} x = y$ ，先解出 $L y = b$ ，再求解 $L^{T} x = y$ 即可，公式为 $y_{i} = \frac{1}{l _{ii}} (b_{i} - \sum_{k = 1}^{i - 1} l_{ik} y_{k})$ （ $i = 1, 2, \dots, n$ ）， $x_{i} = \frac{1}{l _{ii}} (y_{i} - \sum_{k = i + 1}^{n} l_{ik} x_{k})$ （ $i = n, n - 1, \dots, 1$ ）。

改进平方根法（要求 $A$ 对称正定）：

（1）由于使用平方根法分解矩阵 $A$ 时，需要做 $n$ 次开方运算，为避免开方，可以直接用 $A = L D L^{T}$ 分解来运算： $A = L D L^{T} \Rightarrow a_{11} a_{21} ⋮ a_{n 1} a_{12} a_{22} ⋮ a_{n 2} \dots \dots ⋱ \dots a_{1 n} a_{2 n} ⋮ a_{nn} = 1 l_{21} ⋮ l_{n 1} 1 ⋮ l_{n 2} ⋱ \dots 1 d_{11} d_{22} ⋱ d_{nn} 1 l_{21} 1 \dots \dots ⋱ l_{1 n} l_{2 n} ⋮ 1$ ，将 $L D$ 乘到一起得 $a_{11} a_{21} ⋮ a_{n 1} a_{12} a_{22} ⋮ a_{n 2} \dots \dots ⋱ \dots a_{1 n} a_{2 n} ⋮ a_{nn} = d_{11} l_{21} d_{11} ⋮ l_{n 1} d_{11} d_{22} ⋮ l_{n 2} d_{22} ⋱ \dots d_{nn} 1 l_{21} 1 \dots \dots ⋱ l_{1 n} l_{2 n} ⋮ 1$ ，写出 $a_{ij} = l_{i 1} d_{11} l_{1 j} + l_{i 2} d_{22} l_{2 j} + \dots + l_{i, j - 1} d_{j - 1, j - 1} l_{j - 1, j} + l_{ij} d_{jj} = \sum_{k = 1}^{j - 1} l_{ik} d_{kk} l_{jk} + l_{ij} d_{jj}$ ，特别地， $a_{ii} = \sum_{k = 1}^{i - 1} l_{ik}^{2} d_{kk} + d_{ii}$ ，令 $l_{ik} d_{kk} = T_{ik}$ ，则 ${a_{ij} = \sum_{k = 1}^{j - 1} T_{ik} l_{jk} + T_{ij} a_{ii} = \sum_{k = 1}^{i - 1} T_{ik} l_{ik} + d_{ii}$ ，于是 $⎩ ⎨ ⎧ d_{11} = a_{11} T_{ij} = a_{ij} - \sum_{k = 1}^{j - 1} T_{ik} l_{jk} l_{ij} = \frac{T _{ij}}{d _{jj}} d_{ii} = a_{ii} - \sum_{k = 1}^{i - 1} T_{ik} l_{ik}$ ；

（2）对于 $A x = b \to L D L^{T} x = b \to {L y = b D L^{T} x = y$ ；对于 $L y = b$ ，即 $1 l_{21} ⋮ l_{n 1} 1 ⋮ \dots ⋱ \dots 1 y_{1} y_{2} ⋮ y_{n} = b_{1} b_{2} ⋮ b_{n}$ ，则 $y_{i} = b_{i} - \sum_{k = 1}^{i - 1} l_{ik} y_{k}$ （ $i = 1, 2, \dots, n$ ）；对于 $D L^{T} x = y$ ，即 $d_{11} d_{22} ⋱ d_{nn} 1 l_{12} 1 \dots \dots ⋱ l_{1 n} l_{2 n} ⋮ 1 x_{1} x_{2} ⋮ x_{n} = y_{1} y_{2} ⋮ y_{n}$ ，算出 $d_{11} d_{11} l_{12} d_{22} \dots \dots ⋱ d_{11} l_{1 n} d_{22} l_{2 n} ⋮ d_{nn} x_{1} x_{2} ⋮ x_{n} = y_{1} y_{2} ⋮ y_{n}$ ，则 $x_{n} = \frac{y _{n}}{d _{nn}}$ ，而 $y_{i} = d_{ii} x_{i} + d_{ii} l_{i, i + 1} x_{i + 1} + \dots + d_{ii} l_{in} x_{n}$ 故 $x_{i} = \frac{y _{i}}{d _{ii}} - \sum_{k = i + 1}^{n} l_{ik} x_{k}$ （ $i = n, n - 1, \dots, 2, 1$ ，并且由于 $L$ 转置后仍使用原来的下标，所以这里的 $l_{ik} = l_{ki}$ ）。

例：用改进平方根法求解方程组 $A x = b$ ，其中 $A = 121 - 3 250 - 5 10141 - 3 - 5 115$ ， $b = 12168$ .

解：按照公式有 $i = 1 : i = 2 : i = 3 : i = 4 : d_{11} = 1 T_{21} = 2, l_{21} = 2, d_{22} = 1 T_{31} = 1, T_{32} = - 2, l_{31} = 1, l_{32} = - 2, d_{33} = 9 T_{41} = - 3, T_{42} = 1, T_{43} = 6, l_{41} = - 3, l_{42} = 1, l_{43} = \frac{2}{3}, d_{44} = 1$ ，再根据公式可以解出 ${y_{1} = 1, y_{2} = 0, y_{3} = 15, y_{4} = 1 x_{1} = 1, x_{2} = 1, x_{3} = 1, x_{4} = 1$ 。

追赶法-解三对角线方程组

追赶法：对于三对角线方程组 $b_{1} a_{2} c_{1} b_{2} ⋱ c_{2} ⋱ a_{n - 1} ⋱ b_{n - 1} a_{n} c_{n - 1} b_{n} x_{1} x_{2} ⋮ ⋮ x_{n} = d_{1} d_{2} ⋮ ⋮ d_{n}$ ，可以使用追赶法：

（1）先将 $A$ 分解为 $A = LU = α_{1} γ_{1} ⋱ ⋱ ⋱ γ_{n} α_{n} 1 β_{1} ⋱ ⋱ ⋱ β_{n - 1} 1 = α_{1} γ_{2} α_{1} β_{1} γ_{2} β_{1} + α_{2} ⋱ α_{2} β_{2} ⋱ γ_{i} ⋱ γ_{i} β_{i - 1} + α_{i} ⋱ α_{i} β_{i} ⋱ γ_{n - 1} ⋱ γ_{n - 1} β_{n - 2} + α_{n - 1} γ_{n} α_{n - 1} β_{n - 1} γ_{n} β_{n - 1} + α_{n}$ ，对比原来的 $A = b_{1} a_{2} c_{1} b_{2} ⋱ c_{2} ⋱ a_{n - 1} ⋱ b_{n - 1} a_{n} c_{n - 1} b_{n}$ ，可得 $⎩ ⎨ ⎧ γ_{i} = a_{i} β_{i} = c_{i} / α_{i} α_{1} = b_{1} α_{i} = b_{i} - γ_{i} β_{i - 1} i = 2, 3, \dots, n i = 1, \dots, n - 1 i = 2, 3, \dots, n$ ；

（2）追：对于 $A x = LUx = d$ ，求解 $L y = d$ ，即 $α_{1} γ_{1} ⋱ ⋱ ⋱ γ_{n} α_{n} y_{1} y_{2} ⋮ y_{n} = d_{1} d_{2} ⋮ d_{n}$ ，得出 ${y_{1} = \frac{d _{1}}{α _{1}} y_{i} = \frac{( d _{i} - γ _{i} y _{i - 1} )}{α _{i}} i = 2, \dots, n$ ；

（3）赶：再求解 $Ux = y$ ，即 $1 β_{1} ⋱ ⋱ 1 β_{n - 1} 1 x_{1} x_{2} ⋮ x_{n} = y_{1} y_{2} ⋮ y_{n}$ ，得出 ${x_{n} = y_{n} x_{i} = y_{i} - β_{i} x_{i + 1} i = n - 1, \dots, 1$ 。

由第二步可知该方法要求 $α_{i}$ 不等于零，并不是任何三对角阵都可以使用此方法求解。

定理（充分条件）：若 $A$ 为三对角线矩阵，并且满足 $⎩ ⎨ ⎧ ∣ b_{1} ∣ > ∣ c_{1} ∣ > 0 ∣ b_{i} ∣ \geq ∣ a_{i} ∣ + ∣ c_{i} ∣ 且 c_{i} \neq = 0, a_{i} \neq = 0 ∣ b_{n} ∣ > ∣ a_{n} ∣ > 0 i = 2, 3 \dots, n - 1$ ，则追赶法可以求解以 $A$ 为系数矩阵的方程组。

该定理要求三条对角线上都不能有零元素，对于方程组来说太过苛刻，因为它只是充分条件而非必要条件，所以可以有例外的情况：如果 $A$ 是严格对角优势的三对角线矩阵，则不要求三对角线上的所有元素非零。

根据不等式 $∣ β_{i} ∣ < 1$ ， $∣ b_{i} ∣ - ∣ a_{i} ∣ < ∣ b_{i} - γ_{i} β_{i - 1} ∣ < ∣ b_{i} ∣ + ∣ a_{i} ∣$ 可知，分解的过程中，矩阵元素不会过分增大，算法保证稳定。

3. 迭代法求方程组的近似解

迭代原理与误差分析

迭代法思路：与解 $f (x) = 0$ 的不动点迭代类似，将 $A x = b$ 等价改写为 $x = B x + f$ 的形式，建立迭代 $x^{(k + 1)} = B x^{(k)} + f$ ，从初值 $x^{(0)}$ 出发，得到序列 ${x^{(k)}}$ 。迭代法精度可控，特别适用于求解系数为大型稀疏矩阵的方程组。

例：用迭代法解方程组 $⎩ ⎨ ⎧ 8 x_{1} - 3 x_{2} + 2 x_{3} = 20 4 x_{1} + 11 x_{2} - x_{3} = 33 6 x_{1} + 3 x_{2} + 12 x_{3} = 36$ .

解：把 $A x = b$ 改写成 $x = B x + f$ 的形式 $⎩ ⎨ ⎧ x_{1} = \frac{1}{8} (3 x_{2} - 2 x_{3} + 20) x_{2} = \frac{1}{11} (- 4 x_{1} + x_{3} + 33) x_{3} = \frac{1}{12} (- 6 x_{1} - 3 x_{2} + 36)$ ，即 $x_{1} x_{2} x_{3} = 0 - \frac{4}{11} - \frac{1}{2} \frac{3}{8} 0 - \frac{1}{4} - \frac{1}{4} \frac{1}{11} 0 \cdot x_{1} x_{2} x_{3} + \frac{5}{2} 33$ ，建立迭代格式 $x^{(k + 1)} = B x^{(k)} + f$ ，得 $⎩ ⎨ ⎧ x_{1}^{(k + 1)} = \frac{1}{8} (3 x_{2}^{(k)} - 2 x_{3}^{(k)} + 20) x_{2}^{(k + 1)} = \frac{1}{11} (- 4 x_{1}^{(k)} + x_{3}^{(k)} + 33) x_{3}^{(k + 1)} = \frac{1}{12} (- 6 x_{1}^{(k)} - 3 x_{2}^{(k)} + 36)$ （矩阵形式略），初始向量取 $(0, 0, 0)^{T}$ ，迭代得

$x^{(k)}$ $x_{1}$ $x_{2}$ $x_{3}$

$x^{(0)}$ 0 0 0

$x^{(1)}$ 2.5 3 3

$x^{(2)}$ 2.875 2.364 1

$x^{(3)}$ 3.1365 2.0455 0.9715

$x^{(4)}$ 3.02418 1.94777 0.92038

$x^{(5)}$ 3.000318 1.983965 1.000963

$x^{(6)}$ 2.993746 1.999971 1.003849

$x^{(7)}$ 2.999026 2.002624 1.003134

直接解方程得知精确解为 $x^{*} = (3, 2, 1)^{T}$ ，迭代七次后的误差 $e^{(7)} = x^{(7)} - x^{*} = (- 0.000974, 0.002624, 0.003134)^{T}$ ，可见用此迭代格式进行迭代，得到的向量序列是逐步逼近方程组精确解的。

但是不是所有迭代序列都可以逼近精确值的，如果用 $⎩ ⎨ ⎧ x_{1} = 9 x_{1} - 3 x_{2} + 2 x_{3} - 20 x_{2} = 4 x_{1} + 12 x_{2} - x_{3} - 33 x_{3} = 6 x_{1} + 3 x_{2} + 13 x_{3} - 36$ ，取 $x^{(0)} = (0, 0, 0)^{T}$ ，得 $x^{(1)} = (- 20, - 33, - 36)^{T}$ ， $x^{(2)} = (- 173, - 473, - 723)^{T}$ ，显然，这个迭代是发散的。

$x^{(k)}$	$x_{1}$	$x_{2}$	$x_{3}$
$x^{(0)}$	0	0	0
$x^{(1)}$	2.5	3	3
$x^{(2)}$	2.875	2.364	1
$x^{(3)}$	3.1365	2.0455	0.9715
$x^{(4)}$	3.02418	1.94777	0.92038
$x^{(5)}$	3.000318	1.983965	1.000963
$x^{(6)}$	2.993746	1.999971	1.003849
$x^{(7)}$	2.999026	2.002624	1.003134

为了更好地度量误差，需要引入向量范数和矩阵范数的概念。

定义： $R^{n}$ 空间的向量范数 $∣∣ \cdot ∣∣$ 对任意 $x, y \in R^{n}$ 满足下列条件：（1）正定性： $∣∣ x ∣∣ \geq 0$ ， $∣∣ x ∣∣ = 0 \Leftrightarrow x = 0$ ；（2）齐次性： $∣∣ αx ∣∣ = ∣ α ∣ \cdot ∣∣ x ∣∣$ 对任意 $α \in C$ 成立；（3）三角不等式： $∣∣ x + y ∣∣ \leq ∣∣ x ∣∣ + ∣∣ y ∣∣$ 。

最常用的范数是 $L_{p}$ 范数（p-范数）： $∣∣ x ∣ ∣_{p} = (\sum_{i = 1}^{n} ∣ x_{i} ∣^{p})^{\frac{1}{p}}$ ，当 $p$ 取1、2、 $\infty$ 时分别是： $∣∣ x ∣ ∣_{1} = \sum_{i = 1}^{n} ∣ x_{i} ∣$ ， $∣∣ x ∣ ∣_{2} = \sum_{i = 1}^{n} ∣ x_{i} ∣^{2}$ ， $∣∣ x ∣ ∣_{\infty} = 1 \leq i \leq n max ∣ x_{i} ∣$ ，其中 $p \to \infty lim ∣∣ x ∣ ∣_{p} = ∣∣ x ∣ ∣_{\infty}$ 。

定义：向量序列 ${x^{(k)}}$ 收敛于向量 $x^{*}$ 是指对每一个 $1 \leq i \leq n$ 都有 $k \to \infty lim x_{i}^{(k)} = x_{i}^{*}$ ，也可以理解为 $∣∣ x^{(k)} - x^{*} ∣ ∣_{\infty} \to 0$ 。

定义：若存在常数 $C > 0$ 使得对任意 $x \in R^{n}$ 有 $∣∣ x ∣ ∣_{A} \leq C ∣∣ x ∣ ∣_{B}$ ，则称范数 $∣∣ \cdot ∣ ∣_{A}$ 比范数 $∣∣ \cdot ∣ ∣_{B}$ 强。

定义：若范数 $∣∣ \cdot ∣ ∣_{A}$ 比范数 $∣∣ \cdot ∣ ∣_{B}$ 强，同时 $∣∣ \cdot ∣ ∣_{B}$ 也比 $∣∣ \cdot ∣ ∣_{A}$ 强，即存在常数 $C_{1}, C_{2} > 0$ 使得 $C_{1} ∣∣ x ∣ ∣_{B} \leq ∣∣ x ∣ ∣_{A} \leq C_{2} ∣∣ x ∣ ∣_{B}$ ，则称 $∣∣ \cdot ∣ ∣_{A}$ 与 $∣∣ \cdot ∣ ∣_{B}$ 等价。

定理： $R^{n}$ 上一切范数都等价。

定义： $R^{m \times n}$ 空间的矩阵范数 $∣∣ \cdot ∣∣$ 对任意 $A, B \in R^{m \times n}$ 满足：（1）正定性： $∣∣ A ∣∣ \geq 0$ ， $∣∣ A ∣∣ = 0 \Leftrightarrow A = 0$ ；（2）齐次性： $∣∣ α A ∣∣ = ∣ α ∣ \cdot ∣∣ A ∣∣$ 对任意 $α \in C$ 成立；（3）三角不等式： $∣∣ A + B ∣∣ \leq ∣∣ A ∣∣ + ∣∣ B ∣∣$ ；（4）相容性( $m = n$ 时)： $∣∣ A B ∣∣ \leq ∣∣ A ∣∣ \cdot ∣∣ B ∣∣$ 。

要说明的是，不是所有矩阵范数都满足相容性，因为可能部分教材没有此条件，但是仍用前面三个条件定义了一些矩阵范数。

常用的矩阵范数有：

Frobenius范数（F-范数）： $∣∣ A ∣ ∣_{F} = \sum_{i = 1}^{m} \sum_{j = 1}^{n} ∣ a_{ij} ∣^{2}$ ，这是向量范数 $∣∣ \cdot ∣ ∣_{2}$ 的直接推广；如果 $A \in R^{n \times n}$ 且 $x \in R^{n}$ ，则有 $∣∣ A x ∣ ∣_{2} \leq ∣∣ A ∣ ∣_{F} \cdot ∣∣ x ∣ ∣_{2}$ （可以用Cauchy不等式证明）。

算子范数：如果 $A$ 是方阵，则可以由向量范数 $∣∣ \cdot ∣ ∣_{p}$ 导出关于矩阵 $A \in R^{n \times n}$ 的p-范数， $∣∣ A ∣ ∣_{p} = x \neq = 0 max \frac{∣∣ A x ∣ ∣ _{p}}{∣∣ x ∣ ∣ _{p}} = ∣∣ x ∣ ∣_{p} = 1 max ∣∣ A x ∣ ∣_{p}$ ；如果 $A$ 、 $B$ 是同阶方阵，则有 $∣∣ A B ∣ ∣_{p} \leq ∣∣ A ∣ ∣_{p} ∣∣ B ∣ ∣_{p}$ ， $∣∣ A x ∣ ∣_{p} \leq ∣∣ A ∣ ∣_{p} ∣∣ x ∣ ∣_{p}$ ，特别地，当p为 $\infty$ 、1、2时有：（1）行和范数： $∣∣ A ∣ ∣_{\infty} = 1 \leq i \leq n max \sum_{j = 1}^{n} ∣ a_{ij} ∣$ ；（2）列和范数： $∣∣ A ∣ ∣_{1} = 1 \leq j \leq n max \sum_{i = 1}^{n} ∣ a_{ij} ∣$ ；（3）谱范数： $∣∣ A ∣ ∣_{2} = λ_{m a x} (A^{T} A)$ （其中 $λ_{m a x} (A^{T} A)$ 指 $A^{T} A$ 之后的最大的那个特征根）。

注意：Frobenius范数不是算子范数，我们这里只关心有相容性的范数，而算子范数总是相容的。即使 $A$ 中元素全为实数，其特征根和相应的特征向量仍可能是复数，将上述定义中绝对值换成复数模仍均成立。

定义：矩阵 $A$ 的谱半径记为 $ρ (A) = 1 \leq i \leq n max ∣ λ_{i} ∣$ ，其中 $λ_{i}$ 为 $A$ 的特征根。

定理：对任意算子范数 $∣∣ \cdot ∣∣$ 有 $ρ (A) \leq ∣∣ A ∣∣$ 。

证明：由算子范数的相容性可得 $∣∣ A x ∣∣ \leq ∣∣ A ∣∣ \cdot ∣∣ x ∣∣$ ，将任意一个特征根 $λ$ 对应的特征向量 $u$ 代入得 $∣ λ ∣ \cdot ∣∣ u ∣∣ = ∣∣ λ u ∣∣ = ∣∣ A u ∣∣ \leq ∣∣ A ∣∣ \cdot ∣∣ u ∣∣$ 。

误差分析：方程组 ${12 x_{1} + 35 x_{2} = 59 12 x_{1} + 35.000001 x_{2} = 59.000001$ 的解为 ${x_{1} = 2 x_{2} = 1$ ，由于某种原因，第二个方程组的系数有一个小小的扰动（误差），变为 ${12 x_{1} + 35 x_{2} = 59 12 x_{1} + 34.999999 x_{2} = 59.000002$ ，此时方程组的解为 ${x_{1} = 10.75 x_{2} = - 2$ ，可见很小的误差也可能会对结果产生较大的影响，于是我们需要研究求解 $A x = b$ 时， $A$ 和 $b$ 的误差对解 $x$ 有何影响。

情况1：设 $A$ 精确， $b$ 有误差变为 $δ b$ ，则得到的解为 $x + δ x$ ，即 $A (x + δ x) = b + δ b$ ，消去 $A x = b$ 得 $δ x = A^{- 1} δ b$ ，根据相容性可知 $∣∣ δ x ∣∣ = ∣∣ A^{- 1} δ b ∣∣ \leq ∣∣ A^{- 1} ∣∣ \cdot ∣∣ δ b ∣∣$ ，因此将 $∣∣ A^{- 1} ∣∣$ 称为绝对误差放大因子；再看 $∣∣ b ∣∣ = ∣∣ A x ∣∣ \leq ∣∣ A ∣∣ \cdot ∣∣ x ∣∣ \Rightarrow \frac{1}{∣∣ x ∣∣} \leq \frac{∣∣ A ∣∣}{∣∣ b ∣∣}$ ，因此相对误差 $\frac{∣∣ δ x ∣∣}{∣∣ x ∣∣} \leq ∣∣ A ∣∣ \cdot ∣∣ A^{- 1} ∣∣ \cdot \frac{∣∣ δ b ∣∣}{∣∣ b ∣∣}$ ，因此将 $∣∣ A ∣∣ \cdot ∣∣ A^{- 1} ∣∣$ 称为相对误差放大因子。

情况2：设 $b$ 精确， $A$ 有误差 $δ A$ ，得到的解为 $x + δ x$ ，即 $(A + δ A) (x + δ x) = b \Rightarrow A (x + δ x) + δ A (x + δ x)$ ，消去 $A x = b$ 得 $A δ x + δ A (x + δ x) = 0 \Rightarrow δ x = - A^{- 1} δ A (x + δ x)$ ，根据相容性和三角不等式可知 $∣∣ δ x ∣∣ \leq ∣∣ A^{- 1} ∣∣ \cdot ∣∣ δ A ∣∣ \cdot ∣∣ x + δ x ∣∣ \leq ∣∣ A^{- 1} ∣∣ \cdot ∣∣ δ A ∣∣ \cdot ∣∣ x ∣∣ + ∣∣ A^{- 1} ∣∣ \cdot ∣∣ δ A ∣∣ \cdot ∣∣ δ x ∣∣$ ，于是得到 $(1 - ∣∣ A^{- 1} ∣∣ \cdot ∣∣ δ A ∣∣) \cdot ∣∣ δ x ∣∣ \leq ∣∣ A^{- 1} ∣∣ \cdot ∣∣ δ A ∣∣ \cdot ∣∣ x ∣∣$ ，故相对误差 $\frac{∣∣ δ x ∣∣}{∣∣ x ∣∣} \leq \frac{∣∣ A ^{- 1} ∣∣ \cdot ∣∣ δ A ∣∣}{1 - ∣∣ A ^{- 1} ∣∣ \cdot ∣∣ δ A ∣∣} = \frac{∣∣ A ∣∣ \cdot ∣∣ A ^{- 1} ∣∣ \cdot \frac{∣∣ δ A ∣∣}{∣∣ A ∣∣}}{1 - ∣∣ A ∣∣ \cdot ∣∣ A ^{- 1} ∣∣ \cdot \frac{∣∣ δ A ∣∣}{∣∣ A ∣∣}}$ ，即相对误差也与 $∣∣ A ∣∣ \cdot ∣∣ A^{- 1} ∣∣$ 有关。

Jacobi迭代法

Jacobi迭代法：对于方程组 $⎩ ⎨ ⎧ a_{11} x_{1} + a_{12} x_{2} + \dots + a_{1 n} x_{n} = b_{1} a_{21} x_{1} + a_{22} x_{2} + \dots + a_{2 n} x_{n} = b_{2} \dots a_{n 1} x_{1} + a_{n 2} x_{2} + \dots + a_{nn} x_{n} = b_{n}$ ，或写成矩阵的形式 $a_{11} a_{21} ⋮ a_{n 1} a_{12} a_{22} ⋮ a_{n 2} \dots \dots ⋱ \dots a_{1 n} a_{2 n} ⋮ a_{nn} x_{1} x_{2} ⋮ x_{n} = b_{1} b_{2} ⋮ b_{n}$ ，若 $a_{ii} \neq = 0$ ，则有 $⎩ ⎨ ⎧ x_{1} = \frac{1}{a _{11}} (- a_{12} x_{2} - \dots - a_{1 n} x_{n} + b_{1}) x_{2} = \frac{1}{a _{22}} (- a_{21} x_{1} - \dots - a_{2 n} x_{n} + b_{2}) \dots x_{n} = \frac{1}{a _{nn}} (- a_{n 1} x_{1} - \dots - a_{n, n - 1} x_{n - 1} + b_{n})$ 或写成矩阵形式 $x_{1} x_{2} ⋮ x_{n} = - 0 \frac{a _{21}}{a _{22}} ⋮ \frac{a _{n 1}}{a _{nn}} \frac{a _{12}}{a _{11}} 0 ⋮ \frac{a _{n 2}}{a _{nn}} \dots \dots ⋱ \dots \frac{a _{1 n}}{a _{11}} \frac{a _{2 n}}{a _{22}} ⋮ 0 x_{1} x_{2} ⋮ x_{n} + \frac{b _{1}}{a _{11}} \frac{b _{2}}{a _{22}} ⋮ \frac{b _{n}}{a _{nn}}$ ，即迭代原理中的把 $A x = b$ 写成 $x = B x + f$ 的形式，则Jacobi迭代格式为 $⎩ ⎨ ⎧ x_{1}^{(k + 1)} = \frac{1}{a _{11}} (- a_{12} x_{2}^{(k)} - \dots - a_{1 n} x_{n}^{(k)} + b_{1}) x_{2}^{(k + 1)} = \frac{1}{a _{22}} (- a_{21} x_{1}^{(k)} - \dots - a_{2 n} x_{n}^{(k)} + b_{2}) \dots x_{n}^{(k + 1)} = \frac{1}{a _{nn}} (- a_{n 1} x_{1}^{(k)} - \dots - a_{n, n - 1} x_{n - 1}^{(k)} + b_{n})$ （矩阵形式略），可以简写为 $x^{(k + 1)} = B x^{(k)} + f$ 。

对于 $A x = b$ ，分析矩阵形式：如果将 $A$ 分为 $A = L D U$ ，即右上角元素组成矩阵 $U$ 、对角线元素组成矩阵 $D$ 、左下角元素组成矩阵 $L$ ，则 $A x = b \Rightarrow (D + L + U) x = b \Rightarrow x = - D^{- 1} (L + U) x + D^{- 1} b = B x + f$ ，我们把 $- D^{- 1} (L + U)$ 称为Jacobi迭代阵，于是有迭代格式： $x^{(k + 1)} = - D^{- 1} (L + U) x^{(k)} + D^{- 1} b$ 。

例：用Jacobi方法解方程组 $⎩ ⎨ ⎧ 2 x_{1} - x_{2} - x_{3} = - 5 x_{1} + 5 x_{2} - x_{3} = 8 x_{1} + x_{2} + 10 x_{3} = 11$ .（已知精确解为 $x^{*} = (- 1, 2, 1)^{T}$ ）

解：将方程组写成等价形式 $⎩ ⎨ ⎧ x_{1} = 0.5 x_{2} + 0.5 x_{3} - 2.5 x_{2} = - 0.2 x_{1} + 0.2 x_{3} + 1.6 x_{3} = - 0.1 x_{1} - 0.1 x_{2} + 1.1$ ，或写成矩阵形式 $x_{1} x_{2} x_{3} = 0 - 0.2 - 0.1 0.5 0 - 0.1 0.5 0.2 0 x_{1} x_{2} x_{3} + - 2.5 1.6 1.1$ ，于是迭代格式为 $⎩ ⎨ ⎧ x_{1}^{(k + 1)} = 0.5 x_{2}^{(k)} + 0.5 x_{3}^{(k)} - 2.5 x_{2}^{(k + 1)} = - 0.2 x_{1}^{(k)} + 0.2 x_{3}^{(k)} + 1.6 x_{3}^{(k + 1)} = - 0.1 x_{1}^{(k)} - 0.1 x_{2}^{(k)} + 1.1$ ，取初始值 $x^{(0)} = (1, 1, 1)^{T}$ ，计算结果如下表：

$k$ $x_{1}^{(k)}$ $x_{2}^{(k)}$ $x_{3}^{(k)}$ $∣∣ x^{(k)} - x^{(k - 1)} ∣ ∣_{\infty}$

0 1 1 1

1 -1.5 1.6 0.9 2.5

2 -1.25 2.08 1.09 0.48

3 -0.915 2.068 1.017 0.335

4 -0.9575 1.9864 0.9847 0.0816

5 -1.01445 1.9844 0.99711 0.05695

6 -1.00722 2.00231 1.0026 0.01387

7 -0.997543 2.00197 1.00049 0.009687

$k$	$x_{1}^{(k)}$	$x_{2}^{(k)}$	$x_{3}^{(k)}$	$∣∣ x^{(k)} - x^{(k - 1)} ∣ ∣_{\infty}$
0	1	1	1
1	-1.5	1.6	0.9	2.5
2	-1.25	2.08	1.09	0.48
3	-0.915	2.068	1.017	0.335
4	-0.9575	1.9864	0.9847	0.0816
5	-1.01445	1.9844	0.99711	0.05695
6	-1.00722	2.00231	1.0026	0.01387
7	-0.997543	2.00197	1.00049	0.009687

Gauss-Seidel迭代法

引入：在Jacobi迭代格式中，如果算到 $x_{i}^{(k + 1)}$ ，则 $x_{i}^{(k + 1)}$ 完全是根据 $x_{1}^{(k)}, \dots, x_{n}^{(k)}$ 计算出来的，而前面计算出来的 $x_{1}^{(k + 1)}, \dots, x_{i - 1}^{(k + 1)}$ 并没有被利用到，考虑如果把新算出来的 $x_{1}^{(k + 1)}, \dots, x_{i - 1}^{(k + 1)}$ 代进 $x_{i}^{(k + 1)}$ 的公式中，会不会使其收敛得更快。

Gauss-Seidel迭代格式： $⎩ ⎨ ⎧ x_{1}^{(k + 1)} = \frac{1}{a _{11}} (- a_{12} x_{2}^{(k)} - a_{13} x_{3}^{(k)} - a_{14} x_{4}^{(k)} - \dots - a_{1 n} x_{n}^{(k)} + b_{1}) x_{2}^{(k + 1)} = \frac{1}{a _{22}} (- a_{21} x_{1}^{(k + 1)} - a_{23} x_{3}^{(k)} - a_{24} x_{4}^{(k)} - \dots - a_{2 n} x_{n}^{(k)} + b_{2}) x_{3}^{(k + 1)} = \frac{1}{a _{33}} (- a_{31} x_{1}^{(k + 1)} - a_{32} x_{3}^{(k + 1)} - a_{34} x_{4}^{(k)} - \dots - a_{3 n} x_{n}^{(k)} + b_{3}) \dots x_{n}^{(k + 1)} = \frac{1}{a _{nn}} (- a_{n 1} x_{1}^{(k + 1)} - a_{n 2} x^{(k + 1)} - a_{n 3} x_{3}^{(k + 1)} - \dots - a_{n, n - 1} x_{n - 1}^{(k + 1)} + b_{n})$ 。

与Jacobi迭代格式的矩阵表示类似，它可以写成 $x^{(k + 1)} = - D^{- 1} (L x^{(k + 1)} + U x^{(k)}) + D^{- 1} b$ ，进一步可以写成 $(D + L) x^{(k + 1)} = - U x^{(k)} + b \Rightarrow x^{(k + 1)} = - (D + L)^{- 1} U x^{(k)} + (D + L)^{- 1} b = B x + f$ ，我们把 $- (D + L)^{- 1} U$ 称为Gauss-Seidel迭代阵。

例：用Gauss-Seidel法解方程组 $⎩ ⎨ ⎧ 2 x_{1} - x_{2} - x_{3} = - 5 x_{1} + 5 x_{2} - x_{3} = 8 x_{1} + x_{2} + 10 x_{3} = 11$ .

解：迭代格式为 $⎩ ⎨ ⎧ x_{1}^{(k + 1)} = 0.5 x_{2}^{(k)} + 0.5 x_{3}^{(k)} - 2.5 x_{2}^{(k + 1)} = - 0.2 x_{1}^{(k + 1)} + 0.2 x_{3}^{(k)} + 1.6 x_{3}^{(k + 1)} = - 0.1 x_{1}^{(k + 1)} - 0.1 x_{2}^{(k + 1)} + 1.1$ ，取初始值 $x^{(0)} = (1, 1, 1)^{T}$ ，计算结果如下表：

$k$ $x_{1}$ $x_{2}$ $x_{3}$ $∣∣ x^{(k)} - x^{(k - 1)} ∣∣$

0 1 1 1

1 -1.5 2.1 1.04 2.5

2 -0.93 1.994 0.9936 0.57

3 -1.0062 1.99996 1.000624 0.0762

4 -0.999708 2.000064 0.99996416 0.006492

本例题的方程组和Jacobi法中的方程组完全相同，对比可以发现，Gauss-Seidel法确实收敛得更快。

$k$	$x_{1}$	$x_{2}$	$x_{3}$	$∣∣ x^{(k)} - x^{(k - 1)} ∣∣$
0	1	1	1
1	-1.5	2.1	1.04	2.5
2	-0.93	1.994	0.9936	0.57
3	-1.0062	1.99996	1.000624	0.0762
4	-0.999708	2.000064	0.99996416	0.006492

在某些情况下，Jacobi法和Gauss-Seidel法可能其中一种收敛，另一种不收敛，甚至两种都不收敛。比如：

例：分别用Jacobi法和Gauss-Seidel法解方程组 ${x_{1} + 12 x_{2} = 25 15 x_{1} + x_{2} = 17$ .

解：Jacobi法：取初始值 $x^{(0)} = (0, 0)^{T}$ ，得 $x^{(1)} = (17, 20)^{T}$ 、 $x^{(2)} = (- 223, - 235)^{T}$ ，显然是发散的；

Gauss-Seidel法：取初始值 $x^{(0)} = (0, 0)^{T}$ ，得 $x^{(1)} = (17, - 235)^{T}$ 、 $x^{(2)} = (2827, - 42535)^{T}$ ，显然也是发散的。

但是如果交换两个方程的次序后再构造相应的迭代公式，则两种方法都收敛，为此有必要了解迭代法的收敛条件。

迭代法的收敛性

无论是何种迭代方法，迭代格式都可以写成 $x^{(k + 1)} = B x^{(k)} + f$ 的形式，所以第 $k + 1$ 次的迭代误差可以写成 $e^{(k + 1)} = x^{(k + 1)} - x^{*} = (B x^{(k)} + f) - (B x^{*} + f) = B (x^{(k)} - x^{*}) = B e^{(k)}$ ，于是可以得到第k次迭代的误差为 $e^{(k)} = B^{k} e^{(0)}$ ，因此收敛的充要条件为：当 $k \to \infty$ 时 $e^{(k)} \to 0 \Leftrightarrow B^{k} \to 0$ 。

定义（矩阵序列的极限）：设 $A = (a_{ij})_{n \times n}, A_{k} = (a_{ij}^{(k)})_{n \times n} \in R^{n \times n}$ ， $k \to \infty lim A_{k} = A$ 是指 $k \to \infty lim a_{ij}^{(k)} = a_{ij}$ 对所有 $1 \leq i, j \leq n$ 成立，也可以理解为当 $k \to \infty$ 时 $∣∣ A_{k} - A ∣∣ \to 0$ 。

定理：设 $x = B x + f$ 存在唯一解，则从任意 $x^{(0)}$ 出发，迭代 $x^{(k + 1)} = B x^{(k)} + f$ 收敛 $\Leftrightarrow B^{k} \to 0$ 。

定理： $B^{(k)} \to 0 \Leftrightarrow ρ (B) < 1$ 。

证明：由前面的推导知， $∣∣ e^{(k)} ∣∣ \leq ∣∣ B ∣∣ \cdot ∣∣ e^{(k - 1)} ∣∣ \leq \dots \leq ∣∣ B ∣ ∣^{k} \cdot ∣∣ e^{(0)} ∣∣$ ，若 $∣∣ B ∣∣ < 1$ ，则当 $k \to \infty$ 时 $∣∣ B ∣ ∣^{k} \to 0$ ，即 $∣∣ e^{(k)} ∣∣ \to 0$ 。

定理（充分条件）：若存在一个矩阵范数使得 $∣∣ B ∣∣ = q < 1$ ，则迭代收敛，且有下列误差估计：（1） $∣∣ x^{*} - x^{(k)} ∣∣ \leq \frac{q}{1 - q} ∣∣ x^{(k)} - x^{(k - 1)} ∣∣$ ；（2） $∣∣ x^{*} - x^{(k)} ∣∣ \leq \frac{q ^{k}}{1 - q} ∣∣ x^{(1)} - x^{(0)} ∣∣$ 。

根据该定理，若事先给出误差的控制精度，可以估计出迭代次数，即 $∣∣ x^{*} - x^{(k)} ∣∣ \leq \frac{q ^{k}}{1 - q} ∣∣ x^{(1)} - x^{(0)} ∣∣ \leq ε$ ，解出 $k \geq \frac{l n \frac{ε ( 1 - ∣∣ B ∣∣ )}{∣∣ x ^{(1)} - x ^{(0)} ∣∣}}{l n ∣∣ B ∣∣}$ .

例：若用Jacobi法解方程组 $⎩ ⎨ ⎧ 20 x_{1} + 2 x_{2} + 3 x_{3} = 24 x_{1} + 8 x_{2} + x_{3} = 12 2 x_{1} - 3 x_{2} + 15 x_{3} = 30$ ，取 $x^{(0)} = (0, 0, 0)^{T}$ ，问Jacobi法是否收敛，若收敛，需要迭代多少次才能保证各分量误差绝对值小于 $1 0^{- 6}$ ？

解：Jacobi迭代格式为 $x_{1}^{(k + 1)} x_{2}^{(k + 1)} x_{3}^{(k + 1)} = 0 - \frac{1}{8} - \frac{2}{15} - \frac{2}{20} 0 \frac{3}{15} - \frac{3}{20} - \frac{1}{8} 0 x_{1}^{(k)} x_{2}^{(k)} x_{3}^{(k)} + \frac{24}{20} \frac{12}{8} \frac{30}{15} = B x + f$ ，因为 $∣∣ B ∣ ∣_{\infty} = \frac{1}{3} < 1$ ，所以迭代法收敛，迭代一次得到 $x_{1}^{(1)} = \frac{6}{5}$ 、 $x_{2}^{(1)} = \frac{3}{2}$ 、 $x_{3}^{(1)} = 2$ ， $∣∣ x^{(1)} - x^{(0)} ∣∣ = 2$ ，由 $\frac{q ^{k}}{1 - q} ∣∣ x^{(1)} - x^{(0)} ∣∣ \leq 1 0^{- 6}$ ，得到 $k > \frac{l n \frac{1 0 ^{- 6} ( 1 - \frac{1}{2} )}{2}}{l n \frac{1}{3}} \approx 13$ ，因此需要迭代14次。

前面定理中的 $∣∣ B ∣∣ < 1$ 只是迭代序列收敛的充分条件，即使 $∣∣ B ∣∣ < 1$ 对任何范数都不成立，迭代序列仍可能收敛。

例：设 $x^{(k + 1)} = B x^{(k)} + f$ ，其中 $B = [0.9 0.3 0 0.8]$ ， $f = [12]$ ，讨论迭代序列 ${x^{(k)}}$ 的收敛性.

解：显然 $∣∣ B ∣ ∣_{\infty} = 1.1$ 、 $∣∣ B ∣ ∣_{1} = 1.2$ 、 $∣∣ B ∣ ∣_{2} = 1.043$ 、 $∣∣ B ∣ ∣_{F} = 1.54$ ， $B$ 的范数均大于1，但是由于 $ρ (B) = 0.9 < 1$ ，故迭代序列 ${x^{(k)}}$ 仍是收敛的。

定理（充分条件）：若 $A$ 为严格对角优势阵(strictly diagonally dominant matrix)，则解 $A x = b$ 的Jacobi迭代和Gauss-Seidel迭代均收敛。

定义：如果 $A \in R^{n \times n}$ 满足 $∣ a_{ii} ∣ > \sum_{j = 1, j \neq = i}^{n} ∣ a_{ij} ∣$ （ $i = 1, 2, \dots, n$ ），则称 $A$ 为严格对角优势阵。

定义：如果 $A \in R^{n \times n}$ 满足 $∣ a_{ii} ∣ \geq \sum_{j = 1, j \neq = i}^{n} ∣ a_{ij} ∣$ （ $i = 1, 2, \dots, n$ ），且其中至少有一个等号严格成立，则称 $A$ 为弱对角优势阵。

比如 $⎩ ⎨ ⎧ 5 x_{1} + x_{2} + 2 x_{3} = - 3 2 x_{1} - 8 x_{2} + 4 x_{3} = 2 x_{1} + 3 x_{2} + 12 x_{3} = 11$ ， $A = 521 1 - 8 3 2412$ ，显然 $A$ 为严格对角优势阵，因此Jacobi迭代和Gauss-Seidel迭代均收敛。

四、插值与拟合

1. 插值多项式的唯一性

有时我们会遇到这样的问题：由实验或测量得到一批离散样点，要求作出一条通过这些点的光滑曲线，以便满足设计要求或进行加工，此外，还有一些函数虽有表达式，但因式子复杂，不易计算其值和进行理论分析。

当精确函数 $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_{i}) = f (x_{i})$ （ $i = 0, \dots, n$ ），这里的 $g (x)$ 就称为 $f (x)$ 的插值函数。我们常用代数多项式作为插值函数。

插值多项式是指根据给定的 $n + 1$ 个点 $(x_{i}, y_{i})$ （ $i = 0, 1, \dots, n$ ），求一个 $n$ 次多项式 $p (x) = a_{0} + a_{1} x + \dots + a_{n} x^{n}$ 。

定理（插值多项式唯一性）：满足 $P (x_{i}) = y_{i}$ （ $i = 0, 1, \dots, n$ ）的 $n$ 阶插值多项式是唯一存在的，其中 $x_{0}, \dots, x_{n}$ 为互不相同的节点。

证明：（利用范德蒙(Vandermonde)行列式论证）

对于 $⎩ ⎨ ⎧ p (x_{0}) = a_{0} + a_{1} x_{0} + \dots + a_{n} x_{0}^{n} p (x_{1}) = a_{0} + a_{1} x_{1} + \dots + a_{n} x_{1}^{n} \dots p (x_{n}) = a_{0} + a_{1} x_{n} + \dots + a_{n} x_{n}^{n}$ ，其中 $a_{0}, \dots, a_{n}$ 是 $n + 1$ 个待定系数，其系数行列式为 $11 ⋮ 1 x_{0} x_{1} ⋮ x_{n} \dots \dots ⋱ \dots x_{0}^{n} x_{1}^{n} ⋮ x_{n}^{n} \neq = 0$ （由于 $x_{0}, \dots, x_{n}$ 互异，根据范德蒙行列式的计算公式 $\prod_{1 \leq j \leq i \leq n} (x_{i} - x_{j})$ ，必不等于零），由于系数行列式不等于零，故解唯一存在。

2. 拉格朗日插值

求 $n$ 次多项式 $P_{n} (x) = a_{0} + a_{1} x + \dots + a_{n} x^{n}$ 使得 $P_{n} (x_{i}) = y_{i}$ （ $i = 0, 1, \dots, n$ ），其中当 $i \neq = j$ 时 $x_{i} \neq = x_{j}$ ，即无重合节点。

拉格朗日线性插值：当 $n = 1$ 时， $P_{1} (x) = a_{0} + a_{1} x$ ， $P_{1} (x_{0}) = y_{0}$ 、 $P_{1} (x_{1}) = y_{1}$ ，可见 $P_{1}$ 是过点 $(x_{0}, y_{0})$ 和 $(x_{1}, y_{1})$ 两点的直线，根据直线方程， $P_{1} (x) = y_{0} + \frac{y _{1} - y _{0}}{x _{1} - x _{0}} (x - x_{0})$ ，写成一般形式 $P_{1} (x) = (\frac{x - x _{1}}{x _{0} - x _{1}}) y_{0} + (\frac{x - x _{0}}{x _{1} - x _{0}}) y_{1} = l_{0} (x) y_{0} + l_{1} (x) y_{1} = \sum_{i = 0}^{1} l_{i} (x) y_{i}$ ，我们称 $l_{i} (x)$ 为拉格朗日基函数（Lagrange Basis），满足条件 $l_{i} (x_{j}) = δ_{ij}$ （其中 $δ_{ij}$ 是克罗内克函数，Kronecker Delta，当 $i = j$ 时函数值为1，否则为0），称 $P_{1} (x)$ 为一次插值多项式或线性插值多项式（显然 $l_{0} (x)$ 和 $l_{1} (x)$ 也是一次插值多项式）。

可见插值基函数与 $y_{0}$ 、 $y_{1}$ 无关，而由插值点 $x_{0}, x_{1}$ 所决定。

拉格朗日抛物线插值：当 $n = 2$ 时，设插值多项式为 $P_{2} (x) = l_{0} (x) y_{0} + l_{1} (x) y_{1} + l_{2} (x) y_{2}$ ，类似 $n = 1$ 时的情况，要求：①插值基函数 $l_{0}, l_{1}, l_{2}$ 为二次多项式；②它们的函数值满足下表：

	$x_{0}$	$x_{1}$	$x_{2}$
$l_{0} (x)$	1	0	0
$l_{1} (x)$	0	1	0
$l_{2} (x)$	0	0	1

可见 $l_{0} (x)$ 只在 $x_{0}$ 处为1，在 $x_{1}, x_{2}$ 处为0，即 $x_{1}, x_{2}$ 是 $l_{0} (x)$ 的零点，故 $l_{0} (x) = a (x - x_{1}) (x - x_{2})$ ，再根据 $l_{0} (x_{0}) = 1$ 解出 $a = \frac{1}{( x _{0} - x _{1} ) ( x _{0} - x _{2} )}$ ，即 $l_{0} (x) = \frac{( x - x _{1} ) ( x - x _{2} )}{( x _{0} - x _{1} ) ( x _{0} - x _{2} )}$ ；同理可得 $l_{1} (x) = \frac{( x - x _{0} ) ( x - x _{2} )}{( x _{1} - x _{0} ) ( x _{1} - x _{2} )}$ 、 $l_{2} (x) = \frac{( x - x _{0} ) ( x - x _{1} )}{( x _{2} - x _{0} ) ( x _{2} - x _{1} )}$ 。

拉格朗日n次插值：当 $n > 2$ 时，假定 $P_{n} (x_{i}) = y_{i}$ （ $i = 0, 1, \dots, n$ ），设 $P_{n} (x) = \sum_{i = 0}^{n} l_{i} (x) y_{i}$ ，其中 $l_{i} (x_{j}) = δ_{ij}$ ，和 $n = 2$ 的情况类似，推导出 $l_{i} (x) = \frac{( x - x _{0} ) \dots ( x - x _{i - 1} ) ( x - x _{i + 1} ) \dots ( x - x _{n} )}{( x _{i} - x _{0} ) \dots ( x _{i} - x _{i - 1} ) ( x _{i} - x _{i + 1} ) \dots ( x _{i} - x _{n} )} = \prod_{j = 0, j \neq = i}^{n} \frac{( x - x _{j} )}{( x _{i} - x _{j} )}$ 。

有时也将拉格朗日插值多项式 $P_{n} (x)$ 记作 $L_{n} (x)$ 。

例：已知 $l g 10 = 1$ 、 $l g 20 = 1.3010$ ，利用一次插值多项式求 $l g 12$ 的近似值.

解：设 $f (x) = l g x$ ，则 $f (10) = 1$ 、 $f (20) = 1.3010$ ，则 $l_{0} (x) = \frac{x - x _{1}}{x _{0} - x _{1}} = \frac{x - 20}{10 - 20} = - \frac{1}{10} (x - 20)$ ， $l_{1} (x) = \frac{x - x _{0}}{x _{1} - x _{0}} = \frac{1}{10} (x - 10)$ ，故 $L_{1} (x) = l_{0} (x) y_{0} + l_{1} (x) y_{1} = - \frac{1}{10} (x - 20) y_{0} + \frac{1}{10} (x - 10) y_{1} = - \frac{1}{10} (x - 20) + \frac{1.3010}{10} (x - 10)$ ，则 $L_{1} (12) = 1.0602$ 。

例：若已知

$x_{i}$ 10 15 20

$y_{i} = l g x_{i}$ 1 1.1761 1.3010

用二次插值多项式求 $l g 12$ 近似值.

解：设 $x_{0} = 10$ 、 $x_{1} = 15$ 、 $x_{2} = 20$ ，根据公式得到 $l_{0} (x) = \frac{( x - 15 ) ( x - 20 )}{( 10 - 15 ) ( 10 - 20 )} = \frac{1}{50} (x - 15) (x - 20)$ ， $l_{1} (x) = \frac{( x - 10 ) ( x - 20 )}{( 15 - 10 ) ( 15 - 20 )} = - \frac{1}{25} (x - 10) (x - 20)$ ， $l_{2} (x) = \frac{( x - 10 ) ( x - 15 )}{( 20 - 10 ) ( 20 - 15 )} = \frac{1}{50} (x - 10) (x - 15)$ ，故 $L_{2} (x) = y_{0} l_{0} (x) + y_{1} l_{1} (x) + y_{2} l_{2} (x) = \frac{1}{50} (x - 15) (x - 20) - \frac{1.1761}{25} (x - 10) (x - 20) + \frac{1.3010}{50} (x - 10) (x - 15)$ ，所以 $L_{2} (12) = 1.0766$ 。

$l g 12$ 的精确值为 $1.079181246$ ，对比两个例题，发现高次插值的误差要比低次插值的误差小，但这不是绝对的。

$x_{i}$	10	15	20
$y_{i} = l g x_{i}$	1	1.1761	1.3010

定理（拉格朗日插值余项）：设节点 $a \leq x_{0} < x_{1} < \dots < x_{n} \leq b$ ，且 $f$ 满足条件 $f \in C^{n} [a, b]$ ， $f^{(n + 1)}$ 在 $[a, b]$ 存在，则 $\forall x \in [a, b]$ ，其截断误差为 $R_{n} (x) = f (x) - P_{n} (x) = \frac{f ^{(n + 1)} ( ξ _{x} )}{( n + 1 )!} \prod_{i = 0}^{n} (x - x_{i})$ 。

证明：回顾罗尔定理：若 $φ (x)$ 充分光滑，且 $φ (x_{0}) = φ (x_{1}) = 0$ ，则存在 $ξ \in (x_{0}, x_{1})$ 使得 $φ^{'} (ξ) = 0$ .

推广1：若 $φ (x_{0}) = φ (x_{1}) = φ (x_{2}) = 0$ ，则存在 $ξ_{0} \in (x_{0}, x_{1})$ 、 $ξ_{1} \in (x_{1}, x_{2})$ 使得 $φ^{'} (ξ_{0}) = φ^{'} (ξ_{1}) = 0$ ，并且存在 $ξ \in (ξ_{0}, ξ_{1})$ 使得 $φ^{''} (ξ) = 0$ .

推广2：若 $φ (x_{0}) = \dots = φ (x_{n}) = 0$ ，则存在 $ξ \in (a, b)$ 使得 $φ^{(n)} (ξ) = 0$ 。

因为拉格朗日插值多项式 $L_{n} (x)$ 经过 $(x_{0}, y_{0}), (x_{1}, y_{1}), \dots, (x_{n}, y_{n})$ ，故在这些点处的误差 $R_{n} (x_{i}) = 0$ ，因此 $R_{n} (x)$ 至少有 $n + 1$ 个根，可设 $R_{n} (x) = K (x) \prod_{i = 0}^{n} (x - x_{i})$ ；任意固定的 $x \neq = x_{i}$ （ $i = 0, 1, \dots, n$ ），考察 $φ (t) = R_{n} (t) - K (x) \prod_{i = 0}^{n} (t - x_{i})$ ，显然 $φ (t)$ 有 $n + 2$ 个不同的根 $x_{0}, \dots, x_{n}, x$ ，根据前面罗尔定理的推广知，存在 $ξ_{x} \in (a, b)$ 使得 $φ^{(n + 1)} (ξ_{x}) = 0$ （其中 $φ^{(n + 1)}$ 是对 $t$ 求导），即 $φ^{(n + 1)} (ξ_{x}) = R_{n}^{(n + 1)} (ξ_{x}) - K (x) (n + 1)! = f^{(n + 1)} (ξ_{x}) - L_{n}^{(n + 1)} (ξ_{x}) - K (x) (n + 1)! = 0$ ，其中 $L_{n} (x)$ 是 $x$ 的 $n$ 次多项式，求 $n + 1$ 阶导后结果为零，故 $K (x) = \frac{f ^{(n + 1)} ( ξ _{x} )}{( n + 1 )!}$ ，将 $K (x)$ 代入 $R_{n} (x)$ 中得证。

注意：通常不能确定 $ξ_{x}$ ，而是对 $\forall x \in (a, b)$ 估计 $∣ f^{(n + 1)} (x) ∣ \leq M_{n + 1}$ ，得到 $n + 1$ 阶导的上限 $M_{n + 1}$ ，将 $\frac{M _{n + 1}}{( n + 1 )!} \prod_{i = 0}^{n} ∣ x - x_{i} ∣$ 作为误差估计上限。特别地，当 $f (x)$ 为任意一个次数 $\leq n$ 的多项式时， $f^{(n + 1)} (x) \equiv 0$ ，即插值多项式对于次数 $\leq n$ 的多项式是精确的。

3. 牛顿插值

Lagrange插值虽然格式整齐规范，但却没有承袭性质，若要增加一个节点，全部基函数 $l_{i} (x)$ 都要重新算过。如果可以把Lagrange插值改写为 $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})$ 的形式，这样每加一个节点，只需附加一项上去就可以了。

为介绍牛顿插值，需引入差商的概念。

定义：差商用归纳定义法来定义：1阶差商为 $f [x_{i}, x_{j}] = \frac{f ( x _{i} ) - f ( x _{j} )}{x _{i} - x _{j}}$ （ $i \neq = j$ 、 $x_{i} \neq = x_{j}$ ）；2阶差商为 $f [x_{i}, x_{j}, x_{k}] = \frac{f [ x _{i} , x _{j} ] - f [ x _{j} , x _{k} ]}{x _{i} - x _{k}}$ （ $i \neq = k$ ）； $\dots$ ； $k + 1$ 阶差商为 $f [x_{0}, \dots, x_{k + 1}] = \frac{f [ x _{0} , x _{1} , \dots , x _{k} ] - f [ x _{1} , \dots , x _{k} , x _{k + 1} ]}{x _{0} - x _{k + 1}} = \frac{f [ x _{0} , x _{1} , \dots , x _{k} ] - f [ x _{0} , \dots , x _{k - 1} , x _{k + 1} ]}{x _{k} - x _{k + 1}}$ 。

差商的性质：

差商与函数值的关系： $k$ 阶差商 $f [x_{0}, x_{1}, \dots, x_{k}]$ 是由函数值 $f (x_{0}), f (x_{1}), \dots, f (x_{k})$ 的线形组合而成，即 $f [x_{0}, x_{1}, \dots, x_{k}] = \sum_{i = 0}^{k} \frac{1}{( x _{i} - x _{0} ) \dots ( x _{i} - x _{i - 1} ) ( x _{i} - x _{i + 1} ) \dots ( x _{i} - x _{k} )} f (x_{i})$ ；
对称性：差商的值只与节点有关，而与 $x_{i}$ 的顺序无关，即，若 $i_{0}, i_{1}, \dots, i_{k}$ 为 $0, 1, \dots, k$ 的任一排列，则 $f [x_{0}, x_{1}, \dots, x_{k}] = f [x_{i_{0}}, x_{i_{1}}, \dots, x_{i_{k}}]$ ；

我们常把差商列成表，形式如下表（计算时先算一阶，再算二阶，依此类推）：

$i$	$x_{i}$	$f (x_{i})$	一阶差商	二阶差商	$\dots$	$n$ 阶差商
0	$x_{0}$	$f (x_{0})$
1	$x_{1}$	$f (x_{1})$	$f [x_{0}, x_{1}]$
2	$x_{2}$	$f (x_{2})$	$f [x_{1}, x_{2}]$	$f [x_{0}, x_{1}, x_{2}]$
3	$x_{3}$	$f (x_{3})$	$f [x_{2}, x_{3}]$	$f [x_{1}, x_{2}, x_{3}]$
$⋮$	$⋮$	$⋮$	$⋮$	$⋮$	$⋱$
$n$	$x_{n}$	$f (x_{n})$	$f [x_{n - 1}, x_{n}]$	$f [x_{n - 2}, x_{n - 1}, x_{n}]$	$\dots$	$f [x_{0}, \dots, x_{n}]$

例：计算 $(- 2, 17)$ 、 $(0, 1)$ 、 $(1, 2)$ 、 $(2, 19)$ 的一至三阶差商.

解：列出差商表如下

$i$ $x_{i}$ $f (x_{i})$ $f [x_{i - 1}, x_{i}]$ $f [x_{i - 2}, x_{i - 1}, x_{i}]$ $f [x_{i - 2}, x_{i - 1}, x_{i}]$

0 -2 17

1 0 1 $f [x_{0}, x_{1}] = f [- 2, 0] = - 8$

2 1 2 $f [x_{1}, x_{2}] = f [0, 1] = 1$ $f [x_{0}, x_{1}, x_{2}] = f [- 2, 0, 1] = 3$

3 2 19 $f [x_{2}, x_{3}] = f [1, 2] = 17$ $f [x_{1}, x_{2}, x_{3}] = f [0, 1, 2] = 8$ $f [x_{0}, x_{1}, x_{2}, x_{3}] = \frac{5}{4}$

比如 $f [x_{0}, x_{1}] = f [- 2, 0] = \frac{f ( x _{1} ) - f ( x _{0} )}{x _{1} - x _{0}} = \frac{1 - 17}{0 - ( - 2 )} = - 8$ ， $f [x_{1}, x_{2}] = f [0, 1] = \frac{f ( x _{2} ) - f ( x _{1} )}{x _{2} - x _{1}} = \frac{2 - 1}{1 - 0} = 1$ ， $f [x_{2}, x_{3}] = f [1, 2] = \frac{f ( x _{3} ) - f ( x _{2} )}{x _{3} - x _{2}} = \frac{19 - 2}{2 - 1} = 17$ ， $f [x_{0}, x_{1}, x_{2}] = f [- 2, 0, 1] = \frac{f [ x _{1} , x _{2} ] - f [ x _{0} , x _{1} ]}{x _{2} - x _{0}} = \frac{1 - ( - 8 )}{1 - ( - 2 )} = 3$ ， $f [x_{1}, x_{2}, x_{3}] = f [0, 1, 2] = \frac{f [ x _{2} , x _{3} ] - f [ x _{1} , x _{2} ]}{x _{3} - x _{1}} = \frac{17 - 1}{2 - 0} = 8$ ， $f [x_{0}, x_{1}, x_{2}, x_{3}] = f [- 2, 0, 1, 2] = \frac{f [ x _{1} , x _{2} , x _{3} ] - f [ x _{0} , x _{1} , x _{2} ]}{x _{3} - x _{0}} = \frac{8 - 3}{2 - ( - 2 )} = \frac{5}{4}$ 。

$i$	$x_{i}$	$f (x_{i})$	$f [x_{i - 1}, x_{i}]$	$f [x_{i - 2}, x_{i - 1}, x_{i}]$	$f [x_{i - 2}, x_{i - 1}, x_{i}]$
0	-2	17
1	0	1	$f [x_{0}, x_{1}] = f [- 2, 0] = - 8$
2	1	2	$f [x_{1}, x_{2}] = f [0, 1] = 1$	$f [x_{0}, x_{1}, x_{2}] = f [- 2, 0, 1] = 3$
3	2	19	$f [x_{2}, x_{3}] = f [1, 2] = 17$	$f [x_{1}, x_{2}, x_{3}] = f [0, 1, 2] = 8$	$f [x_{0}, x_{1}, x_{2}, x_{3}] = \frac{5}{4}$

牛顿插值多项式：形式为 $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})$ ，其中 $a_{0} = f (x_{0})$ 、 $a_{i} = f [x_{0}, \dots, x_{i}]$ 。

证明：根据差商公式有Ⅰ. $f (x) = f (x_{0}) + (x - x_{0}) f [x, x_{0}]$ 、Ⅱ. $f [x, x_{0}] = f [x_{0}, x_{1}] + (x - x_{1}) f [x, x_{0}, x_{1}]$ 、 $\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}]$ ，将 $式 1 + (x - x_{0}) \times 式 2 + \dots + (x - x_{0}) \dots (x - x_{n - 1}) \times 式 3$ ，整理后得到 $f (x) = 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 + f [x_{0}, \dots, x_{n}] (x - x_{0}) \dots (x - x_{n - 1}) + f [x, x_{0}, \dots, x_{n}] (x - x_{0}) \dots (x - x_{n - 1}) (x - x_{n})$ ，其中前面的 $n + 1$ 项是牛顿插值多项式，记作 $N_{n} (x)$ ，最后一项是牛顿插值的余项，记作 $R_{n} (x)$ 。

注意：由唯一性可知 $N_{n} (x) = L_{n} (x)$ ，只是计算方式不同，故其余项也相同，即 $f [x, x_{0}, \dots, x_{n}] ω_{n + 1} (x) = \frac{f ^{(n + 1)} ( ξ _{x} )}{( n + 1 )!} ω_{n + 1} (x)$ ，因此 $f [x, x_{0}, \dots, x_{n}] = \frac{f ^{(n + 1)} ( ξ _{x} )}{( n + 1 )!}$ 。

例：已知 $x = 1, 4, 9$ 的平方根为 $1, 2, 3$ ，利用牛顿差商公式求 $7$ 的近似值.

解：先根据给定的节点及函数值，构造出差商表

$x$ $y$ 一阶差商二阶差商

1 1

4 2 $\frac{2 - 1}{4 - 1} = \frac{1}{3}$

9 3 $\frac{3 - 2}{9 - 4} = \frac{1}{5}$ $\frac{\frac{1}{5} - \frac{1}{3}}{9 - 1} = - \frac{1}{60}$

则二次牛顿插值多项式 $N_{2} (x) = 1 = \frac{1}{3} (x - 1) + \frac{- 1}{60} (x - 1) (x - 4)$ ，因此 $7$ 的近似值为 $N_{2} (7) = 2.69992$ 。

$x$	$y$	一阶差商	二阶差商
1	1
4	2	$\frac{2 - 1}{4 - 1} = \frac{1}{3}$
9	3	$\frac{3 - 2}{9 - 4} = \frac{1}{5}$	$\frac{\frac{1}{5} - \frac{1}{3}}{9 - 1} = - \frac{1}{60}$

4. 分段低次插值

为了提高插值多项式对函数的逼近程度，自然希望增加节点个数，即提高插值多项式的次数，特别地，当 $n \to \infty$ 时，期望 $P_{n} (x) \to f (x)$ ，但是事实并非如此，比如：

例：在 $[- 5, 5]$ 上考察 $f (x) = \frac{1}{1 + x ^{2}}$ 的 $L_{n} (x)$ .

解：取 $x_{i} = - 5 + \frac{10}{n} i$ ，其中 $i = 0, \dots, n$ ，则 $L_{n} (x) = \sum_{i = 0}^{n} f (x_{i}) \cdot l_{i} (x) = \sum_{i = 0}^{n} f (x_{i}) \cdot \prod_{j = 0, j \neq = i}^{n} \frac{( x - x _{j} )}{( x _{i} - x _{j} )}$ ，我们考察最后两个节点的中点，记作 $x_{n - \frac{1}{2}} = \frac{1}{2} (x_{n - 1} + x_{n}) = 5 - \frac{5}{n}$ ，如图所示红、绿、橙分别是 $n = 2, 5, 10$ 的曲线：

当 $n = 5$ 时 $L_{5} (x)$ 的最后两个点的中点是最接近原函数的，当 $n = 10$ 的时候偏差反而比 $n = 5$ 的时候大了。

也就是说，当 $n \to \infty$ 时 $L_{n} (x)$ 并不能逼近 $f (x)$ ，我们把这种 $n$ 越大，端点附近抖动越大的现象称为Runge现象。

由于高次插值收敛性无法保证，因此当插值节点 $n$ 较大时，通常不采用高次多项式插值，而改用分段低次插值。

分段线性插值：在每个区间 $[x_{i}, x_{i + 1}]$ 上，用1阶多项式（直线）逼近 $f (x)$ ，即 $f (x) \approx P_{1} (x) = \frac{x - x _{i + 1}}{x _{i} - x _{i + 1}} y_{i} + \frac{x - x _{i}}{x _{i + 1} - x _{i}} y_{i + 1}$ ，这种插值方法可以逼近函数，但是失去了原函数的光滑性（左右导数不相等）。

例：已知函数 $y = f (x) = \frac{1}{1 + x ^{2}}$ 在区间 $[0, 5]$ 上取等距插值节点（如下表），求区间上分段线性插值函数，并用它求出 $f (4.5)$ 的近似值.

$x_{i}$ 0 1 2 3 4 5

$y_{i}$ 1 0.5 0.2 0.1 0.05882 0.03846

解：利用公式 $P (x) = \frac{x - ( k + 1 )}{k - ( k + 1 )} y_{k} + \frac{x - k}{( k + 1 ) - k} y_{k + 1} = - y_{k} (x - k - 1) + y_{k + 1} (x - k)$ ，得到 $⎩ ⎨ ⎧ - (x - 1) + 0.5 x - 0.5 (x - 2) + 0.2 (x - 1) - 0.2 (x - 3) + 0.1 (x - 2) - 0.1 (x - 4) + 0.05882 (x - 3) - 0.05882 (x - 5) + 0.03846 (x - 4), x \in [0, 1], x \in [1, 2], x \in [2, 3], x \in [3, 4], x \in [4, 5]$ ，于是 $P (4.5) = 0.04864$ 。

$x_{i}$	0	1	2	3	4	5
$y_{i}$	1	0.5	0.2	0.1	0.05882	0.03846

5. 埃尔米特(Hermite)插值

有时我们不仅要求函数值重合，而且要求若干阶导数也重合。

提示：根据解方程的思想，如果题目给了 $N$ 个条件，则可以构造 $N - 1$ 阶多项式。

下面通过例题了解埃尔米特插值多项式如何构造。

例：设 $x_{0} \neq = x_{1} \neq = x_{2}$ ，已知 $f (x_{0})$ 、 $f (x_{1})$ 、 $f (x_{2})$ 和 $f^{'} (x_{1})$ ，求多项式 $P (x)$ 满足 $P (x_{i}) = f (x_{i})$ ，且 $P^{'} (x_{1}) = f^{'} (x_{1})$ .

解：题目给出了 $4$ 个条件，故 $P (x)$ 的阶数为 $3$ ，模仿Lagrange多项式的思想，将 $P (x)$ 写成用插值基函数表示的形式， $P_{3} (x) = \sum_{i = 0}^{2} f (x_{i}) α_{i} (x) + f^{'} (x_{1}) β_{1} (x)$ ，其中 $α_{i} (x)$ 和 $β_{1} (x)$ 都是埃尔米特插值基函数，都是三次多项式，并且它们满足 $α_{i} (x_{j}) = δ_{ij}$ 、 $a_{i}^{'} (x_{1}) = 0$ 、 $β_{1} (x_{j}) = 0$ 、 $β_{1}^{'} (x_{1}) = 1$ ，这里的 $δ_{ij}$ 是克罗内克函数，列出函数值如下表：

$α_{0} (x)$ $α_{1} (x)$ $α_{2} (x)$ $β_{1} (x)$

函数值 $x_{0}$ 1 0 0 0

函数值 $x_{1}$ 0 1 0 0

函数值 $x_{2}$ 0 0 1 0

导数值 $x_{1}$ 0 0 0 1

比如 $α_{0} (x_{0}) = 1$ ，而 $α_{0} (x_{1}) = α_{0} (x_{2}) = 0$ ，因此 $α_{0} (x)$ 有 $x_{1}, x_{2}$ 两个根，而且 $a_{0}^{'} (x_{1}) = 0$ ，因此 $x_{1}$ 是重根，故可以设 $α_{0} (x) = C_{0} (x - x_{1})^{2} (x - x_{2})$ ，其中常数 $C_{0}$ 可以通过 $α_{0} (x) = 1$ 算出，最后得到 $α_{0} (x) = \frac{( x - x _{1} ) ^{2} ( x - x _{2} )}{( x _{0} - x _{1} ) ^{2} ( x _{0} - x _{2} )}$ ；（ $α_{2}$ 同理）

而 $α_{1} (x)$ 有 $x_{0}, x_{2}$ 两个根，但 $x_{1}$ 不是重根，因此可设 $α_{1} (x) = (A x + B) (x - x_{0}) (x - x_{2})$ ，根据 $α_{1} (x_{1}) = 1$ 和 $α_{1}^{'} (x_{1}) = 0$ 可解出 $A$ 和 $B$ ；

而 $β_{1} (x)$ 有 $x_{0}, x_{1}, x_{2}$ 三个根，因此可设 $β_{1} (x) = C_{1} (x - x_{0}) (x - x_{1}) (x - x_{2})$ ，其中常数 $C_{1}$ 可以根据 $β_{1} (x_{1}) = 1$ 解出。

		$α_{0} (x)$	$α_{1} (x)$	$α_{2} (x)$	$β_{1} (x)$
函数值	$x_{0}$	1	0	0	0
函数值	$x_{1}$	0	1	0	0
函数值	$x_{2}$	0	0	1	0
导数值	$x_{1}$	0	0	0	1

埃尔米特插值多项式的一般形式：已知 $x_{0}, \dots, x_{n}$ 处有 $y_{0}, \dots, y_{n}$ 和 $y_{0}^{'}, \dots, y_{n}^{'}$ ，求 $H_{2 n + 1} (x)$ 满足 $H_{2 n + 1} (x_{i}) = y_{i}$ ， $H_{2 n + 1}^{'} (x_{i}) = y_{i}^{'}$ ，则 $H_{2 n + 1} (x) = \sum_{i = 0}^{n} [1 - 2 l_{i}^{'} (x_{i}) (x - x_{i})] l_{i}^{2} (x) y_{i} + \sum_{i = 0}^{n} (x - x_{i}) l_{i}^{2} (x) y_{i}$ ，其中 $l_{i} (x) = \prod_{j \neq = i} \frac{( x - x _{j} )}{( x _{i} - x _{j} )}$ 。

推导（和例题类似）：设 $H_{2 n + 1} (x) = \sum_{i = 0}^{n} y_{i} α_{i} (x) + \sum_{i = 0}^{n} y_{i} β_{i} (x)$ ，其中 $α_{i} (x)$ 和 $β_{i} (x)$ 都是 $2 n + 1$ 次多项式，并且满足 $α_{i} (x_{j}) = δ_{ij}$ 、 $α_{i}^{'} (x_{j}) = 0$ 、 $β_{i} (x_{j}) = 0$ 、 $β_{i}^{'} (x_{j}) = δ_{ij}$ ；

根据这些条件可知，对于 $α_{i} (x)$ ，除 $x_{i}$ 外其余的 $x_{0}, \dots, x_{n}$ 都是二重根，因此 $α_{i} (x) = (A_{i} x + B_{i}) l_{i}^{2} (x)$ ，其中 $l_{i} (x) = \prod_{j \neq = i} \frac{( x - x _{j} )}{( x _{i} - x _{j} )}$ ，根据 $α_{i} (x_{i}) = 1$ 和 $α_{i}^{'} (x_{i}) = 0$ 可求出 $A_{i}$ 和 $B_{i}$ ，最后得到 $α_{i} (x) = [1 - 2 l_{i}^{'} (x_{i}) (x - x_{i})] l_{i}^{2} (x)$ ；

对于 $β_{i} (x)$ ， $x_{0}, \dots, x_{n}$ 都是它的根，而且除 $x_{i}$ 外其余的 $x_{0}, \dots, x_{n}$ 都是二重根，因此 $β_{i} (x) = C_{i} (x - x_{i}) l_{i}^{2} (x)$ ，其中 $l_{i} (x) = \prod_{j \neq = i} \frac{( x - x _{j} )}{( x _{i} - x _{j} )}$ ，根据 $β_{i}^{'} (x_{i}) = 1$ 可算出 $C_{i} = 1$ ，最后得到 $β_{i} (x) = (x - x_{i}) l_{i}^{2} (x)$ 。

（ $β_{i} (x)$ 图像特点：与 $x$ 轴交点为 $x_{0}, \dots, x_{n}$ ，且在 $x_{i}$ 处斜率为 $1$ ，其余 $x_{j}$ ( $j \neq = i$ )处斜率为 $0$ ）

埃尔米特余项：若 $f (x)$ 满足 $f^{(2 n + 1)} (x)$ 在 $[a, b]$ 上连续， $f^{(2 n + 2)} (x)$ 在 $(a, b)$ 上存在，又 $x_{0}, x_{1}, \dots, x_{n}$ 是 $n + 1$ 个互异节点，则满足插值条件 $H_{2 n + 1} (x_{i}) = y_{i}$ 、 $H_{2 n + 1}^{'} (x_{i}) = y_{i}^{'}$ 的埃尔米特插值多项式 $H_{2 n + 1} (x)$ 唯一存在，且插值余项为 $R_{n} (x) = \frac{f ^{(2 n + 2)} ( ξ _{x} )}{( 2 n + 2 )!} [\prod_{i = 0}^{n} (x - x_{i})]^{2}$ 。

6. 最小二乘法

最小二乘法是一种曲线拟合的方法。

对于一组数据 $y_{0} = f (x_{0})$ 、 $\dots$ 、 $y_{n} = f (x_{n})$ ，要确定 $y$ 与 $x$ 之间的近似表达式：（1）方法一：插值，在几何上插值曲线经过所有点；（2）方法二：曲线拟合，求一连续曲线 $y = φ (x)$ 使得误差 $θ = \sum_{i = 0}^{n} [φ (x_{i}) - y_{i}]^{2}$ （连续）或 $θ = 0 \leq i \leq n max ∣ φ (x_{i}) - y_{i} ∣$ （离散）达到最小（本来应该用误差的绝对值之和的，但是由于绝对值函数不便于分析，故改为平方和）。

插值法适用于数据精确或可靠度较高的情况，而曲线拟合法适用于数据本身就有误差的情况。当数据量特别大时一般不用插值法，因为高次插值效果不理想，而分段低次插值精度又不高，另外，由于测量数据本身就有误差，使插值曲线刻意经过这些点也不必要；而曲线拟合是根据草图画出一条拟合曲线（或使用低次多项式），然后根据最小二乘法求出该曲线，这条曲线未必经过所以已知点，但能反映处数据的基本趋势。

定义： $δ_{i} = φ (x_{i}) - y_{i}$ ，叫做偏差，在回归分析中又称为残差，其中 $φ_{i} (x)$ 是拟合函数， $y_{i}$ 是测量值。

最小二乘原则：使偏差的平方和最小。（顾名思义最小二乘法就是根据最小二乘原则选择拟合曲线 $y = φ (x)$ 的方法）

数学描述如下：设 $(x_{i}, y_{i})$ （ $i = 1, 2, \dots, m$ ）为给定的一组数据，要求在函数类 $Φ = {φ_{0} (x), φ_{1} (x), \dots, φ_{n} (x)}$ （ $n < m$ ）中，求一函数 $ϕ^{*} (x) = \sum_{j = 0}^{n} a_{j}^{*} φ_{j} (x)$ 满足 $∣∣ δ ∣ ∣_{2}^{2} = \sum_{i = 1}^{m} (φ^{*} (x_{i}) - y_{i})^{2} = min_{φ \in Φ} \sum_{i = 1}^{m} (φ (x_{i}) - y_{i})^{2}$ ，其中 $φ (x) = \sum_{j = 0}^{n} a_{j} φ_{j} (x)$ 为 $Φ$ 中任一函数。我们称 $φ^{*} (x)$ 为最小二乘解， $φ (x)$ 为拟合函数。

定义：对于给定的一组实验数据 $(x_{i}, y_{i})$ （ $x_{i}$ 互异， $i = 0, 1, \dots, m$ ），在函数类 $Φ = {φ_{0} (x), φ_{1} (x), \dots, φ_{n} (x)}$ （ $n < m$ 且 $φ_{1} (x), \dots, φ_{n} (x)$ 之间线性无关），存在唯一的函数 $φ^{*} (x) = a_{0}^{*} φ_{0} (x) + a_{1}^{*} φ_{1} (x) + \dots + a_{n}^{*} φ_{n} (x) = \sum_{j = 0}^{n} a_{j}^{*} φ_{j} (x)$ 使得关系式 $\sum_{i = 0}^{m} ω_{i} (φ^{*} (x_{i}) - y_{i})^{2} = min_{φ \in Φ} \sum_{i = 0}^{m} ω_{i} (φ (x_{i}) - y_{i})^{2}$ 成立（ $ω_{i}$ 是权系数，如果不是加权最小二乘法则等于1），并且其系数 $a_{k}^{*}$ （ $k = 0, 1, \dots, n$ ）可以通过解 $(φ_{0}, φ_{0}) (φ_{1}, φ_{0}) \dots (φ_{n}, φ_{0}) (φ_{0}, φ_{1}) (φ_{1}, φ_{1}) \dots (φ_{n}, φ_{1}) \dots \dots \dots \dots (φ_{0}, φ_{n}) (φ_{1}, φ_{n}) \dots (φ_{n}, φ_{n}) a_{0} a_{1} \dots a_{n} = (φ_{0}, f) (φ_{1}, f) \dots (φ_{n}, f)$ 得到，记作 $(h, g) = \sum_{i = 1}^{m} h (x_{i}) g (x_{i})$ 。

最小二乘法的多项式拟合形式的解法：对于给定的一组数据 $(x_{i}, y_{i})$ （ $i = 0, 1, \dots, m$ ），求作 $n$ 次多项式（ $n \leq m$ ） $φ^{*} (x) = a_{0}^{*} + a_{1}^{*} x + \dots + a_{n}^{*} x^{n} = \sum_{j = 0}^{n} a_{j}^{*} x^{k}$ 使得 $\sum_{i = 0}^{m} (φ^{*} (x_{i}) - y_{i})^{2} = min_{φ \in Φ} \sum_{i = 0}^{m} (φ (x_{i}) - y_{i})^{2}$ 成立，即拟合函数取为多项式，这样的曲线拟合问题又叫做多项式拟合问题，相应的法方程为 $m \sum_{i = 1}^{m} x_{i} \dots \sum_{i = 1}^{m} x_{i}^{n} \sum_{i = 1}^{m} x_{i} \sum_{i = 1}^{m} x_{i}^{2} \dots \sum_{i = 1}^{m} x_{i}^{n + 1} \dots \dots \dots \dots \sum_{i = 1}^{m} x_{i}^{n} \sum_{i = 1}^{m} x_{i}^{n + 1} \dots \sum_{i = 1}^{m} x_{i}^{2 n} a_{0} a_{1} \dots a_{n} = \sum_{i = 1}^{m} y_{i} \sum_{i = 1}^{m} x_{i} y_{i} \dots \sum_{i = 1}^{m} x_{i}^{n} y_{i}$ ，且该方程组存在唯一的一组解。

比如：直线拟合（一次多项式拟合）：若 $y = φ (x) = a_{0} + a_{1} x$ ，则 $a_{0}, a_{1}$ 可通过方程组 ${n a_{0} + (\sum_{k = 1}^{n} x_{k}) a_{1} = \sum_{k = 1}^{n} y_{k} (\sum_{k = 1}^{n} x_{k}) a_{0} + (\sum_{k = 1}^{n} x_{k}^{2}) a_{1} = \sum_{k = 1}^{n} x_{k} y_{k}$ 解出；

二次多项式拟合：若 $y = φ (x) = a_{0} + a_{1} x + a_{2} x^{2}$ ，则 $a_{0}, a_{1}, a_{2}$ 可通过方程组 $⎩ ⎨ ⎧ n a_{0} + (\sum_{k = 1}^{n} x_{k}) a_{1} + (\sum_{k = 1}^{n} x_{k}^{2}) a_{2} = \sum_{k = 1}^{n} y_{k} (\sum_{k = 1}^{n} x_{k}) a_{0} + (\sum_{k = 1}^{n} x_{k}^{2}) a_{1} + (\sum_{k = 1}^{n} x_{k}^{3}) a_{2} = \sum_{k = 1}^{n} x_{k} y_{k} (\sum_{k = 1}^{n} x_{k}^{2}) a_{0} + (\sum_{k = 1}^{n} x_{k}^{3}) a_{1} + (\sum_{k = 1}^{n} x_{k}^{4}) a_{2} = \sum_{k = 1}^{n} x_{k}^{2} y_{k}$ 解出。

例：已知一组实验数据：

$x_{k}$ $2$ $2.5$ $3$ $4$ $5$ $5.5$

$y_{k}$ $4$ $4.5$ $6$ $8$ $8.5$ $9$

试用直线拟合这组数据（计算过程保留3位小数）

解：设直线 $y = a_{0} + a_{1} x$ ，那么法方程公式为 ${n a_{0} + (\sum_{k = 1}^{n} x_{k}) a_{1} = \sum_{k = 1}^{n} y_{k} (\sum_{k = 1}^{n} x_{k}) a_{0} + (\sum_{k = 1}^{n} x_{k}^{2}) a_{1} = \sum_{k = 1}^{n} x_{k} y_{k}$ ，代入数据得 ${6 a_{0} + 22 a_{1} = 40 22 a_{0} + 90.5 a_{1} = 161.25$ ，解出 $a_{0} = 1.229$ 、 $a_{1} = 1.483$ ，故所求直线方程为 $y = 1.229 + 1.483 x$ 。

$x_{k}$	$2$	$2.5$	$3$	$4$	$5$	$5.5$
$y_{k}$	$4$	$4.5$	$6$	$8$	$8.5$	$9$

实际应用中，拟合函数可能是其他类型的函数，不一定用多项式进行拟合。对于非线性最小而二乘问题，一般有如下解法（1）化为线性最小二乘问题，部分可化为线性拟合问题的常见函数见下表；（2）马奎特(Marquardt)方法，计算机上求解非线性最小二乘拟合问题的最常用的方法之一（略）。

加权最小二乘法：有时数据中每一点的重要性可能是不一样的，因此设 $ω_{i}$ 表示数据点 $(x_{i}, y_{i})$ 的重度（重度：即权重或密度，统称为权系数），定义加权平方误差为 $∣∣ δ ∣ ∣_{2}^{2} = \sum_{i = 0}^{m} ω_{i} δ_{i}^{2} = \sum_{i = 0}^{m} ω_{i} [φ^{*} (x_{i}) - y_{i}]^{2}$ ，如果选用的拟合曲线为 $φ^{*} (x) = a_{0}^{*} + a_{1}^{*} x + \dots + a_{n}^{*} x^{n}$ ，那么相应的法方程为 $\sum_{i = 1}^{m} ω_{i} \sum_{i = 1}^{m} ω_{i} x_{i} \dots \sum_{i = 1}^{m} ω_{i} x_{i}^{n} \sum_{i = 1}^{m} ω_{i} x_{i} \sum_{i = 1}^{m} ω_{i} x_{i}^{2} \dots \sum_{i = 1}^{m} ω_{i} x_{i}^{n + 1} \dots \dots \dots \dots \sum_{i = 1}^{m} ω_{i} x_{i}^{n} \sum_{i = 1}^{m} ω_{i} x_{i}^{n + 1} \dots \sum_{i = 1}^{m} ω_{i} x_{i}^{2 n} a_{0} a_{1} \dots a_{n} = \sum_{i = 1}^{m} ω_{i} y_{i} \sum_{i = 1}^{m} ω_{i} x_{i} y_{i} \dots \sum_{i = 1}^{m} ω_{i} x_{i}^{n} y_{i}$ 。

例：已知一组数据 $(x_{i}, y_{i})$ 及权 $W_{i}$ 如下表，若 $x$ 与 $y$ 之间有线性关系 $y = a + b x$ ，试用最小二乘法确定系数 $a$ 和 $b$ .

$i$ 1 2 3 4

$W_{i}$ 14 27 12 1

$x_{i}$ 2 4 6 8

$y_{i}$ 2 11 28 40

解：将数据代入相应的法方程公式得 ${54 a + 216 b = 701 216 a + 984 b = 3580$ ，解出 $a = - 12.885$ 、 $b = 6.467$ 。

$i$	1	2	3	4
$W_{i}$	14	27	12	1
$x_{i}$	2	4	6	8
$y_{i}$	2	11	28	40

五、数值积分

1. 数值积分的概念

要近似计算积分 $I = \int_{a}^{b} f (x) d x$ ，可以用一个近似 $f (x)$ 且易于求原函数的函数代替 $f (x)$ 。

比如用拉格朗日插值多项式来近似 $f (x)$ ：在 $[a, b]$ 上取 $a \leq x_{0} < x_{1} < \dots < x_{n} \leq b$ ，做 $f$ 的 $n$ 次插值多项式 $L_{n} (x) = \sum_{k = 0}^{n} f (x_{k}) l_{k} (x)$ ，则 $\int_{a}^{b} f (x) d x \approx \int_{a}^{b} \sum_{k = 0}^{n} f (x_{k}) l_{k} (x) d x = \sum_{k = 0}^{n} f (x_{k}) \int_{a}^{b} l_{k} (x) d x = \sum_{k = 0}^{n} f (x_{k}) A_{k}$ ，其中 $A_{k} = \int_{a}^{b} \prod_{j = 0, j \neq = k}^{n} \frac{( x - x _{j} )}{( x _{k} - x _{j} )} d x$ 由节点决定，与 $f (x)$ 无关；

此时误差 $R_{n} [f] = \int_{a}^{b} f (x) d x - \sum_{k = 0}^{n} A_{k} f (x_{k}) = \int_{a}^{b} [f (x) - L_{n} (x)] d x = \int_{a}^{b} R_{n} (x) d x$ ，通过拉格朗日插值余项可得插值型求积公式的误差 $R_{n} [f] = \int_{a}^{b} \frac{f ^{(n + 1)} ( ξ )}{( n + 1 )!} \prod_{k = 0}^{n} (x - x_{k}) d x$ 。

例：对于 $[a, b]$ 上1次插值有 $L_{1} (x) = \frac{x - b}{a - b} f (a) + \frac{x - a}{b - a} f (b)$ ，则 $\int_{a}^{b} L_{1} (x) d x = \int_{a}^{b} \frac{x - b}{a - b} d x \cdot f (a) + \int_{a}^{b} \frac{x - a}{b - a} d x \cdot f (b)$ ， $A_{0} = \int_{a}^{b} \frac{x - b}{a - b} d x = \frac{b - a}{2}$ 、 $A_{1} = \int_{a}^{b} \frac{x - a}{b - a} d x = \frac{b - a}{2}$ ，故 $\int_{a}^{b} f (x) d x \approx \frac{b - a}{2} [f (a) + f (b)]$ 。

观察发现，上述积分的近似值其实就是由两个端点及坐标轴围成的梯形的面积，故我们将 $\frac{b - a}{2} [f (a) + f (b)]$ 称为梯形公式。

定义：若某个求积公式所对应的误差 $R [f]$ 满足：对任意 $k \leq n$ 阶的多项式有 $R [P_{k}] = 0$ 成立，且对某个 $n + 1$ 阶多项式有 $R [P_{n + 1}] \neq = 0$ 成立，则称此求积公式的代数精度为 $n$ 。

考察梯形公式 $\int_{a}^{b} f (x) d x \approx \frac{b - a}{2} [f (a) + f (b)]$ 的精度：

逐次检查公式是否精确成立，取 $f (x) = 1$ ，则左端为 $\int_{a}^{b} 1 d x = b - a$ ，右端 $\frac{b - a}{2} [1 + 1] = b - a$ ，左右两端相等，公式精确成立；

取 $f (x) = x$ ，则左端 $\int_{a}^{b} x d x = \frac{b ^{2} - a ^{2}}{2}$ ，右端 $\frac{b - a}{2} [a + b]$ ，左右两端仍相等，公式精确成立；

取 $f (x) = x^{2}$ ，则左端 $\int_{a}^{b} x^{2} d x = \frac{b ^{3} - a ^{3}}{3}$ ，右端 $\frac{b - a}{2} [a^{2} + b^{2}]$ ，此时左右两端不相等，故代数精度 $= 1$ 。

定理：求积公式 $\int_{a}^{b} f (x) d x \approx \sum_{k = 0}^{n} A_{k} f (x_{k})$ 至少有 $n$ 次代数精度的充要条件是该公式为插值型求积公式（即 $A_{k} = \int_{a}^{b} l_{k} (x) d x$ ，其中 $l_{i} (x)$ 是插值基函数）。

2. 牛顿-柯特斯公式及其复合

Newton-Cotes公式

Newton-Cotes(牛顿-柯特斯)公式：当节点等距分布时，即 $x_{k} = a + kh$ （其中 $h = \frac{b - a}{n}$ ， $k = 0, 1, \dots, n$ ），则 $A_{k} = \int_{x_{0}}^{x_{n}} \prod_{j = 0, j \neq = k}^{n} \frac{( x - x _{j} )}{( x _{k} - x _{j} )} d x = 令 x = a + t h \int_{0}^{n} \prod_{j = 0, j \neq = k}^{n} \frac{( t - j ) h}{( k - j ) h} \times h d t = \frac{( b - a ) ( - 1 ) ^{n - k}}{n k ! ( n - k )!} \int_{0}^{n} \prod_{j = 0, j \neq = k}^{n} (t - j) d t$ ，其中 $\frac{( - 1 ) ^{n - k}}{n k ! ( n - k )!} \int_{0}^{n} \prod_{j = 0, j \neq = k}^{n} (t - j) d t$ 称为Cotes系数，记作 $C_{k}^{(n)}$ ，因此有 $\int_{a}^{b} f (x) d x \approx \sum_{k = 0}^{n} f (x_{k}) A_{k} = \sum_{k = 0}^{n} f (x_{k}) (b - a) C_{k}^{(n)}$ ，该公式称为Newton-Cotes公式，简记为N-C公式。

显然Newton-Cotes公式是一种插值型求积公式。（其复合形式亦然，因为它是通过拉格朗日插值求积通过换元化简得到的）

Cotes系数仅取决于 $n$ 和 $k$ ，与 $f (x)$ 及区间 $[a, b]$ 均无关，其值可查表得到（如下表）。

当 $n$ 比较大时，Cotes公式会开始出现负项，计算过程的稳定性没有保证，因此一般只使用 $n \leq 4$ 的公式。

Newton-Cotes公式的截断误差可以通过拉格朗日插值多项式的截断误差推导出， $R_{n} [f] = \int_{a}^{b} R_{n} (x) d x = \int_{a}^{b} \frac{f ^{(n + 1)} ( ξ )}{( n + 1 )!} \prod_{k = 0}^{n} (x - x_{k}) d x = \frac{h ^{n + 2}}{( n + 1 )!} \int_{0}^{n} f^{(n + 1)} (ξ) [\prod_{j = 0}^{n} (t - j)] d t$ ， $ξ \in [a, b]$ ,，该误差公式表面 $h$ 越小，误差就越小。

Cotes系数的性质：

（1）权性： $\sum_{k = 0}^{n} C_{k}^{(n)} = 1$ ，即每一行的Cotes系数之和等于1；

（2）对称性： $C_{k}^{(n)} = C_{n - k}^{(n)}$ ；

一些比较常用的Newton-Cotes公式：

$n = 1$ 时， $C_{0}^{(1)} = \frac{1}{2}$ 、 $C_{1}^{(1)} = \frac{1}{2}$ ，故 $A_{0} = A_{1} = \frac{b - a}{2}$ ，因此 $\int_{a}^{b} f (x) d x \approx \frac{b - a}{2} [f (a) + f (b)]$ ，即前面的梯形公式，代数精度为1，误差 $R_{n} [f] = \int_{a}^{b} \frac{f ^{''} ( ξ _{x} )}{2 !} (x - a) (x - b) d x = 令 x = a + t h, 根据中值定理 - \frac{1}{12} h^{3} f^{''} (ξ)$ ，其中 $ξ \in [a, b]$ 、 $h = \frac{b - a}{1}$ ；
$n = 2$ 时，查表得 $\int_{a}^{b} f (x) d x \approx \frac{b - a}{6} [f (a) + 4 f (\frac{a + b}{2}) + f (b)]$ ，该公式称为Simpson公式，代数精度为3，误差 $R [f] = - \frac{1}{90} h^{5} f^{(4)} (ξ)$ ，其中 $ξ \in [a, b]$ 、 $h = \frac{b - a}{2}$ ；
$n = 3$ 时，查表得 $\int_{a}^{b} f (x) d x \approx \frac{b - a}{8} [f (x_{0}) + 3 f (x_{1}) + 3 f (x_{2}) + f (x_{3})]$ ，该公式称为Simpson $\frac{3}{8}$ 公式（或第二Simpson公式），代数精度为3，误差 $R [f] = - \frac{3}{80} h^{5} f^{(5)} (ξ)$ ；
$n = 4$ 时，得到Cotes求积公式，代数精度为5，误差 $R [f] = - \frac{8}{945} h^{7} f^{(6)} (ξ)$ 。

定理： $n$ 为偶数阶的Newton-Cotes公式至少有 $n + 1$ 次代数精度。

例：试分别用梯形公式和Simpson公式计算积分 $I = \int_{0}^{1} e^{x} d x$ 的近似值，并估计截断误差.

解：梯形公式： $I \approx \frac{b - a}{2} [f (a) - f (b)] = \frac{e ^{0} + e ^{1}}{2} \approx 1.85914$ ，截断误差 $∣ R_{1} [f] ∣ = \frac{1}{12} ∣ h^{3} f^{''} (ξ) ∣ \leq \frac{( 1 - 0 ) ^{3}}{12} max_{0 \leq x \leq 1} ∣ f^{''} (x) ∣ = \frac{e}{12} \approx 0.22652$ （ $f^{''} (ξ)$ 一般用二阶导绝对值的最大值替代）；

Simpson公式： $I \approx \frac{b - a}{6} [f (a) + 4 f (\frac{a + b}{2}) + f (b)] = \frac{1}{6} (e^{0} + 4 e^{\frac{1}{2}} + e^{1}) \approx 1.71886$ ，截断误差 $R [f] = - \frac{( b - a ) ^{5}}{2880} max_{0 \leq x \leq 1} ∣ f^{(4)} (x) ∣ \leq \frac{e}{2880} \approx 0.00095$ 。

显然Simpson公式的结果优于梯形公式的结果。

Newton-Cotes复合求积公式

根据Newton-Cotes公式的截断误差公式，步长 $h$ 越小，截断误差就越小，为了得到比较精确的积分值，必须用高阶求积公式，但高阶（ $n \geq 8$ ）的Newton-Cotes公式是不稳定的，因此不能通过提高Newton-Cotes公式阶数的方法来提高精度。

我们考虑对被积函数采用分段低次多项式插值（分段通常是等分），这样就得到了分段低次合成的Newton-Cotes复合求积公式，比如：

复合梯形公式：令 $h = \frac{b - a}{n}$ 、 $x_{k} = a + kh$ （ $k = 0, \dots, n$ ），在每个 $[x_{k - 1}, x_{k}]$ 上用梯形公式，得到 $\int_{x_{k - 1}}^{x_{k}} f (x) d x \approx \frac{x _{k} - x _{k - 1}}{2} [f (x_{k - 1}) + f (x_{k})]$ ，其中 $k = 1, \dots, n$ ；再把每个区间的值加起来得到复合梯形公式 $\int_{a}^{b} f (x) d x \approx \sum_{k = 1}^{n} \frac{h}{2} [f (x_{k - 1}) + f (x_{k})] = \frac{h}{2} [f (a) + 2 \sum_{k = 1}^{n - 1} f (x_{k}) + f (b)]$ ，记作 $T_{n}$ ，误差 $R [f] = - \frac{h ^{2}}{12} (b - a) f^{''} (ξ)$ ，其中 $ξ \in (a, b)$ ；
复合Simpson公式： $x_{k}$ 与 $x_{k + 1}$ 的中点用 $x_{k + \frac{1}{2}}$ 表示，则 $\int_{x_{k}}^{x_{k + 1}} f (x) d x \approx \frac{h}{6} [f (x_{k}) + 4 f (x_{k + \frac{1}{2}}) + f (x_{k + 1})]$ ，把每个区间的值加起来得复合Simpson公式 $\int_{a}^{b} f (x) d x \approx \frac{h}{6} [f (a) + 4 \sum_{k = 0}^{n - 1} f (x_{k + \frac{1}{2}}) + 2 \sum_{k = 0}^{n - 1} f (x_{k + 1}) + f (b)]$ ，记作 $S_{n}$ ，误差 $R [f] = - \frac{b - a}{180} (\frac{h}{2})^{4} f^{(4)} (ξ)$ ，其中 $ξ \in (a, b)$ ；

由于公式中出现了 $x_{\frac{1}{2}}$ 、 $x_{\frac{3}{2}}$ 这样的半节点，为方便编程实现，令 $n^{'} = 2 n$ 为偶数，这时 $h^{'} = \frac{b - a}{n ^{'}} = \frac{h}{2}$ ， $x_{k} = a + k h^{'}$ ，于是 $S_{n} = \frac{h ^{'}}{3} [f (a) + 4 \sum_{odd k} f (x_{k}) + 2 \sum_{even k} f (x_{k}) + f (b)]$ 。
复合Cotes公式：类似地， $\int_{x_{k}}^{x_{k + 1}} f (x) d x \approx \frac{h}{90} [7 f (x_{k}) + 32 f (x_{k + \frac{1}{4}}) + 12 f (x_{k + \frac{1}{2}}) + 32 f (x_{k + \frac{3}{4}}) + 7 f (x_{k + 1})]$ ，加起来得复合Cotes公式 $\int_{a}^{b} f (x) d x \approx \frac{h}{90} [7 f (a) + 14 \sum_{k = 1}^{n - 1} f (x_{k}) + 32 \sum_{k = 0}^{n - 1} f (x_{k + \frac{1}{4}}) + 12 \sum_{k = 0}^{n - 1} f (x_{k + \frac{1}{2}}) + 32 \sum_{k = 0}^{n - 1} f (x_{k + \frac{3}{4}}) + 7 f (b)]$ ，记作 $C_{n}$ ，误差 $R_{c} [f] = \int_{a}^{b} f (x) d x - C_{n} = - \frac{2 ( b - a )}{945} (\frac{h}{4})^{6} f^{(6)} (ξ)$ ，其中 $ξ \in (a, b)$ 。

定义：若一个积分公式的误差满足 $h \to 0 lim \frac{R [ f ]}{h ^{p}} = C < \infty$ ，且 $C \neq = 0$ 是常数，则称该公式是p阶收敛的（收敛阶越高，收敛速度越快）。

前述三种复合积分公式的收敛阶数： $T_{n} \sim O (h^{2})$ 、 $S_{n} \sim O (h^{4})$ 、 $C_{n} \sim O (h^{6})$ 。

例：计算 $π = \int_{0}^{1} \frac{4}{1 + x ^{2}} d x$ ，函数值如下表

$x$ $0$ $\frac{1}{8}$ $\frac{2}{8}$ $\frac{3}{8}$ $\frac{4}{8}$ $\frac{5}{8}$ $\frac{6}{8}$ $\frac{7}{8}$ 1

$f (x)$ 4 3.938461538 3.764705882 3.506849315 3.2 2.876404494 2.56 2.265486726 2

解：令 $x_{k} = \frac{k}{8}$ ，则 $T_{8} = \frac{1}{16} [f (0) + 2 \sum_{k = 1}^{7} f (x_{k}) + f (1)] = 3.1 38988494$ ，只有2位有效数字，而如果使用复合Simpson公式 $S_{4} = \frac{1}{24} [f (0) + 4 \sum_{odd} f (x_{k}) + 2 \sum_{even} f (x_{k}) + f (1)] = 3.141592 502$ ，差不多的计算量却有7位有效数字，实际上 $T_{215} = 3.141592 02$ 才有7位有效数字。

通过上例的分析，有必要讨论给定精度 $ε$ ，要如何取 $n$ 才能达到精度要求。

$x$	$0$	$\frac{1}{8}$	$\frac{2}{8}$	$\frac{3}{8}$	$\frac{4}{8}$	$\frac{5}{8}$	$\frac{6}{8}$	$\frac{7}{8}$	1
$f (x)$	4	3.938461538	3.764705882	3.506849315	3.2	2.876404494	2.56	2.265486726	2

复合梯形公式的误差分析：对于复合梯形公式，要求 $∣ I - T_{n} ∣ < ε$ ，则 $R [f] = - \frac{h ^{2}}{12} (b - a) f^{''} (ξ) = - \frac{h ^{2}}{12} \sum_{k = 1}^{n} [f^{n} (ξ_{k}) \cdot h] \approx - \frac{h ^{2}}{12} \int_{a}^{b} f^{''} (x) d x = - \frac{h ^{2}}{12} [f^{'} (b) - f^{'} (a)]$ 。（其他公式类似）

上例若要求 $∣ I - T_{n} ∣ < 1 0^{- 6}$ ，则 $∣ R_{n} [f] ∣ \approx \frac{h ^{2}}{12} ∣ f^{'} (1) - f^{'} (0) ∣ = \frac{h ^{2}}{6} < 1 0^{- 6}$ ，解出 $h < 0.00244949$ ，即 $n = \frac{b - a}{h}$ 应取409。

通常采用将区间不断对分的方法，即取 $n = 2^{k}$ 。上例中 $2^{k} \geq 409$ ，即 $k = 9$ 时， $T_{512} = 3.141592 02$ 满足精度要求。

根据误差公式 $R [f]$ ，不难发现每次对分一次，误差大概会变成原来的 $\frac{1}{4}$ （ $n$ 变为 $2 n$ ，则 $h$ 变为 $\frac{b - a}{2 n} = \frac{h}{2}$ ，代入 $R [f]$ 中），即 $\frac{I - T _{2 n}}{I - T _{n}} \approx \frac{1}{4} \Rightarrow I - T_{2 n} \approx \frac{1}{3} (T_{2 n} - T_{n})$ ，如果采用区间一直对分的方法，可以使用此公式判断是否停止。

3. 变步长求积

前面介绍的Newton-Cotes公式是固定步长的方法，但是往往题目中只会给出精度要求，要根据精度要求求出步长的话要通过 $R [f]$ ，这时就会出现 $f (x)$ 的高阶导数，不便于计算。

变步长求积思想：变步长方法就是根据规定的精度要求，在计算过程中将积分区间逐次分半，每对分一次就用同一种复合求积公式计算出相应的积分近似值，并同时查看相继两次计算结果的误差是否达到要求，直到所求得的积分近似值满足精度要求为止。（也就是说，点仍然是等距分布的，只是每次迭代都再对分一次区间）

以梯形公式为例，将其递推化为变步长的梯形公式：

（1）将积分区间 $[a, b]$ 等分，等分点 $x_{k} = a + kh$ （其中 $k = 0, 1, \dots, n$ 、 $h = \frac{b - a}{n}$ ），复合梯形公式为 $T_{n} = \frac{h}{2} [f (a) + 2 \sum_{k = 1}^{n - 1} f (x_{k}) + f (b)]$ ，画图可知是两个端点的函数值加一次，中间点的函数值全部加两次，最后再乘以 $\frac{h}{2}$ ；

（2）若区间精度达不到要求，则将 $[x_{k}, x_{k + 1}]$ 每个小区间二分一次，即令区间中点为 $x_{k + \frac{1}{2}} = \frac{x _{k} + x _{k + 1}}{2}$ ，该小区间二分前后的两个积分值分别记为 $T_{k 1}$ 和 $T_{k 2}$ ，代入梯形公式结合图像可知 $T_{k 1} = \frac{h}{2} [f (x_{k}) + f (x_{k + 1})]$ 、 $T_{k 2} = \frac{\frac{h}{2}}{2} [f (x_{k}) + 2 f (x_{k + \frac{1}{2}}) + f (x_{k + 1})]$ ，整理后得 $T_{k 2} = \frac{1}{2} T_{k 1} + \frac{h}{2} f (x_{k + \frac{1}{2}})$ ；

（3）将二分后的所有小区间积分累加 $\sum_{k = 0}^{n - 1} T_{k 2} = \sum_{k = 0}^{n - 1} [\frac{1}{2} T_{k 1} + \frac{h}{2} f (x_{k + \frac{1}{2}})]$ ，于是得到公式 $T_{2 n} = \frac{1}{2} T_{n} + \frac{h}{2} \sum_{k = 0}^{n - 1} f (x_{k + \frac{1}{2}})$ ，这就是梯形公式的递推化公式，也称为梯形公式的变步长法，式中 $h$ 表示二分前的步长，即 $h = \frac{b - a}{n}$ ；

（4）计算过程中，常用 $∣ T_{2 n} - T_{n} ∣ < ε$ 来判断精度是否满足要求。

可见，计算 $T_{2 n}$ （分成 $2 n$ 份时的积分）时用到了 $T_{n}$ 的结果，只需对新的中间节点求和即可，这样可以使计算量减小很多。

例：用变步长法计算 $\int_{0}^{1} \frac{s i n x}{x} d x$ ，使误差不超过 $1 0^{- 6}$ （假设 $\frac{s i n x}{x}$ 的函数值已知）.

解：先对整个积分区间 $[0, 1]$ 使用梯形公式，得 $T_{1} = \frac{1}{2} [f (0) + f (1)] = 0.920735492$ ，通过 $T_{2 n} = \frac{1}{2} T_{n} + \frac{h}{2} \sum_{k = 0}^{n - 1} f (x_{k + \frac{1}{2}})$ 有：

$k$ $n$ $T_{n}$ $∣ T_{n} - T_{\frac{n}{2}} ∣$

0 1 0.920735492

1 2 0.939793285 0.019057793

2 4 0.944513522 0.004720237

3 8 0.945690864 0.001177342

4 16 0.94598503 0.000294166

5 32 0.946058561 7.3531 $\times 1 0^{- 5}$

6 64 0.946076943 1.8382 $\times 1 0^{- 5}$

7 128 0.946081539 4.596 $\times 1 0^{- 6}$

8 256 0.946082687 1.148 $\times 1 0^{- 6}$

9 512 0.946082975 2.88 $\times 1 0^{- 7}$

故 $T_{512} = 0.946082975$ 。

$k$	$n$	$T_{n}$	$∣ T_{n} - T_{\frac{n}{2}} ∣$
0	1	0.920735492
1	2	0.939793285	0.019057793
2	4	0.944513522	0.004720237
3	8	0.945690864	0.001177342
4	16	0.94598503	0.000294166
5	32	0.946058561	7.3531 $\times 1 0^{- 5}$
6	64	0.946076943	1.8382 $\times 1 0^{- 5}$
7	128	0.946081539	4.596 $\times 1 0^{- 6}$
8	256	0.946082687	1.148 $\times 1 0^{- 6}$
9	512	0.946082975	2.88 $\times 1 0^{- 7}$

4.龙贝格(Romberg)积分

变步长梯形求积虽然算法简单，但是收敛速度仍较慢。（因为它本质上仍是梯形法，只是不需要每次分区间都再从头算一遍）

前面提到了梯形公式的误差有如下结论： $\frac{I - T _{2 n}}{I - T _{n}} \approx \frac{1}{4} \Rightarrow I - T_{2 n} \approx \frac{1}{3} (T_{2 n} - T_{n})$ 。如果用 $\frac{1}{3} (T_{2 n} - T_{n})$ 对 $T_{2 n}$ 进行误差补偿，即 $\overline{T} = T_{2 n} + \frac{1}{3} (T_{2 n} - T_{n}) = \frac{4}{3} T_{2 n} - \frac{1}{3} T_{n}$ ，那么由 $I \approx \frac{4 T _{2 n} - T _{n}}{4 - 1} = \frac{4}{3} T_{2 n} - \frac{1}{3} T_{n}$ 来计算 $I$ 效果是否好一些？事实上，容易验证 $\frac{4 T _{2 n} - T _{n}}{4 - 1} = S_{n}$ ，即用梯形法二分前后的两个积分值 $T_{n}$ 和 $T_{2 n}$ 按上式作线性组合，却得到了复合Simpson公式 $S_{n}$ ，精度提高了。

再考察Simpson公式，其截断误差与 $h^{4}$ 成正比，因此如果将步长折半，则误差减至 $\frac{1}{16}$ ，即 $\frac{I - S _{2 n}}{I - S _{n}} \approx \frac{1}{16}$ ，由此可得 $I \approx \frac{16 S _{2 n} - S _{n}}{15}$ ，容易验证， $\frac{4 ^{2} S _{2 n} - S _{n}}{4 ^{2} - 1} = C_{n}$ ，即右端项就是复合Cotes公式。依此类推， $\frac{4 ^{3} C _{2 n} - C _{n}}{4 ^{3} - 1} = R_{n}$ ，其中 $R_{n}$ 是Romberg序列。

根据前面分析，得到公式一： $\frac{4 T _{2 n} - T _{n}}{4 - 1} = S_{n}$ ，公式二： $\frac{4 ^{2} S _{2 n} - S _{n}}{4 ^{2} - 1} = C_{n}$ ，公式三： $\frac{4 ^{3} C _{2 n} - C _{n}}{4 ^{3} - 1} = R_{n}$ ，于是有Romberg算法：

（1）用梯形公式计算积分近似值 $T_{1}$ ，令 $T_{1} = T_{0}^{(0)}$ ；

（2）按变步长梯形公式计算 $T_{2}$ ，令 $T_{2} = T_{0}^{(1)}$ ；

（3）将 $T_{1}$ 、 $T_{2}$ 按公式一得到 $S_{1} = T_{1}^{(0)}$ ，判断相邻两次积分值 $∣ S_{1} - T_{2} ∣$ 是否小于 $ε$ ，若满足精度要求则停止计算并输出 $S_{1}$ ，否则进行下一步；

（4）按变步长梯形公式计算 $T_{4}$ ，令 $T_{4} = T_{0}^{(2)}$ ；

（5）将 $T_{2}$ 和 $T_{4}$ 按公式一得到 $S_{2} = T_{1}^{(1)}$ ；

（6）将 $S_{1}$ 、 $S_{2}$ 按格式二得到 $C_{1} = T_{2}^{(0)}$ ，判断 $∣ C_{1} - S_{2} ∣$ 是否小于 $ε$ ，若满足则停止计算并输出 $C_{1}$ ，否则进行下一步；

（7）按变步长梯形公式计算 $T_{8}$ ，令 $T_{8} = T_{0}^{(3)}$ ；

（8）将 $T_{4}$ 、 $T_{8}$ 按公式一得到 $S_{4} = T_{1}^{(2)}$ ；

（9）将 $S_{2}$ 、 $S_{4}$ 按公式二得到 $C_{2} = T_{2}^{(1)}$ ；

（10）将 $C_{1}$ 、 $C_{2}$ 按格式三得到 $R_{1} = T_{3}^{(0)}$ ；判断 $∣ R_{1} - C_{2} ∣$ 是否小于 $ε$ ，若满足则停止计算并输出 $R_{1}$ ，否则类似进行下去；

可画出过程如下：

例：用Romber求积算法计算 $I = \int_{0}^{1} \frac{s i n x}{x} d x$ 的近似值，使其具有6位有效数字.

解：仿照过程的画图表示，列表如下

$k$ $T_{2}^{k}$ $S_{2}^{k - 1}$ $C_{2}^{k - 2}$ $R_{2}^{k - 3}$

0 0.9207355

1 0.9397933 0.9461459

2 0.9445135 0.9460869 0.9400830

3 0.94456909 0.9460833 0.9460831 0.9460831

故 $I \approx 0.9460831$ 。（精确值为 $0.946083070367 \dots$ ，即计算结果前6位确实是有效数字）

$k$	$T_{2}^{k}$	$S_{2}^{k - 1}$	$C_{2}^{k - 2}$	$R_{2}^{k - 3}$
0	0.9207355
1	0.9397933	0.9461459
2	0.9445135	0.9460869	0.9400830
3	0.94456909	0.9460833	0.9460831	0.9460831

六、常微分方程的数值解法

1. 初值问题介绍

微分方程：有一个或多个导数及其函数的方程称为微分方程，它分为①常微分方程-只有一个自变量；②偏微分方程-一个以上自变量。

实际中求解常微分方程的定解问题有两类：初值问题和边值问题。（定解指已知因变量或其导数在某些点上是已知的）

边解问题（边值问题）：约束条件已知，在自变量的任一非初值上，已知函数值或其导数值，如 ${y^{''} = f (x, y, y^{'}) y (a) = α, y (b) = β$ ，常常可以将边解问题转化为初值问题求解。

初值问题：约束条件为在自变量的初值上已知函数值，如 ${y^{'} = f (x, y) y (x_{0}) = y_{0} \Rightarrow {\frac{d y}{d x} = f (x, y) \frac{x}{x _{0}} y (x_{0}) = y_{0}$ ，要求解 $y (x)$ ，可在 $a \leq x_{0}^{'} \leq x_{1}^{'} \leq \dots \leq x_{n}^{'} \leq b$ 上的 $y (x_{i})$ 的近似值 $y_{i} = y (x_{i})$ （ $i = 0, 1, \dots, n$ ），通常取等距节点，令 $h = x_{i + 1} - x_{i}$ ，则 $x_{i} = x_{0} + ih$ （ $i = 0, 1, \dots, n$ ）。初值问题的数值解法特点：按节点顺序依次推进，若已知 $y_{0}, y_{1}, \dots, y_{i}$ ，可根据递推公式求出 $y_{i + 1}$ 。

初值问题的常用解法：（1）单步法：利用前一个单步的信息（即一个点），在 $y = f (x)$ 上找下一点 $y_{i}$ ，欧拉法、龙格-库塔法都属于单步法；（2）预测校正法：又称多步法，利用一个以上的节点信息求 $f (x)$ 的下一个 $y_{i}$ ，常用迭代法，如改进欧拉法、阿姆斯特当法。

单步递推法的基本思想是从 $(x_{i}, y_{i})$ 出发，以某一斜率沿直线达到 $(x_{i + 1}, y_{i + 1})$ 。

比如，方程 $x y^{'} - 2 y = 4 x \Rightarrow y^{'} = \frac{2 y}{x} + 4$ ，令 $f (x, y) = \frac{2 y}{x} + 4$ ，且给出初值 $y (1) = - 3$ ，就可以得到一阶常微分方程的初值问题 ${\frac{d y}{d x} = f (x, y) = \frac{2 y}{x} + 4 y (1) = - 3$ 。

只要函数 $f (x, y)$ 适当光滑连续，且关于 $y$ 满足李普希兹(Lipschitz)条件（即存在常数 $L$ ，使得 $∣ f (x, y) - f (x, \overline{y}) ∣ \leq L ∣ y - \overline{y} ∣$ ），则由常微分方程理论知，初值问题的解必存在且唯一。

数值解法的含义：设法将常微分方程离散化，建立差分方程，给出解在一些离散点上的近似值。

数学描述如下：设方程问题的解 $y (x)$ 的存在区间是 $[a, b]$ ，令 $a = x_{0} < x_{1} < \dots < x_{n} = b$ ，其中 $h_{k} = x_{k + 1} - x_{k}$ （如果是等距节点 $h = \frac{b - a}{n}$ 则称为步长），由于 $y (x)$ 的解析表达式不容易得到或无法得到，我们用数值方法求得 $y (x)$ 在每个节点 $x_{k}$ 上 $y (x_{k})$ 的近似值，用 $y_{k}$ 表示，即 $y_{k} \approx y (x_{k})$ ，这样 $y_{0}, y_{1}, \dots, y_{n}$ 称为微分方程的数值解。

2. 欧拉(Euler)方法、欧拉预校法

欧拉显格式计算公式（显式欧拉公式或显式欧拉格式）： $y_{i + 1} = y_{i} + h f (x_{i}, y_{i})$ ，它可循环求解。

推导：Taylor展开法， $y (x_{i + 1}) = y (x_{i}) + y^{'} (x_{i}) (x_{i + 1} - x_{i}) + \frac{y ^{''} ( ξ )}{2 !} (x_{i + 1} - x_{i})^{2} = y (x_{i}) + h f (x_{i}, y (x_{i})) + \frac{h ^{2}}{2} y^{''} (ξ_{i})$ ，其中二阶导的部分可忽略，于是得到 $y_{i + 1} = y_{i} + h f (x_{i}, y_{i})$ （其中 $i = 0, 1, 2, \dots, n - 1$ ，且根据初值问题的定义有 $\frac{d y}{d x} = f (x, y)$ ）。

或者用向前差商近似导数 $y^{'} (x_{0}) \approx \frac{y ( x _{1} ) - y ( x _{0} )}{h}$ ，即 $y (x_{1}) \approx y (x_{0}) + h y^{'} (x_{0}) = y_{0} + h f (x_{0}, y_{0})$ ，将其记作 $y_{1}$ ，依此类推得到 $y_{i + 1} = y_{i} + h f (x_{i}, y_{i})$ ，这是显式格式的欧拉计算公式。

隐式欧拉法（隐式欧拉公式或隐式欧拉格式）： $y_{i + 1} = y_{i} + h f (x_{i + 1}, y_{i + 1})$ 。

推导：用向后差商近似导数 $y^{'} (x_{1}) \approx \frac{y ( x _{1} ) - y ( x _{0} )}{h}$ ，即 $y (x_{1}) \approx y_{0} + h f (x_{1}, y (x_{1}))$ ，依此类推得到 $y_{i + 1} = y_{i} + h f (x_{i + 1}, y_{i + 1})$ ，其中 $i = 0, 1, 2, \dots, n - 1$ 。由于未知数 $y_{i + 1}$ 同时出现在等式两边，不能直接得到，故称为隐式欧拉公式。使用时常用显式计算一个初值，在迭代求解。

例：用欧拉法解初值问题 ${\frac{d y}{d x} = 1 - x y y (0) = 0, 0 < x < 1$ ，取步长为 $h = 0.2$ ，计算过程保留4位小数.

解：令 $f (x, y) = 1 - x y$ ，根据欧拉公式 $y (x_{i + 1}) \approx y_{i + 1} = y_{i} + h f (x_{i}, y_{i}) = y_{i} + 0.2 (1 - x_{i} y_{i})$ ，其中 $i = 0, 1, 2, 3, 4, 5$ ，得

$i$ 0 1 2 3 4 5

$x_{i}$ 0 0.2 0.4 0.6 0.8 1

$y_{i}$ 0 0.2000 0.3920 0.5606 0.6934 0.7824

$i$	0	1	2	3	4	5
$x_{i}$	0	0.2	0.4	0.6	0.8	1
$y_{i}$	0	0.2000	0.3920	0.5606	0.6934	0.7824

除了欧拉格式外，还可以用梯形格式： $y_{i + 1} = y_{i} + \frac{h}{2} [f (x_{i}, y_{i}) + f (x_{i + 1}, y_{i + 1})]$ （其中 $i = 0, 1, 2, \dots, n - 1$ ），它是隐格式的。

观察可知，它是显、隐式两种算法的平均，它有两种计算思路（1）化为显格式；（2）用迭代法求 $y_{i + 1}$ ，初值用Euler显格式确定。

中点欧拉公式： $y_{i + 1} = y_{i - 1} + 2 h f (x_{i}, y_{i})$ 。

使用中心差商近似导数 $y^{'} (x_{1}) \approx \frac{y ( x _{2} ) - y ( x _{0} )}{2 h}$ ，即 $y (x_{2}) \approx y (x_{0}) + 2 h f (x_{1}, y (x_{1}))$ ，依此类推得 $y_{i + 1} = y_{i - 1} + 2 h f (x_{i}, y_{i})$ ，其中 $i = 1, \dots, n - 1$ 。

其实不同的迭代公式就是在每个小区间上用不同的形状的面积来近似积分值。

它们的优缺点如下：显式欧拉简单但精度低，隐式欧拉稳定性最好但精度低且计算量大，梯形公式精度高但计算量大，中点公式精度较高但需要多一个初值，可能影响精度。

为加快收敛速度，将显式欧拉格式与梯形格式结合起来就得到了欧拉预估-校正法。

欧拉预估-校正法（改进欧拉法）步骤：

（1）用显式欧拉公式作预测，算出 $\overline{y_{n + 1}} = y_{i} + h f (x_{i}, y_{i})$ ；

（2）将 $\overline{y_{i + 1}}$ 代入隐式梯形公式的右半部分作校正，得到 $y_{i + 1} = y_{i} + \frac{h}{2} [f (x_{i}, y_{i}) + f (x_{i + 1}, \overline{y_{i + 1}})]$ 。

综合以上两步，可得到 $y_{i + 1} = y_{i} + \frac{h}{2} [f (x_{i}, y_{i}) + f (x_{i + 1}, y_{i} + h f (x_{i}, y_{i}))]$ ，其中 $i = 0, \dots, n - 1$ 。易证，当 $\frac{n}{2} ∣ \frac{\partial f}{\partial y} ∣ \leq 2 < 1$ 时，该方法关于 $i$ 的迭代收敛。

欧拉预估-校正法也可以写成 $⎩ ⎨ ⎧ y_{i + 1} = y_{i} + \frac{h}{2} (K_{1} + K_{2}) K_{1} = f (x_{i}, y_{i}) K_{2} = f (x_{i} + h, y_{i} + h K_{1}) y_{0} = α$ ，其中 $i = 0, \dots, n - 1$ 。

例：求初值问题 ${y^{'} = y - \frac{2 x}{y} y (0) = 1, 0 \leq x \leq 1$ 的数值解，取步长 $h = 0.1$ 。（精确解为 $y (x) = (1 + 2 x)^{\frac{1}{2}}$ ）

解： $⎩ ⎨ ⎧ y_{i + 1} = y_{i} + 0.05 (K_{1} + K_{2}) K_{1} = y_{i} - \frac{2 x _{i}}{y _{i}} K_{2} = y_{i} + 0.1 K_{1} - \frac{2 ( x _{i} + 0.1 )}{y _{i} + 0.1 K _{1}} y_{0} = 1$ ，其中 $i = 0, 1, \dots, 9$ ，计算结果如下表（为对比不同方法，加入了欧拉格式和精确解）

$i$ $x_{i}$ 欧拉法 $y_{i}$ 欧拉预校法 $y_{i}$ 精确解 $y (x_{i})$

0 0 1 1 1

1 0.1 1.1 1.095909 1.095445

2 0.2 1.191818 1.184096 1.183216

3 0.3 1.277438 1.266201 1.264991

4 0.4 1.358213 1.343360 1.341641

5 0.5 1.435133 1.416402 1.414214

6 0.6 1.508966 1.485956 1.483240

7 0.7 1.580338 1.552515 1.549193

8 0.8 1.649783 1.616476 1.612452

9 0.9 1.717779 1.678168 1.673320

10 1 1.784770 1.737869 1.732051

$i$	$x_{i}$	欧拉法 $y_{i}$	欧拉预校法 $y_{i}$	精确解 $y (x_{i})$
0	0	1	1	1
1	0.1	1.1	1.095909	1.095445
2	0.2	1.191818	1.184096	1.183216
3	0.3	1.277438	1.266201	1.264991
4	0.4	1.358213	1.343360	1.341641
5	0.5	1.435133	1.416402	1.414214
6	0.6	1.508966	1.485956	1.483240
7	0.7	1.580338	1.552515	1.549193
8	0.8	1.649783	1.616476	1.612452
9	0.9	1.717779	1.678168	1.673320
10	1	1.784770	1.737869	1.732051

3. 误差估计、收敛性、稳定性

定义：假设 $y_{i} = y (x_{i})$ ，按某种方法由 $y_{i}$ 严格算出 $y_{i + 1}$ 这一步的误差叫做局部截断误差，记为 $R_{i + 1} = y (x_{i + 1}) - y_{i + 1}$ 。

定义：由 $y_{0}$ 逐步严格算出 $y_{i + 1}$ 的误差称为整体截断误差，记作 $ε_{i + 1} = y (x_{i + 1}) - y_{i + 1}$ 。

定义：若某种方法的局部截断误差为 $R_{i + 1} = O (h^{p + 1})$ ，则称该方法为p阶方法。（显然 $p$ 越大越好）

比如，欧拉显格式的局部误差：完整的Taylor展开为 $y (x_{i + 1}) = y (x_{i}) + h f (x_{i}, y (x_{i})) + \frac{h ^{2}}{2} y^{''} (ξ_{i})$ ，假设 $y_{i} = y (x_{i})$ ，则根据欧拉显格式 $y_{i + 1} = y_{i} + h f (x_{i}, y_{i}) = y (x_{i}) + h f (x_{i}, y (x_{i}))$ ，两式相减得 $R_{i + 1} = y (x_{i + 1}) - y_{i + 1} = \frac{h ^{2}}{2} y^{''} (ε_{n}) = O (h^{2})$ ，所以欧拉显格式是一阶方法。

同理，梯形法 $R_{i + 1} = O (h^{3})$ ，欧拉预校法 $R_{i + 1} = O (h^{3})$ 。

定理：设 $f (x, y)$ 满足关于第2个变元的Lipschitz条件 $∣ f (x_{1}, y_{1}) - f (x_{2}, y_{2}) ∣ \leq L ∣ y_{1}, y_{2} ∣$ ，且欧拉显格式的局部截断误差满足 $∣ R_{i} ∣ \leq \frac{i ^{2}}{2} M_{2}$ ，则 $∣ ε_{i} ∣ \leq e^{(b - a)} ∣ ε_{0} ∣ + \frac{M _{2}}{2 L} [e^{(b - a)} - 1] h$ ，其中 $i = 1, 2, \dots$ ， $M_{2} = a \leq x \leq b max ∣ y^{''} (x) ∣$ ， $[a, b]$ 为求解区间。

若初始误差为 $ε_{0} = 0$ ，则 $∣ ε_{i} ∣$ 比 $∣ R_{i} ∣$ 低一阶（具有普遍意义）。该定理揭示了局部与整体截断误差的关系，而讨论局部截断误差相对简单，所以我们这里只讨论局部截断误差。

定义：对某方法，任意固定 $x = x_{i} = x_{0} + ih$ ，当 $h \to 0$ （同时 $i \to \infty$ ）时，若有 $ε_{i} = y (x_{i}) - y_{i} \to 0$ ，则称该方法收敛。

比如，初值问题 ${y^{'} = λ y y (0) = y_{0}$ ，考察欧拉显格式的收敛性：该问题的精确解为 $y (x) = y_{0} e^{λ x}$ ，欧拉公式为 $y_{i + 1} = y_{i} + hλ y_{i} = (1 + λh) y_{i}$ ，对任意固定的 $x_{i} = ih$ ，有 $y_{i} = y_{0} (1 + λh)^{\frac{x _{i}}{h}} = y_{0} [(1 + λh)^{\frac{1}{λh}}]^{λ x_{i}} = 重要极限 y_{0} e^{λ x_{i}} = y (x_{i})$ ，因此欧拉显式法收敛（当 $h \to 0$ 时 $ε_{n} = O (h) \to 0$ ）。

同理，梯形法、欧拉预校法均收敛（当 $h \to 0$ 时 $ε_{n} = O (h^{2}) \to 0$ ）。

定义：用某方法固定步长 $h$ 作计算，由初值 $y_{0}$ 严格得出 $y_{i}$ ，现假设初值有微小误差 $δ_{0}$ （即实际初值为 $y_{0} + δ_{0}$ ），则引起第 $i$ 步计算值有误差 $δ_{i}$ （实际为 $y_{i} + S_{i}$ ），若 $∣ δ_{i} ∣ \leq ∣ δ_{0} ∣$ （ $i = 1, 2, \dots$ ），即 $∣ δ_{i} ∣$ 不随 $i$ 无限扩大，则称该方法关于步长 $h$ 绝对稳定。一般分析时为简单起见，只考虑试验方程 $y^{'} = λ y$ （可以是常数或复数）。

比如，对于欧拉法， $y_{i + 1} = y_{i} + h (λ y_{i}) = (1 + λh) y_{i} = \dots = (1 + λh)^{i + 1} y_{0}$ ，设初值有小扰动 $δ_{0}$ ，则 $y_{i + 1} + δ_{i + 1} = (1 + λh)^{i + 1} (y_{0} + δ_{0})$ ，两式相减得 $δ_{i + 1} = (1 + λh)^{i + 1} δ_{0}$ ，显然，要使欧拉法绝对稳定，就必须保证 $∣1 + λh ∣ \leq 1$ ，即 $h \leq - \frac{1}{λ}$ 时欧拉法绝对稳定。

4.龙格-库塔(R-K)方法

前面介绍的欧拉法的各种变形所能达到的最高精度只有2阶，我们有必要构造更高精度的单步递推格式。

二阶龙格-库塔格式： $⎩ ⎨ ⎧ y_{i + 1} = y_{i} + h (λ_{1} K_{1} + λ_{2} K_{2}) K_{1} = f (x_{i}, y_{i}) K_{2} = f (x_{i} + p h, y_{i} + p h K_{1})$ ，且 ${λ_{1} + λ_{2} = 1 λ_{2} p = \frac{1}{2}$ .

考察改进欧拉法（即欧拉预校法）： $⎩ ⎨ ⎧ y_{i + 1} = y_{i} + \frac{h}{2} (K_{1} + K_{2}) K_{1} = f (x_{i}, y_{i}) K_{2} = f (x_{i} + h, y_{i} + h K_{1}) y_{0} = y (x_{0})$ ，若考虑把 $K_{1}$ 、 $K_{2}$ 改为其他不同的斜率，步长也不再一定是一个 $h$ ，则可将其推广为 $⎩ ⎨ ⎧ y_{i + 1} = y_{i} + h (λ_{1} K_{1} + λ_{2} K_{2}) K_{1} = f (x_{i}, y_{i}) K_{2} = f (x_{i} + p h, y_{i} + p h K_{1})$ ，我们希望能确定系数 $λ_{1}$ 、 $λ_{2}$ 、 $p$ ，使得算法格式有2阶精度，即在 $y_{i} = y (x_{i})$ 的前提下，使得 $R_{i} = y (x_{i + 1}) - y_{i + 1} = O (h^{3})$ ，方法如下：

$K_{2}$ 在 $(x_{i}, y_{i})$ Taylor展开 $K_{2} = f (x_{i} + p h, y_{i} + p h K_{1}) = f (x_{i}, y_{i}) + p h f_{x} (x_{i}, y_{i}) + p h K_{1} f_{y} (x_{i}, y_{i}) + O (h^{2}) = y^{'} (x_{i}) + p h y^{''} (x_{i}) + O (h^{2})$ ，而 $y^{''} (x) = \frac{d}{d x} f (x, y) = f_{x} (x, y) + f_{y} (x, y) \frac{d y}{d x} = f_{x} (x, y) + f_{y} (x, y) f (x, y)$ （其中 $y^{'} = f (x, y)$ 是我们定义的函数），将 $K_{2}$ 代入第一式得 $y_{i + 1} = y_{i} + h {λ_{1} y^{'} (x_{i}) + λ_{2} [y^{'} (x_{i}) + p h y^{''} (x_{i}) + O (h^{2})]} = y_{i} + (λ_{1} + λ_{2}) h y^{'} (x_{i}) + λ_{2} p h^{2} y^{''} (x_{i}) + O (h^{3})$ ，对比Taylor展开 $y (x_{i + 1}) = y (x_{i}) + h y^{'} (x_{i}) + \frac{h ^{2}}{2} y^{''} (x_{i}) + O (h^{3})$ ，若要求 $R_{i} = y (x_{i + 1}) - y_{i + 1} = O (h^{3})$ ，则必须有 $λ_{1} + λ_{2} = 1$ ， $λ_{2} p = \frac{1}{2}$ ，满足这两个条件的有无穷多个解，所有满足上式的格式统称为二阶龙格-库塔格式。

特别地，当 $p = 1$ 、 $λ_{1} = λ_{2} = \frac{1}{2}$ 时就是改进欧拉法。

高阶龙格-库塔格式：在二阶的基础上进行推广，得到 $⎩ ⎨ ⎧ y_{i + 1} = y_{i} + h (λ_{1} K_{1} + λ_{2} K_{2} + \dots + λ_{m} K_{m}) K_{1} = f (x_{i}, y_{i}) K_{2} = f (x_{i} + α_{2} h, y_{i} + β_{21} h K_{1}) K_{3} = f (x_{i} + α_{3} h, y_{i} + β_{31} h K_{1} + β_{32} h K_{2}) \dots K_{m} = f (x_{i} + α_{m} h, y + β_{m 1} h K_{1} + β_{m 2} h K_{2} + \dots + β_{m, m - 1} h K_{m - 1})$ ，其中 $λ_{i}$ 、 $α_{i}$ 与 $β_{ij}$ 均为待定系数，其确定方法和二阶时类似。

比如，三阶龙格-库塔法： $⎩ ⎨ ⎧ y_{i + 1} = y_{i} + \frac{h}{6} (K_{1} + 4 K_{2} + K_{3}) K_{1} = f (x_{i}, y_{i}) K_{2} = f (x_{i} + \frac{h}{2}, y_{i} + \frac{h}{2} K_{1}) K_{3} = f (x_{i} + h, y_{i} - h K_{1} + 2 h K_{2})$ ；四阶龙格-库法： $⎩ ⎨ ⎧ y_{i + 1} = y_{i} + \frac{h}{6} (K_{1} + 2 K_{2} + 2 K_{3} + K_{4}) K_{1} = f (x_{i}, y_{i}) K_{2} = f (x_{i} + \frac{h}{2}, y_{i} + \frac{h}{2} K_{1}) K_{3} = f (x_{i} + \frac{h}{2}, y_{i} + \frac{h}{2} K_{2}) K_{4} = f (x_{i} + h, y_{i} + h K_{3})$ .

龙格-库塔法的计算量与最高精度阶数的关系如下表：

每步须算 $K$ 的个数	2	3	4	5	6	7	$n \geq 8$
可达到的最高精度	$O (h^{2})$	$O (h^{3})$	$O (h^{4})$	$O (h^{4})$	$O (h^{5})$	$O (h^{6})$	$O (h^{n - 2})$

由于龙格-库塔法的导出基于泰勒展开，故精度主要受函数的光滑性影响。对于光滑性不太好的解，最好采用低阶算法而将步长 $h$ 取小。

隐式龙格-库塔法： ${y_{i + 1} = y_{i} + h (λ_{1} K_{1} + \dots + λ_{m} K_{m}) K_{j} = f (x_{i} + α_{j} h, y_{i} + β_{j 1} h K_{1} + \dots + β_{jm} h K_{m})$ ，其中 $j = 1, 2, \dots, m$ 。

对于显式格式的龙格-库塔法，阶数越高，其稳定的区域就越大。而对于隐式龙格-库塔法，其绝对稳定区域远高于显式格式。