科学计算/计算方法备考纲要

科学计算/计算方法 备考纲要

目录

• 科学计算/计算方法 备考纲要
  • 误差
  • 插值与拟合
    • 多项式插值
    • 误差估计
    • 三次样条插值
    • 数据拟合
  • 数值微分与数值积分
    • 数值微分
    • 数值积分的精度
    • Newton-Cotes积分
    • 复合求积与逐次分半
    • 高斯求积公式
  • 线性方程组
    • 消元法
    • 三角分解
    • 误差分析
    • 迭代法——基础迭代格式
    • 迭代法的收敛
  • 非线性方程组
    • 迭代法与收敛速度
    • 加速迭代
    • 非线性方程组牛顿迭代法
  • 常微分方程初值问题
    • 一步解法
    • 误差估计与收敛性
    • 稳定性
    • R-K法
  • 求矩阵的特征值
    • 幂法
    • 反幂法
    • Jacobi方法

声明：本备考纲要以科学计算为主，如有计算方法的相关内容将额外标注。

误差

关于误差的分类，需掌握绝对误差与相对误差。

• 绝对误差：若用\(x^*\)表示\(x\)的准确值的一个近似值，则称\(e^*=x^*-x\)为绝对误差。

  注意：在计算方法中，绝对误差表示为\(e^*(x)=x-x^*\)，与科学计算中的定义正好相反。

• 绝对误差限：若存在某正数\(\varepsilon^*\)使得\(|e^*(x)|<\varepsilon^*\)，则称之为近似值\(x^*\)的绝对误差限。

• 相对误差：记绝对误差与真值的比值为相对误差，即\(e_r(x)=\dfrac{e^*(x)}{x}\)。

• 相对误差限：类似地，若存在\(\varepsilon_r^*\)使\(|e_r(x)|<\varepsilon^*\)，则称之为近似值\(x^*\)的相对误差限。

有效数字：若某近似值\(x^*\)的误差限为其某位数的半个单位，且该位数字到\(x^*\)最左边的非零数字共有\(n\)位，则称\(x^*\)具有\(n\)位有效数字。

误差传播规律：设\(x_1^*,x_2^*\)分别是\(x_1,x_2\)的近似，且\(y=f(x_1,x_2)\)，则用\(y^*=f(x_1^*,x_2^*)\)作为\(y\)的近似时，有

\[y-y^*=f(x_1,x_2)-f(x_1^*,x_2^*)\approx \frac{\partial f(x_1^*,x_2^*) }{\partial x_1}e^*(x_1)+\frac{\partial f(x_1^*,x_2^*)}{\partial x_2}e^*(x_2). \]

事实上，可将相应的误差看作相应的微分，即\(\mathrm{d}x\approx x^*-x\)，\(\mathrm{d}y\approx y^*-y\)，于是

\[\mathrm{d}(x\pm y)=\mathrm{d}x\pm\mathrm{d}y,\\ \mathrm{d}(xy)=y\mathrm{d}x+x\mathrm{d}y,\\ \mathrm{d}\left(\frac{x}{y} \right)=\frac{y\mathrm{d}x-x\mathrm{d}y}{y^2}. \]

此时相对误差为\(\mathrm{d}_r(x)=\left|\dfrac{\mathrm{d}x}{x}\right|=|\mathrm{d}(\ln x)|\)。

例题：设\(x_1^*=0.9863\)，\(x_2^*=0.0062\)分别是\(x_1,x_2\)的近似值，求\(\dfrac{1}{x_1^*}\)的相对误差限以及\(x_1^*x_2^*\)的相对误差限。

解：\(\mathrm{d}x_1=\mathrm{d}x_2=5\times 10^{-5}\)，故

\[\left|\mathrm{d}\left(\ln\frac{1}{x_1} \right)\right|\approx \frac{\mathrm{d}x_1}{x_1^*}=5.0695\times 10^{-5},\\ \left|\mathrm{d}(\ln x_1x_2) \right|\approx \left|\frac{\mathrm{d}x_1}{x_1}+\frac{\mathrm{d}x_2}{x_2} \right|=8.1152\times 10^{-3}. \]

插值与拟合

多项式插值

Lagrange插值的思想为构造基函数，即在某一点处取\(1\)其他点处取\(0\)的函数，最终求和即得到插值多项式。Lagrange插值只适用于最基础的插值构造，当限制条件中有导数时，基函数的形式可能很复杂。

例题：设\(p\in\mathcal P_2\)，且\(p(0)=1\)，\(p(1)=2\)，\(p(2)=1\)，求\(p\)。

解：取

\[l_1(x)=\frac{(x-1)(x-2)}{(0-1)(0-2)}=\frac{1}{2}(x-1)(x-2),\\ l_2(x)=\frac{(x-0)(x-2)}{(1-0)(1-2)}=-x(x-2),\\ l_3(x)=\frac{(x-0)(x-1)}{(2-0)(2-1)}=\frac{1}{2}x(x-1). \]

则

\[p(x)=l_1(x)+2l_2(x)+l_3(x)=-x^2+2x+1. \]

欲构造稍复杂一些的插值多项式，以便处理复杂的条件，一般使用Newton多项式，需要用到差商。一般地，有

\[f[x_i,x_{i+1},\cdots,x_{i+k}]=\frac{f[x_1,\cdots,x_{i+k}]-f[x_0,\cdots,x_{i+k-1}]}{x_{i+k}-x_{i}}. \]

且重结点差商与导数之间具有如下的关系：下设\(k\)为结点重数，则

\[f[x,\cdots,x]=\frac{1}{k!}f^{(k)}(x). \]

例题：设\(p\in\mathcal P_4\)，且

\[p(0)=5,\quad p'(0)=4,\quad p''(0)=6,\quad p(1)=15,\quad p'(1)=20, \]

求\(p\)。

解：列其差商表，有

\[\begin{array}{|c|c|} \hline x & f(x) & \partial_1 & \partial_2 & \partial _3 & \partial _4 \\ \hline 0 & 5 \\ & &4 \\ 0 &5 &&[3]\\ &&4&&3 \\ 0&5&&6&&1 \\ &&10&&4 \\ 1&15&&10 \\ && 20 \\ 1&20\\ \hline \end{array} \]

尤其注意方括号中的\(3\)，它是重结点处的二阶差商，因此应有\(f[0,0,0]=\dfrac{1}{2}p''(0)\)。从而

\[p(x)=5+4x+3x^2+3x^3+x^3(x-1)=x^4+2x^3+3x^2+4x+5. \]

误差估计

代数精度：

以下定义\(\omega_n(x)=\displaystyle{\prod_{i=0}^{n}(x-x_i)}\)。若使用\(\mathcal P_n\)插值来估计一个足够光滑的函数\(f(x)\)，则多项式插值的误差为

\[R_n(x)=f(x)-p_n(x)=\frac{f^{(n+1)}(\xi)}{(n+1)!}\omega_n(x),\quad \xi \in(x_0,x_n). \]

证明即对任何\(x\in[x_0,x_n]\)，构造\(\phi(t)=R(t)-\dfrac{R(x)}{\omega_n(x)}\omega_n(t)\)，它具有至少\(n+2\)个零点，可用\(n+1\)次Rolle定理找到使得\(\phi^{(n+1)}(\xi)=0\)的点。

若使用Hermite插值\(\mathcal P_{2n+1}\)来估计一个足够光滑的函数\(f(x)\)，则插值多项式的误差为

\[R_n(x)=f(x)-p_{2n+1}(x)=\frac{f^{(2n+2)}(\xi)}{(2n+2)!}\omega_n^2(x). \]

类似构造\(\phi(t)=R(t)-\dfrac{R(x)}{\omega_n^2(x)}\omega_n^2(t)\)，这里\(\phi'(t)\)至少具有\(2n+1\)个零点，因此对\(\phi'(t)\)用\(2n+1\)次Rolle定理使得找到\(\phi^{(2n+2)}(\xi)=0\)的点。

始终记住，求导次数等于条件数。

三次样条插值

本部分只针对计算方法，因科学计算对样条插值的计算不作要求。

三次样条插值需要增加两个边界条件，常见的有：

1. 端点处一阶导数值：\(S'(x_0)\)和\(S'(x_n)\)已知。
2. 端点处二阶导数值：\(S''(x_0)\)和\(S''(x_n)\)已知。
3. 周期条件：\(S'(x_0)=S'(x_n)\)，\(S''(x_0)=S''(x_n)\)。

为计算，假设\(S(x)\)在每个结点处有\(S''(x_i)=M_i\)。在每个小区间\([x_{i-1},x_i]\)上\(S(x)=S_i(x)\)是不高于三次的多项式，故其二阶导数不超过一次，可构造插值基函数为

\[S_i''(x)=M_{i-1}\frac{x-x_i}{x_{i-1}-x_i}+M_i\frac{x-x_{i-1}}{x_i-x_{i-1}}\xlongequal{def}\frac{M_{i-1}}{h_{i}}(x_i-x)+\frac{M_i}{h_i}(x-x_{i-1}). \]

称\(h_i=x_i-x_{i-1}\)为步长。利用基函数思想对此函数二次积分（自行体会系数的配凑），得到

\[\begin{aligned} S_i(x)&=\frac{M_{i-1}(x_i-x)^3}{6h_i}+\frac{M_i(x-x_{i-1})^3}{6h_i}\\ &+\left(y_{i-1}-\frac{M_{i-1}h_i^2}{6} \right)\frac{x_i-x}{h_i}+\left(y_i-\frac{M_ih_i^2}{6} \right)\frac{x-x_{i-1}}{h_i}. \end{aligned} \]

这样，每一个方程都被一个\(M_i\)唯一确定，为求\(M_i\)，需用一阶导数条件，即\(S_i'(x_{i})=S_{i+1}'(x_i)\)，所以

\[S_i'(x_i)=\frac{M_ih_i}{2}+\frac{y_i-y_{i-1}}{h_i}-\frac{h_i(M_i-M_{i-1})}{6}=\frac{y_i-y_{i-1}}{h_i}+\frac{h_iM_{i-1}}{6}+\frac{h_iM_i}{3};\\ S'_{i+1}(x_i)=-\frac{M_{i}h_{i+1}}{2}+\frac{y_{i+1}-y_{i}}{h_{i+1}}-\frac{h_{i+1}(M_{i+1}-M_{i})}{6}=\frac{y_{i+1}-y_{i}}{h_{i+1}}-\frac{h_{i+1}M_{i}}{3}-\frac{h_{i+1}M_{i+1}}{6}. \]

利用它们相等，\(i=1,2,\cdots,n-1\)，可得

\[\frac{h_i}{6}M_{i-1}+\frac{h_i+h_{i+1}}{3}M_i+\frac{h_{i+1}}{6}M_{i+1}=\frac{y_{i+1}-y_{i}}{h_{i+1}}-\frac{y_i-y_{i-1}}{h_i}.\\ \]

为求解，需将边界条件予以应用。最简单的边界条件是给出\(S''(x_0)\)和\(S''(x_n)\)时的边界，此时\(M_0\)和\(M_n\)直接给出；而如果给出\(S'(x_0)\)和\(S'(x_n)\)，则要对\(S_1(x)\)和\(S_n(x)\)求导，得到

\[-\frac{h_1M_0}{2}+\frac{y_1-y_0}{h_1}+\frac{h_1(M_1-M_0)}{6}=0,\\ \frac{h_nM_n}{2}+\frac{y_{n}-y_{n-1}}{h_n}+\frac{h_n(M_n-M_{n-1})}{6}=0. \]

数据拟合

拟合即使用最小二乘法计算，由于正交多项式不作要求，因此只考察基础最小二乘法，至多扩展至多项式拟合。回归方程为\(y=(a_0,a_1,\cdots,a_k)(1,x,\cdots,x^k)'\)，考察残差平方和为

\[e'e=\sum_{i=1}^{n}\left(y_i-\sum_{j=0}^{k}a_jx^{j} \right)^2, \]

对\(k+1\)个参数分别求偏导，就得到法方程组为

\[\sum_{i=1}^{n}x^{l}\left(y_i-\sum_{j=0}^{k}a_jx^{j} \right)=0,\quad l=0,1,\cdots,k. \]

由此可以解出最优统计量。特别当回归方程为直线时，有

\[\hat a_1=\frac{\sum(x_i-\bar{x})(y_i-\bar{y})}{\sum(x_i-\bar{x})^2},\quad \hat a_0=\bar{y}-\hat a_1\bar{x}. \]

这是回归分析的基础结果。

例：给定以下超定方程组

\[2x+3y=8,\\ 3x+6y=6,\\ x+y=5. \]

求最小二乘解。

解：正规方程为

\[39-14x-25y=0,\\ 65-25x-46y=0, \]

解得

\[\hat x=8.894737,\\ \hat{y}=-3.421053. \]

数值微分与数值积分

数值微分

我们所用的数值微分基于插值多项式构造，即用插值多项式的微商来替代函数的微商，且一般用等距节点构造的插值多项式求微分，常用的有两点式和三点式。对于两点式，有

\[\varphi(x)=\frac{x_1-x}{h}y_0+\frac{x-x_0}{h}y_1,\\ \varphi'(x)=\frac{y_1-y_0}{h}, \]

事实上就是两点之间的斜率。对于三点式，有

\[\varphi(x)=\frac{(x-x_1)(x-x_2)}{2h^2}y_0+\frac{(x-x_0)(x-x_2)}{-h^2}y_1+\frac{(x-x_0)(x-x_1)}{2h^2}y_2,\\ \varphi'(x)=\frac{(x-x_1)+(x-x_2)}{2h^2}y_0-\frac{(x-x_0)(x-x_2)}{h^2}y_1+\frac{(x-x_0)+(x-x_1)}{2h^2}y_2,\\ \varphi'(x_0)=\frac{-3y_0+4y_1-y_2}{2h},\\ \varphi'(x_1)=\frac{y_2-y_0}{2h},\\ \varphi'(x_2)=\frac{y_0-4y_1+3y_2}{2h}. \]

对其误差估计，可以直接由插值多项式的误差估计推得，有

\[f(x)-\varphi(x)=\frac{f^{(n+1)}(\xi)}{(n+1)!}\omega_n(x),\\ f'(x)-\varphi'(x)=\frac{\omega_n'(x)}{(n+1)!}f^{(n+1)}(\xi)+\frac{\mathrm{d}f^{(n+1)}(\xi(x))}{(n+1)!}\omega_n(x). \]

注意\(f^{(n+1)}(\xi)\)的导数一般无法求得，但若在插值结点处，则由\(\omega_n(x_i)=0\)可以忽略这一项，此时

\[f'(x_i)-\varphi'(x_i)=\frac{f^{(n+1)}(\xi)}{(n+1)!}\frac{\omega_n(x)}{(x-x_i)}. \]

对两点式，有

\[R_1(x_0)=-\frac{hf^{(2)}(\xi)}{2}=O(h),\\ R_1(x_1)=\frac{hf^{(2)}(\xi)}{2}=O(h). \]

对三点式，则可以类似推导，得到其截断误差为\(O(h^2)\)。

例：对三点估计的\(\varphi''(x)\)，证明\(\varphi''(x_1)\)的截断误差是\(O(h^2)\)。

解：由

\[f(x)-\varphi(x)=\frac{f^{(3)}(\xi)}{3!}(x-x_0)(x-x_1)(x-x_2), \]

有

\[\begin{aligned} &\quad f'(x)-\varphi'(x)\\ &=\frac{f^{(3)}(\xi)}{6}\left[(x-x_1)(x-x_2)+(x-x_0)(x-x_2)+(x-x_0)(x-x_1) \right]\\ & \quad +\frac{\omega_2(x)}{6}\frac{\mathrm{d}f^{(3)}(\xi)}{\mathrm{d}x},\\ &\quad f''(x_1)-\varphi''(x_1)\\ &=\frac{f^{(3)}(\xi)}{6}[-h+h-h+h]+\frac{-h^2}{6}\frac{\mathrm{d}f^{(3)}(\xi)}{\mathrm{d}x}+\frac{-h^2}{6}\frac{\mathrm{d}f^{(3)}(\xi)}{\mathrm{d}x}\\ &=-\frac{h^2}{6}f^{(4)}(\eta)\\ &=O(h^2). \end{aligned} \]

数值积分的精度

在区间\([a,b]\)上，不同的求积公式表现为

\[\int_{a}^{b}f(x)\mathrm{d}x\stackrel{I}=\sum_{k=0}^{n}A_kf(x_k). \]

这里\(A_k\)是不依赖于\(f(x)\)的常数。对正常积分区间，一般将\(x\in[a,b]\)作变换\(t=\dfrac{x-a}{b-a}\in[0,1]\)，这样对求积公式的计算将变得简单；而如果是反常积分，则一般要用一定的换元手段将积分换成正常区间，且使得换元后的积分在区间内有界（即非瑕积分）。

为衡量不同求积公式的表现，引入代数精度概念。

• 代数精度：若求积公式对任意一个次数不高于\(m\)次的多项式\(f(x)\)都能精确求积，而当\(f(x)\)为\(m+1\)阶多项式时不精确成立，则求积公式具有\(m\)次代数精度。

为验证求积公式具有\(m\)阶代数精度，一般对\(f(x)=1,x,\cdots,x^{m}\)分别验证，因为\(\mathcal P_m\)构成一个线性空间，所以只需取其一组基。

例：设\(f\in C^{4}[a,b]\)，求积公式为

\[\int_{a}^{b}f(x)\mathrm{d}x\approx \frac{b-a}{8}\left[f(a)+3f\left(\frac{2a+b}{3} \right)+3f\left(\frac{a+2b}{3} \right)+f(b) \right]. \]

求该求积公式的代数精度。

解：如果该公式在\([0,1]\)上成立，则由换元，有

\[\begin{aligned} \int_{a}^{b}f(x)\mathrm{d}x&=(b-a)\int_{0}^{1}f[(b-a)t+a]\mathrm{d}t\\ &=\frac{b-a}{8}\left[f(a)+3f\left(\frac{2a+b}{3} \right)+3f\left(\frac{a+2b}{3} \right)+f(b) \right]. \end{aligned} \]

代入\(f(x)=1,x,x^2,x^3\)，有

\[\int_{0}^{1}1\mathrm{d}x=\frac{1}{8}\cdot 8=1,\\ \int_{0}^{1}x\mathrm{d}x=\frac{1}{8}\left[0+3\frac{1}{3}+3\frac{2}{3}+1 \right]=\frac{1}{2},\\ \int_{0}^{1}x^2\mathrm{d}x=\frac{1}{8}\left[0+3\frac{1}{9}+3\frac{4}{9}+1 \right]=\frac{1}{3},\\ \int_{0}^{1}x^{3}\mathrm{d}x=\frac{1}{8}\left[0+3\frac{1}{27}+3\frac{8}{27}+1 \right]=\frac{1}{4}. \]

但当\(f(x)=x^4\)时，

\[\int_{0}^{1}x^{4}\mathrm{d}x\ne \frac{1}{8}\left[0+3\frac{1}{81}+3\frac{16}{81}+1 \right]=\frac{11}{54}. \]

所以该求积公式具有\(3\)阶代数精度。

Newton-Cotes积分

Newton-Cotes积分公式基于等距插值多项式构造，即对插值多项式积分来近似原积分。现讨论\([0,1]\)区间，在结点等距的情况下，诸\(A_k\)分别为

\[\begin{array}{c|c} \hline n & A_k \\ \hline 1 & \frac{1}{2}\quad \frac{1}{2} \\ 2 & \frac{1}{6} \quad \frac{4}{6} \quad \frac{1}{6} \\ 3 & \frac{1}{8} \quad \frac{3}{8} \quad \frac{3}{8} \quad \frac{1}{8} \\ \hline \end{array} \]

一般只会用到至多四个结点的Newton-Cotes公式，其中尤以\(n=1\)，\(n=2\)应用广泛，分别称为梯形公式与抛物线公式。可以证明：

• \(n\)为奇数时，Newton-Cotes公式具有\(n\)次代数精度。
• \(n\)为偶数时，Newton-Cotes公式具有\(n+1\)次代数精度。

对\(n\)为奇数时的证明比较容易，因\(f\in\mathcal P_n\)时\(f^{(n+1)}(\eta)=0\)，所以

\[R(x)=\int_{1}^{0}\frac{f^{(n+1)}(\xi(x))}{(n+1)!}\omega_n(x)\mathrm{d}x=0. \]

但当\(f\in\mathcal P_{n+1}\)时容易计算\(\displaystyle{\int_{0}^{1}\omega_n(x)\mathrm{d}x\ne 0}\)。

\(n\)为偶数时，一个关键的结果是\(\displaystyle{\int_{0}^{1}\omega_n(x)\mathrm{d}x=0}\)，因可通过换元将其变为\([-1,1]\)上奇函数的积分。

对梯形公式和抛物线公式，下面计算其误差估计，先给出积分中值定理：对闭区间\([a,b]\)上连续的函数\(f(x),g(x)\)，若\(g(x)\)在\([a,b]\)上不变号，则\(\displaystyle{\int_{a}^{b}f(x)g(x)\mathrm{d}x=f(\xi)\int_{a}^{b}g(x)\mathrm{d}x}\)。

现在通过插值的余项估计，我们能给出梯形公式与求积公式的误差估计。首先是梯形公式，因\(\omega(x)=x(x-1)\)在\([0,1]\)上不变号，所以可以用积分中值定理，得到

\[R(f)=\int_{0}^{1}\frac{f^{(2)}(\xi(x))}{2}x(x-1)\mathrm{d}x=-\frac{f^{(2)}(\eta)}{12}. \]

而对抛物线公式则有一个小技巧，因为它具有\(4\)阶代数精度，因此可以构造一个\(4\)次多项式来作为插值多项式，不妨选择\(x=\dfrac{1}{2}\)处的重结点，这样做的好处是\(x(x-\dfrac{1}{2})^2(x-1)\)会在\([0,1]\)上不变号，从而可以使用积分中值定理，有

\[R(f)=\int_{0}^{1}\frac{f^{(4)}(\xi(x))}{4!}x\left(x-\frac{1}{2} \right)^2(x-1)\mathrm{d}x=-\frac{f^{(5)}(\eta)}{24\times 120}=-\frac{f^{(4)}(\eta)}{2880}. \]

需要注意，如果是在\([a,b]\)上作积分，误差估计需要乘上\((b-a)^{n+2}\)，这里\(n\)是方法的代数精度。通过积分中值定理的做法，只能处理梯形公式与抛物线公式的误差估计，对于更多结点数的Cotes系数，则需要更高级的方法。

复合求积与逐次分半

复合求积指的是在小区间上作积分，常常在小区间上使用梯形求积公式或者抛物线求积公式。这里给出复合梯形、抛物线求积公式的误差估计（\(n\)均代表复合求积的复合区间数）：

\[R(f,T_n)=-\frac{b-a}{12}h^2f''(\eta),\quad h=\frac{b-a}{n},\\ R(f,S_n)=-\frac{b-a}{2880}h_1^{4}f^{(4)}(\eta)=-\frac{b-a}{180}h^4f^{(4)}(\eta),\quad h=\frac{b-a}{2n},h_1=2h. \]

为了保证计算速度的同时保证计算精度，往往不用事前估计来确定需要使用的结点数，而是使用逐次分半法。我们知道\(h\)较小时\(R(f,T_n)=O(h^2)\)，因此当\(h\)减半时，\(R(f,T_{2n})\approx R(f,T_n)/4\)，也就是\(I-R(f,T_n)\approx 4[I-R(f,T_{2n})]\)，于是

\[T_{2n}-T_{n}=(I-T_{2n})-(1-T_{n})=R(f,T_{2n})-R(f,T_{n})=3[I-R(f,T_{2n})]=3\varepsilon,\\ \]

对梯形公式，其逐次分半法表现为

\[T_{2^0}=\frac{b-a}{2}\left[f(a)+f(b) \right];\\ T_{2^1}=\frac{T_{2^{0}}}{2}+\frac{b-a}{2}\left[f\left(\frac{a+b}{2} \right) \right];\\ T_{2^2}=\frac{T_{2^{1}}}{2}+\frac{b-a}{4}\left[f\left(\frac{3a+b}{4} \right)+f\left(\frac{a+3b}{4} \right) \right];\\ \cdots \\ T_{2^{k}}=\frac{T_{2^{k-1}}}{2}+\frac{b-a}{2^{k}}\sum_{i=1}^{2^{k-1}}f\left(a+(2i-1)\frac{b-a}{2^{k}} \right). \]

这样，\(T_{2^{k}}\)序列将趋近于积分值\(I\)，当\(|T_{2^{k}}-T_{2^{k-1}}|<3\varepsilon\)就停止分半。

下面介绍Romberg算法。在前面的推导中，用\(T_{2n}\)作为积分的近似时，其误差约为\(\dfrac{1}{3}(T_{2n}-T_{n})\)，所以如果将这个值作为补偿，得到一个新的序列，就有

\[\tilde{T}_n=T_{2n}+\frac{1}{3}(T_{2n}-T_{n})=\frac{4}{3}T_{2n}-\frac{1}{3}T_{n}. \]

事实上，我们有\(\tilde{T}_n=S_n\)，然而利用同样的外推手段，可以获得

\[\tilde{S}_n=\frac{16}{15}S_{2n}-\frac{1}{15}S_{n}={C}_n,\\ \tilde{C}_n=\frac{64}{63}C_{2n}-\frac{1}{63}C_n=R_n,\\ \cdots \]

这就是Romberg外推，这里每一个\(\{T_n\}\)、\(\{S_n\}\)、\(\{C_n\}\)、\(\{R_n\}\)都将收敛于\(I\)。

高斯求积公式

Gauss求积公式相对于Cotes公式就没有那么好记，它的构造基于一个非常自然的想法：最大化代数精度。因为\(\displaystyle{\sum_{k=0}^{n}f(x_k)A_k}\)实际上有\(2n+2\)个参数，所以它理论上应当能具有\(2n+1\)阶代数精度，这远比Cotes公式提供的\(n\)或\(n+1\)阶代数精度高。

为了凑出这样的\(x_k\)与\(A_k\)，先将积分区域变换到\([-1,1]\)，从而要使\(2n+2\)个方程精确成立：

\[\sum_{k=0}^{n}A_kf(x_k)=\int_{-1}^{1}x^{l}\mathrm{d}x,\quad l=0,1,\cdots,2n+1. \]

然而，由于此时的\(\{x_k\}\)是一组方程的解，所以它不再是等距的，且没有规律可循。因此，Gauss公式虽然具有较高代数精度，但形式繁琐，没有Cotes公式简洁。

线性方程组

消元法

解线性方程组\(Ax=b\)，最原始的方法是高斯消元法，即通过初等行变换将系数矩阵变成上三角的。实际操作中，视初始矩阵为\(A^{(1)}\)，先将第一行乘以\(-\dfrac{a_{i}^{(1)}}{a_1^{(1)}}\)加到第\(i=2:n\)行，得到新的系数矩阵\(A^{(2)}\)和\(b^{(2)}\)；再将第二行乘以\(-\dfrac{a_{i}^{(2)}}{a_2^{(2)}}\)加到第\(i=3:n\)行，得到新的系数矩阵\(A^{(3)}\)和\(b^{(3)}\)。以此类推，一共需作\(n-1\)次乘法，记最后的方程组为\(A^{(n)}x=b^{(n)}\)，它与原方程组通解，且\(A^{(n)}\)是上三角的。

用矩阵的观点看高斯消元，并记\(m_{i,j}=\dfrac{a_i^{(j)}}{a_j^{(j)}}(i>j)\)，那么每一步操作相当于用变换矩阵\(M_i\)来左乘\(A^{(i)}\)（行变换为左乘），并且

\[M_i=\begin{bmatrix} 1 \\ \vdots & \ddots \\ 0 & \cdots & 1 \\ 0 & \cdots & -m_{i+1,i} & \ddots \\ \vdots & & \vdots & 0& \ddots \\ 0 & \cdots & -m _{n,i} & 0 & \cdots& 1 \end{bmatrix},\quad \begin{array}{} A^{(n)}=M_{n-1}M_{n-2}\cdots M_1A^{(1)},\\ b^{(n)}=M_{n-1}M_{n-2}\cdots M_1b^{(1)}. \end{array} \]

也就是\(M_i\)只有第\(i\)列有元素，主对角线上元素为\(1\)，是一个多重倍加阵，由倍加的性质显然

\[M^{-1}=\begin{bmatrix} 1 \\ \vdots & \ddots \\ 0 & \cdots & 1 \\ 0 & \cdots & m_{i+1,i} & \ddots \\ \vdots & & \vdots & 0& \ddots \\ 0 & \cdots & m _{n,i} & 0 & \cdots& 1 \end{bmatrix}, \]

即反向倍加即可。并且我们有

\[A^{(1)}=M_1^{-1}M_2^{-1}\cdots M_{n-1}^{-1}A^{(n)}. \]

显然\(\det(A^{(1)})=\det(A^{(n)})\)。

为使计算精度增大，常使用列主元消去法，它基于除数不能太小的想法，每次将当前主元位置中，系数最大的那个方程移到当前行。如果当前主元是\(x_i\)，\(A^{(i)}\)中第\(i\)列系数最大的是\(a_{k,i}\)，则需将第\(k\)行与第\(i\)行互换。在矩阵的操作中，相当于在\(M_i\)处理之前，先用对换矩阵\(P_i\)处理一遍，也就是

\[A^{(n)}=M_{i-1}P_{i-1}M_{i-2}P_{i-2}\cdots M_1P_1A^{(1)},\\ b^{(n)}=M_{i-1}P_{i-1}M_{i-2}P_{i-2}\cdots M_1P_1b^{(1)}. \]

三角分解

三角分解基于如下的想法：如果线性方程组的系数矩阵是上三角或者下三角的，则可以分别通过前代与后代，递推地给出方程的解。将一个矩阵分成两个三角矩阵的乘积，就能通过简单的代入求解方程，而\(LU\)分解就是指将矩阵分解成下三角阵\(L\)与上三角阵\(U\)的乘积。

\(LU\)分解中较常用的是杜利特尔分解，将矩阵分成主对角线为\(1\)的下三角矩阵，与一般的上三角矩阵的乘积。事实上，在高斯消元法中，我们有\(A^{(1)}=M_1^{-1}M_2^{-1}\cdots M_{n-1}^{-1}A^{(n)}\)，如令\(L=M_1^{-1}\cdots M_{n-1}^{-1}\)，则\(L\)是主对角线为\(1\)的下三角阵，\(A^{(n)}\)是一般的上三角阵，已经达到目的。这说明顺序主子式不为\(0\)的矩阵均可以进行\(LU\)分解。

现在我们直观地从代数上进行杜利特尔\(LU\)分解，设

\[L=\begin{pmatrix} 1 & \\ l_{2,1} & 1 \\ l_{3,1} & l_{3,2} & 1 \\ \vdots & \vdots & \vdots & \ddots \\ l_{n,1} & l_{n,2} & l_{n,3} & \cdots & 1 \end{pmatrix},\quad U=\begin{pmatrix} u_{1,1} & u_{1,2} & u_{1,3} & \cdots & u_{1,n} \\ & u_{2,2} & u_{2,3} & \cdots & u_{2,n} \\ && \ddots & & \vdots \\ &&& u_{n-1,n-1} & n_{n-1,n} \\ &&&& u_{n,n} \end{pmatrix}. \]

那么计算过程中，我们有

\[u_{1,i}=a_{1,i},\\ l_{i,1}=\frac{a_{i,1}}{u_{1,1}};\\ u_{2,i}=a_{2,i}-l_{2,1}u_{1,i},\\ l_{i,2}=\frac{a_{i,2}-l_{i,1}u_{1,2}}{u_{2,2}}. \]

随着\(k\)增大，逐层先\(u_{k,i}\)后\(l_{i,k}\)计算下去，最终有

\[u_{k,i}=a_{k,i}-\sum_{j=1}^{k-1}l_{k,j}u_{j,i},\\ l_{i,k}=\frac{a_{i,k}-\displaystyle{\sum_{j=1}^{k-1}l_{i,j}u_{j,k}}}{u_{k,k}}. \]

然后，方程变为\(Ax=LUx=b\)，令\(Ux=y\)，可以先从\(Ly=b\)前代解出\(y\)，再从\(Ux=y\)后代解出\(x\)。

从\(LU\)分解可以推广出\(LDR\)分解，它使\(L\)是单位下三角阵，\(R\)是单位上三角阵。因\(LU\)分解中的\(L\)已经是单位下三角的，所以只要把\(U\)分解成单位上三角的，即每一行乘以其主对角线元素的倒数即可，乘积矩阵记作\(D\)，它是对角阵。

例：\(LU\)分解，求解

\[\begin{pmatrix} 3 & 2 & 1 \\ 2 & 4 & 1 \\ 1 & 2 & 4 \end{pmatrix}\begin{pmatrix} x_1 \\ x_2 \\ x_3 \end{pmatrix}=\begin{pmatrix} 2 \\ -1 \\ 3 \end{pmatrix}. \]

解：对系数矩阵进行\(LU\)分解，有

\[u_{1,1}=3,\quad u_{1,2}=2, \quad u_{1,3}=1,\\ l_{2,1}=\frac{2}{3},\quad l_{3,1}=\frac{1}{3},\\ u_{2,2}=\frac{8}{3},\quad u_{2,3}=\frac{1}{3},\\ l_{3,2}=\frac{1}{2},\\ u_{3,3}=\frac{7}{2}. \]

即

\[L=\begin{pmatrix} 1 \\ \frac{2}{3} & 1 \\ \frac{1}{3} & \frac{1}{2} & 1 \end{pmatrix},\quad U=\begin{pmatrix} 3 & 2 & 1 \\ & \frac{8}{3} & \frac{1}{3} \\ && \frac{7}{2} \end{pmatrix},\quad A=LU. \]

先解\(Ly=b\)，得到\(y=(2,-7/3,7/2)'\)，再解\(Ux=y\)，得\(x=(-1,-1,1)'\)。

对对称正定阵\(A\)，存在一个实的非奇异下三角阵\(L\)使得\(A=LL'\)，从而将\(A\)分成下三角矩阵与上三角矩阵的乘积，这种分解称为Cholesky分解。从代数上看Cholesky分解，有

\[L=\begin{pmatrix} l_{1,1} \\ l_{1,2} & l_{2,2} \\ \vdots & \vdots & \ddots \\ l_{1,n} & l_{2,n} & \cdots & l_{n,n} \end{pmatrix},\quad L'=\begin{pmatrix} l_{1,1} & l_{1,2} & \cdots & l_{1,n} \\ & l_{2,2} & \cdots & l_{2,n} \\ && \ddots & \vdots \\ &&& l_{n,n} \end{pmatrix}, \]

分别计算，就有

\[l_{1,1}=\sqrt{a_{1,1}},\\ l_{1,i}=\frac{a_{i,1}}{l_{1,1}},\\ l_{2,2}=\sqrt{a_{2,2}-l_{1,2}^2},\\ l_{2,i}=\frac{a_{2,i}-l_{1,1}l_{1,i}}{l_{2,2}}, \]

类推计算，可得

\[l_{i,i}=\sqrt{a_{i,i}-\displaystyle{\sum_{j=1}^{i-1}a_{j,i}^2}},\\ l_{i,k}=\frac{a_{i,k}-\displaystyle{\sum_{j=1}^{j-1}l_{i,j}l_{j,k}}}{l_{i,i}}. \]

为避免开方运算，可以使用\(LDL'\)分解，即额外用一个对角矩阵来存储一个需要开方的结果。

误差分析

先给出向量范数的定义：

• 向量范数：满足以下三个性质的映射称为向量范数。
  1. 正定性：\(\forall x\in\mathbb{R}^{n}\)，\(\|x\|\ge 0\)，等号成立当期仅当\(x=0\)。
  2. 齐次性：\(\forall x\in\mathbb{R}^{n}\)，\(\lambda\in\mathbb{R}\)，成立\(\|\lambda x\|=\lambda\|x\|\)。
  3. 三角不等式：\(\forall x,y\in\mathbb{R}\)，成立\(\|x+y\|\le \|x\|+\|y\|\)。

常用的向量范数一般是Minkowski范数如下：

• \(1\)范数：\(\|x\|_1=\displaystyle{\sum_{i=1}^{n}|x_i|}\)。
• \(2\)范数：\(\|x\|_2=\displaystyle{\sqrt{\sum_{i=1}^{n}x_i^2}}\)。
• \(\infty\)范数：\(\|x\|_\infty=\displaystyle{\lim_{p\to \infty}\left(\sum_{i=1}^{n}x_i^{p} \right)^{1/p}=\max_i|x_i|}\)。

根据向量范数，可以推广出矩阵范数，它也要满足正定性、齐次性与三角不等式三个性质，还需要满足相容性：\(\|AB\|\le \|A\|\cdot \|B\|\)，\(\|Ax\|\le \|A\|\cdot \|x\|\)。常用的矩阵范数有以下几种：

• \(1\)范数：\(\displaystyle{\|A\|_1=\max_{\|x\|_1=1}\|Ax\|_1=\max_{1\le j\le n}\sum_{i=1}^{n}|a_{ij}|}\)，即各列元素绝对值之和中的最大者。

  取\(x=e_i\)，分别计算\(\|Ae_i\|_1\)即可推知。

• \(2\)范数：\(\|A\|_2=\displaystyle{\max_{\|x\|_2=1}\|Ax\|_2=\sqrt{\lambda_1}}\)，\(\lambda_1\)为矩阵\(A'A\)的最大特征值。

• \(\infty\)范数：\(\|A\|_{\infty}=\displaystyle{\max_{\|x\|_{\infty}=1}\|x\|_\infty=\max_{1\le i\le n}\sum_{j=1}^{n}|a_{ij}|}\)，即各行元素绝对值中最大者。

  取\(x=(1/n,\cdots,1/n)\)即可验证。

• \(F\)范数：\(\|A\|_{r}=\displaystyle{\sqrt{\sum_{i,j=1}^{n}|a_{ij}|^2}}\)。

下给出矩阵谱半径与条件数的定义，从而用它进行误差分析。

• 谱半径：设\(A\)的特征值为\(\lambda_i(i=1,2,\cdots,n)\)，则称矩阵\(A\)的谱半径为

  \[\rho(A)=\max_{1\le i\le n}|\lambda_i|. \]

• 条件数：称矩阵\(A\)的条件数为

  \[\mathrm{Cond}(A)=\|A\|\cdot\|A^{-1}\|. \]

  条件数与范数有关，其中最常用的为谱条件数，即\(\mathrm{Cond}(A)_2=\|A\|_2\cdot\|A^{-1}\|_2=\sqrt{\lambda_1/\lambda_n}\)，这里\(\lambda_1,\lambda_n\)分别是\(A'A\)的最大、最小特征值。可以验证有

  1. \(\forall A\)，\(\mathrm{Cond}(A)\ge 1\)。
  2. 单位阵，置换阵，正交阵的谱条件数\(\mathrm{Cond}_2=1\)。
  3. 对任意非零实数\(\alpha\)，有\(\mathrm{Cond}(A)=\mathrm{Cond}(\alpha A)\)。

对于某些线性方程组\(Ax=b\)，对\(A\)和\(b\)进行微小扰动将极大影响解的精确度，事实上这主要由\(A\)影响。若\(b\)受到微小扰动，则

\[A(x+\Delta x)=b+\Delta b\Rightarrow A\Delta x=\Delta b, \]

所以

\[\Delta x=A^{-1}\Delta b,\\ \|\Delta x\|\le \|A^{-1}\Delta b\|\le \|A^{-1}\|\cdot \|\Delta b\|,\\ \|\Delta x\|\cdot\|b\|\le \|A^{-1}\|\cdot\|A\|\cdot\|x\|\cdot\|\Delta b\|, \]

即相对误差是

\[\frac{\|\Delta x\|}{\|x\|}\le \mathrm{Cond}(A)\frac{\|\Delta b\|}{\|b\|}. \]

对\(A\)发生扰动也可以证明类似结果，这说明如果矩阵的\(\mathrm{Cond}(A)\)很大，则方程组为病态方程组。

例：设

\[A=\begin{bmatrix} 3 & 1.001 \\ 6 & 1.997 \end{bmatrix}, \]

求\(A\)的谱条件数。

解：计算得

\[A'A=\begin{bmatrix} 3 & 6 \\ 1.001 & 1.997 \end{bmatrix}\begin{bmatrix} 3 & 1.001 \\ 6 & 1.997 \end{bmatrix}=\begin{bmatrix} 45 & 14.985 \\ 14.985 & 4.99001 \end{bmatrix}, \]

有

\[\lambda_1=49.99,\quad \lambda_2=4.5009\times 10^{-6}, \]

所以

\[\mathrm{Cond}(A)_2=\sqrt{\lambda_1/\lambda_2}\approx3333. \]

迭代法——基础迭代格式

迭代法即对\(Ax=b\)，构造某个等价迭代序列\(x=Bx+f\)，如果\(\{x_k\}\)收敛，则其极限必为\(x\)。需掌握Jacobi迭代法与Gauss-Seidel迭代法。

Jacobi迭代法相对最为直观，对方程组

\[\begin{pmatrix} a_{1,1} & a_{1,2} & \cdots & a_{1,n} \\ a_{2,1} & a_{2,2} & \cdots & a_{2,n} \\ \vdots & \vdots & & \vdots \\ a_{n,1} & a_{n,2} & \cdots & a_{n,n} \end{pmatrix}\begin{pmatrix} x_1 \\ x_2 \\ \vdots \\x_n \end{pmatrix}=\begin{pmatrix} b_1 \\ b_2 \\ \vdots \\ b_n \end{pmatrix}, \]

先对角元素归一化，即第\(i\)行除以\(a_{i,i}\)，得到

\[\begin{pmatrix} 1 & -b_{1,2} & \cdots & -b_{1,n} \\ -b_{2,1} & 1 & \cdots & -b_{2,n} \\ \vdots & \vdots & & \vdots \\ -b_{n,1} & -b_{n,2} & \cdots & 1 \end{pmatrix}\begin{pmatrix} x_1 \\ x_2 \\ \vdots \\x_n \end{pmatrix}=\begin{pmatrix} g_1 \\ g_2 \\ \vdots \\ g_n \end{pmatrix},\\ b_{i,j}=-\frac{a_{i,j}}{a_{i,i}},\quad g_i=\frac{b_i}{a_{i,i}}. \]

此时的系数矩阵出现了一个单位阵\(I\)，由于\(Ix=x\)，可以构造如此的迭代序列：

\[\begin{pmatrix} x_1 \\ x_2 \\ \vdots \\ x_n \end{pmatrix}=\begin{pmatrix} 0 & b_{1,2} & \cdots & b_{1,n} \\ b_{2,1} & 0 & \cdots & b_{2,n} \\ \vdots & \vdots & & \vdots \\ b_{n,1} & b_{n,2} & \cdots & b_{n,n} \end{pmatrix}\begin{pmatrix} x_1 \\ x_2 \\ \vdots \\ x_n \end{pmatrix}+\begin{pmatrix} g_1 \\ g_2 \\ \vdots \\ g_n \end{pmatrix},\\ x=B_{J}x+f_{J}. \]

从上面的操作中，我们可以提取\(A\)的对角部分为\(D\)，这样，第一步我们用\(D^{-1}\)左乘\(Ax+b\)的两边，第二步从\(D^{-1}A\)中分离了一个\(I\)，并将其移到等式的另一边，因此Jacobi迭代的步骤是

\[Ax=b\\ D^{-1}Ax=D^{-1}b,\\ x=(I-D^{-1}A)x+D^{-1}b\xlongequal{def}B_{J}x+f_{J}. \]

同时，可以写出其分量形式为

\[x_i^{(k+1)}=\sum_{j\ne i}^{n}b_{i,j}x_{j}^{(k)}+g_{i}=-\sum_{j\ne i}^{n}\frac{a_{i,j}}{a_{i,i}}+\frac{b_{i}}{a_{i,i}}. \]

欲理解Gauss-Seidel迭代法，应从Jacobi迭代法的分量形式看起。对Jacobi迭代，\(n\)个分量是用上一次的数值\(x^{(k)}\)同时迭代的，这样需要额外的\(n\)个存储空间；而如果每次计算完\(x^{(k+1)}_1\)，就将\(x^{(k+1)}_1\)的结果用于\(x_2^{(k+1)}\)的计算，也就是用\(x_1^{(k+1)},\cdots,x_i^{(k+1)}\)的结果计算\(x_{i+1}^{(k+1)}\)，这时候，分量形式就是

\[\left\{\begin{array}{} x_1^{(k+1)}&=&b_{1,2}x_2^{(k)}&+b_{1,3}x_3^{(k)}&+\cdots&+b_{1,n}x_n^{(k)}&+g_{1},\\ x_2^{(k+1)}&=b_{2,1}x_1^{(k+1)}&&+ b_{2,3}x_3^{(k)}&+\cdots&+b_{2,n}x_n^{(k)}&+g_2,\\ \cdots \\ x_n^{(k+1)}&=b_{n,1}x_{1}^{(k+1)}&+b_{n,2}x_2^{(k+1)}&+b_{n,3}x_3^{(k+1)}&+\cdots&+b_{n,n-1}x_{n-1}^{(k+1)}&+g_n. \end{array}\right. \]

如果将\(I-D^{-1}A\)，具体地，令\(A=D-L-U\)，\(I-D^{-1}A=D^{-1}(L+U)\)，就将原先的迭代矩阵拆成两个严格三角矩阵\(D^{-1}L\)和\(D^{-1}U\)，即

\[D^{-1}L=\begin{pmatrix} 0 \\ b_{2,1} & 0 \\ b_{3,1} & b_{3,2} & 0 \\ \vdots & \vdots & & \ddots \\ b_{n,1} & b_{n,2} & b_{n,3}&\cdots & 0 \end{pmatrix},\quad D^{-1}U=\begin{pmatrix} 0 & b_{1,2} & b_{1,3} & \cdots & b_{1,n} \\ & 0 & b_{2,3} & \cdots & b_{2,n} \\ & & \ddots & & \vdots\ \\ & & & 0 & b_{n-1,n} \\ & & & & 0 \end{pmatrix}. \]

就可以将分量迭代写成向量形式，即

\[x^{(k+1)}=D^{-1}Lx^{(k+1)}+D^{-1}Ux^{(k)}+D^{-1}b,\\ x^{(k+1)}=(D-L)^{-1}Ux^{(k)}+(D-L)^{-1}b. \]

所以\(B_{G}=(D-L)^{-1}U\)，\(f_{G}=(D-L)^{-1}b\)，就得到Gauss-Seidel迭代\(x=B_{G}x+f_{G}\)。注意，矩阵迭代格式只是便于表达，实际运用时，选择分量形式会更合适，也更便于计算。

例：写出用两种迭代法计算下述方程的过程。

\[A=\begin{bmatrix} 1 & 2 & -2 \\ 1 & 1 & 1 \\ 2 & 2 & 1 \end{bmatrix},\quad b=\begin{bmatrix} 1 \\ 1 \\ 1 \end{bmatrix}. \]

解：改写方程为便于迭代的形式：

\[\left\{\begin{array}{} x_1+2x_2-2x_3=1,\\ x_1+x_2+x_3=1,\\ 2x_1+2x_2+x_3=1, \end{array}\right.\Rightarrow \left\{\begin{array}{} x_1=-2x_2+2x_3+1,\\ x_2=-x_1-x_3+1,\\ x_3=-2x_1-2x_2+1. \end{array}\right. \]

从而Jacobi迭代法是

\[\left\{\begin{array}{} x_1^{(k+1)}=-2x_2^{(k)} +2x_3^{(k)} +1,\\ x_2^{(k+1)} =-x_1^{(k)} -x_3^{(k)} +1,\\ x_3^{(k+1)} =-2x_1^{(k)} -2x_2^{(k)} +1. \end{array}\right.\quad B_J=\begin{bmatrix} 0 & -2 & 2 \\ -1 & 0 & -1 \\ -2 & -2 & 0 \end{bmatrix}, \]

Gauss-Seidel迭代法是

\[\left\{\begin{array}{} x_1^{(k+1)} = -2x_2^{(k)} + 2x_3^{(k)} +1,\\ x_2^{(k+1)} = -x_1^{(k+1)} -x_3^{(k)} +1,\\ x_3^{(k+1)} =-2x_1^{(k+1)} -2x_2^{(k+1)} +1. \end{array}\right. \]

松弛法(SOR)可以看作Gauss-Seidel迭代法的加速，由\(x^{(k+1)}=D^{-1}Lx^{(k+1)}+D^{-1}Ux^{(k)}+D^{-1}b\)，有

\[\Delta x = x^{(k+1)}-x^{(k)}=D^{-1}L^{(k+1)}+D^{-1}Ux^{(k)}+D^{-1}b-x^{(k)}, \]

可以看作\(x^{(k+1)}+\Delta x\)。如果在误差前面加上一个松弛因子\(\omega\)，记\(x^{(k+1)}=x^{(k)}+\omega\Delta x\)，就有

\[x^{(k+1)}=x^{(k)}+\omega\Delta x=(1-\omega)x^{(k)}+\omega[D^{-1}Lx^{(k+1)}+D^{-1}Ux^{(k)}+D^{-1}b]. \]

用分量形式将更直观，有

\[x_i^{(k+1)}=(1-\omega)x_i^{(k)}+\omega\left(\sum_{j=0}^{i-1}b_{i,j}x_j^{(k+1)}+\sum_{j=i+1}^{n}b_{i,j}x_j^{(k)}+g_i \right). \]

这是上一个分量与新算出分量的加权平均，且同样具有Gauss-Seidel法的节省空间优点。可以计算其迭代矩阵，有

\[(D-\omega L)x^{(k+1)}=[(1-\omega)D+\omega U]x^{(k)}+\omega b,\\ x^{(k+1)}=(D-\omega L)^{-1}[(1-\omega)D+\omega U]+\omega(D-\omega L)^{-1}b.\\ x=B_{\omega}x+f_{\omega}. \]

迭代法的收敛

迭代法的迭代序列不一定收敛，能否收敛主要由其迭代矩阵决定。

• 迭代收敛定理：对任何初始向量\(x^{(0)}\)和常数项\(f\)，由迭代格式\(x^{(k+1)}=Mx^{(k)}+f\)产生的迭代向量\(\{x_k\}\)收敛且极限与初值无关的充要条件是
  \[\rho(M)<1. \]
  其中\(\rho(M)\)是\(M\)的谱半径，即特征值绝对值的最大值。

为证明必要性，假设\(\{x_k\}\)会最终收敛到\(x\)，令\(\varepsilon_k=x^{(k)}-x\)，有

\[x^{(k+1)}-x=Mx^{(k)}-Mx=M(x^{k}-x),\\ \varepsilon_{k+1}=M\varepsilon_k=\cdots=M^{k+1}\varepsilon_0. \]

欲使向量序列\(\{M^{k}\varepsilon_0\}\)收敛于\(0\)，必有\(M^{k}\to O\)，也就是\(\rho(M)<1\)。

再证充分性，若\(\rho(M)<1\)，则\(I-M\)可逆，从而\((I-M)x=f\)有唯一解\(x\)，上述误差收敛仍成立。

• 收敛速度：若\(\rho(M)<1\)，称\(-\ln \rho(M)\)为迭代法的收敛速度。

对于Jacobi迭代、Gauss-Seidel迭代与SOR迭代法，有一些收敛定理。

1. 对\(Ax=b\)，若系数矩阵\(A\)具有严格对角优势，则Jacobi迭代法和Gauss-Seidel迭代法收敛。
2. 对\(Ax=b\)，若\(A\)对称正定，则Gauss-Seidel迭代法收敛。
3. 为使松弛法收敛，必要条件是\(|1-\omega|<1\)，即\(\omega\in(0,1)\)。
4. 对\(Ax=b\)，若\(A\)具有严格对角优势，则松弛法收敛的充分条件是\(\omega\in (0,1]\)。
5. 对\(Ax=b\)，若\(A\)是对角线元素均为正实数，则\(\omega\in(0,2)\)时松弛法收敛的充要条件是\(A\)正定。

最后，由于松弛法的松弛因子是可选的，因此存在一个最优松弛因子使松弛法收敛最快。在\(A\)正定时，定理保证了\(\omega\in(0,2)\)时最优松弛因子是

\[\omega_{opt}=\frac{2}{1+\sqrt{1-\rho^2(B)}}, \]

这里\(\rho(B)\)是相应Jacobi方法迭代矩阵\(B\)的谱半径。

例：设

\[A=\begin{pmatrix} 4 & 0 & 3 \\ 0 & 3 & 0 \\ 1 & 0 & 2 \end{pmatrix},\quad b=\begin{pmatrix} 1 \\ 1 \\ 1 \end{pmatrix}, \]

松弛迭代格式为\(x^{(k+1)}=x^{(k)}-\omega(Ax^{(k)}-b)\)，求迭代法收敛时\(\omega\)的取值范围与最优松弛因子。

解：迭代矩阵为\(I-\omega A\)，由于\(A\)的特征值分别为\(1,3,5\)，所以\(I-\omega A\)的特征值为

\[1-\omega,\quad 1-3\omega,\quad 1-5\omega. \]

为使迭代法收敛，有

\[-1< 1-\omega<1,\\ -1<1-3\omega<1,\\ -1<1-5\omega<1, \]

解得\(\omega\in(0,\dfrac{2}{5})\)。对最优松弛因子，要使其谱半径最小，应有

\[1-\omega=-(1-5\omega),\\ \omega_{opt}=\frac{1}{3}. \]

非线性方程组

迭代法与收敛速度

对线性方程组，最常使用的求解方法是迭代法，即将\(f(x)=0\)通过某种方法转化为\(x=g(x)\)，然后取一个初值一直迭代，如果迭代序列收敛，就得到一个解。

• 收敛定理：若\(x=g(x)\)中，\(\exists L<1\)使\(g(x)\)满足

  \[|g(x_1)-g(x_2)|\le L|x_1-x_2|,\quad \forall x_1,x_2\in I. \]

  则在区间\(I\)上任取初值\(x_0\)，它将收敛到附近的解\(\xi\)。为验证此条件，往往验证\(|g'(x)|\le L<1\)。

• 误差估计：在满足上述定理要求时，有

  \[|x_n-\xi|\le \frac{L^{n}}{1-L}|x_1-x_0|. \]

收敛速度：为求非线性方程\(f(x)=0\)的根，往往用某个解序列\(\{x_n\}\)来逼近方程，设解序列收敛于\(\alpha\)。对一种特定的求解方法，若存在常数\(p\ge 1\)，有\(c>0\)使

\[\lim_{n\to \infty}\frac{|x_{n+1}-\alpha|}{|x_n-\alpha|^{p}}=c, \]

则称该解序列\(\{x_n\}\)是\(p\)阶收敛的，该求解方法是一个\(p\)阶方法，特别当\(p=1\)时称为线性收敛，\(p\ge 1\)时称为超线性收敛。以下给出一些迭代收敛速度的相关定理，对它们的证明只需不断使用Lagrange中值定理：

• 设\(\varphi(\xi)=\xi\)，且\(|\varphi'(\xi)|<1\)，\(\varphi'(x)\)在\(\xi\)附近连续，则迭代法\(x=\varphi(x)\)局部收敛。
• 设\(\varphi(\xi)=\xi\)，且\(\varphi'(\xi)=\cdots=\varphi^{(p-1)}(\xi)=0\)，\(\varphi^{(p)}(\xi)\ne 0\)，则迭代法\(x=\varphi(x)\)在\(\xi\)附近\(p\)阶收敛。

例：\(x=\varphi(x)\)的迭代法中，取\(\varphi(x)=Ax+Bx^2+Cx^3\)，方程有一个根为\(\xi=\dfrac{1}{7}\)。为使迭代具有三阶收敛速度，求\(A,B,C\)的值。

解：我们有

\[\varphi'(x)=A+2Bx+3Cx^2,\\ \varphi''(x)=2B+6Cx. \]

为使迭代法有三阶收敛速度，应有\(\varphi'(\xi)=\varphi''(\xi)=0\)，结合\(\varphi(\xi)=\xi\)，有

\[\frac{1}{7}A+\frac{1}{49}B+\frac{1}{343}C=\frac{1}{7},\\ A+\frac{2}{7}B+\frac{3}{49}C=0,\\ 2B+\frac{6}{7}C=0. \]

解得

\[(A,B,C)=(3,-21,49). \]

二分法较为特殊，它能够保证每次二分都缩小绝对误差限到原先的一半，因而按照收敛速度的预期，它是线性收敛的；但是二分法不能保证上述极限是收敛的，即

\[\left\{\frac{|x_{n+1}-\alpha|}{|x_{n}-\alpha|}\right\} \]

可能不存在极限，甚至可能是无界数列。不过，我们仍称其为线性收敛的方法。

加速迭代

松弛加速法：对一个特定的收敛迭代\(x=\varphi(x)\)，两边加上\(\lambda x\)将其改写为

\[x=\frac{\lambda}{1+\lambda}x+\frac{1}{1+\lambda}\varphi(x):=\psi(x),\\ \psi'(x)=\frac{1}{1+\lambda}(\lambda+\varphi'(x)) \]

为使迭代方程更快，应取\(\lambda=-\varphi'(\xi)\)，用\(x_k\)替代，只要\(\varphi'(x_k)\ne -1\)，就令

\[\omega_k=\frac{1}{1+\lambda_k}=\frac{1}{1-\varphi'(x_k)},\\ x_{n+1}=(1-\omega_n)x_n+\omega_n\varphi(x_n). \]

Aitken加速法避免了松弛法要计算导数的缺点，步骤如下：

1. 根据\(x_n\)求得\(y_n=\varphi(x_n)\)和\(z_n=\varphi(y_{n})\)。

2. 估计误差，令\(c=\varphi'(\xi)\)，\(\varepsilon=x_n-x\)，有

   \[y_n-\xi=\varphi(x_n)-\varphi(\xi) \approx \varphi'(\xi)\varepsilon=c\varepsilon,\\ z_n-\xi=\varphi(y_n)-\varphi(\xi)\approx \varphi'(\xi)(y_n-\xi)=c^2\varepsilon,\\ y_n-x_n\approx (c-1)\varepsilon,\\ z_n-y_n\approx c(c-1)\varepsilon, \]

   从而

   \[c^2\varepsilon=\frac{[c(c-1)\varepsilon]^2}{(c-1)^2\varepsilon}=\frac{(y_n-z_n)^2}{(z_n-y_n)-(y_n-x_n)}=\frac{(y_n-z_n)^2}{z_n-2y_n+x_n}, \]

3. 由\(\xi\approx z_n-c^2\varepsilon\)，得到新解：

   \[x_{n+1}=z_n-\frac{(y_n-z_n)^2}{z_n-2y_n+x_n}. \]

对初始迭代方程\(x=\varphi(x)\)的选择，常使用Newton迭代法，取\(f(x)\)的一阶Taylor级数改写而成：

\[0=f(x)=f(x_0)+f'(x_0)(x-x_0),\\ x_1=x_0-\frac{f(x_0)}{f'(x_0)}=\varphi(x_0). \]

Newton法在单根附近是二阶收敛的，在重根附近一阶收敛，如果想使得Newton迭代法在\(m\)重根附近也具有二阶收敛速度，需将迭代方程改写为

\[x_{n+1}=x_{n}-m\frac{f(x_n)}{f'(x_n)}. \]

非线性方程组牛顿迭代法

对\(k\)个未知数\(X=(x_1,x_2,\cdots,x_k)'\)以及\(k\)个方程\(F(X)=(u_1(X),\cdots,u_k(X))’\)，记

\[F'(X)=\frac{\partial (u_1,\cdots,u_k)}{\partial(x_1,\cdots,x_k)}, \]

类似推广得到

\[X^{(k+1)}=X^{(k)}-[F'(X^{(k)})]^{-1}F(X^{(k)}), \]

这就是非线性方程组的牛顿迭代法。

常微分方程初值问题

一步解法

考虑常微分方程初值问题：

\[\frac{\mathrm{d}y}{\mathrm{d}x}=f(x,y),\quad x\in[a,b];\\ y(x_0)=y_0. \]

Euler法即节选向量场中的部分点，用一条折线来替代函数在该点附近的性状，折线的终点即下一步计算的起点，即

\[y(x_0)=y_0,\quad y_{i+1}=y_i+hf(x_i,y_i). \]

注意到Euler法中的核心关系式为\(\displaystyle{\frac{y_{i+1}-y_i}{h}\approx f(x_i,y_i)}\)，即用向前差商代替导数真值，如用向后差商来替代，即\(\displaystyle{\frac{y_{i+1}-y_i}{h}\approx f(x_{i+1},y_{i+1})}\)，则得到向后Euler方法：

\[y(x_0)=y_0,\quad y_{i+1}=y_{i}+f(x_{i+1},y_{i+1}). \]

但要注意此方程无法直接代入得到\(y_{i+1}\)，因为求\(y_{i+1}\)的过程中用到了\(y_{i+1}\)，因此应将其视作一个关于\(y_{i+1}\)方程，使用迭代法求解，初值自然地使用向前差分的结果，即

\[y_{i+1}^{(0)}=y_i+hf(x_i,y_i),\\ y_{i+1}^{(k)}=y_i+hf(x_{i+1},y_{i+1}^{(k-1)}). \]

迭代法收敛的条件是\(\{y_{i-1}^{(k)}\}\)是一个Cauchy序列，不妨假设\(f(x,y)\)满足Lipschitz条件（这是ode初值问题收敛的充分条件），从而

\[|y_{i+1}^{(k+1)}-y_{i+1}^{(k)}|=h|f(x_{i+1},y_{i+1}^{(k)})-f(x_{i+1},y_{i+1}^{(k-1)})|\le hL|y_{i+1}^{(k)}-y_{i+1}^{(k-1)}|, \]

这说明只要\(hL< 1\)，上述迭代即收敛。实际应用时，如果\(h\)足够小，往往仅迭代一次。

一个自然的想法是，将向前Euler方法与向后Euler方法相结合，用它们的平均得到一个新的求解方法，这便是改进Euler方法，具体表现为\(y_{i+1}=y_i+\displaystyle{\frac{h}{2}[f(x_i,y_i)+f(x_{i+1},y_{i+1})]}\)。只迭代一次，就是

\[\tilde y_{i+1}=y_i+hf(x_i,y_i),\\ y_{i+1}=y_i+\frac{h}{2}[f(x_i,y_i)+f(x_{i+1},\tilde y_{i+1})]. \]

称改进Euler方法的公式为梯形公式，是基于下式：

\[y(x_{i+1})=y(x_i)+\int_{x_i}^{x_{i+1}}f(t,y(t))\mathrm{d}t\stackrel{*}\approx y(x_i)+\frac{h[f(x_i,y_i)+f(x_{i+1},y_{i+1})]}{2}, \]

其中\((*)\)部分是使用梯形公式对积分进行估计的结果。

例：用改进Euler方法计算\(y(0)=1\)，\(y'=x+y+1\)，取步长\(h=1\)，计算\(3\)个点。

解：令\(f(x,y)=x+y+1\)，有

\[\tilde y_1=1+f(0,1)=3,\\ y_1=1+\frac{[f(0,1)+f(1,3)]}{2}=4.5;\\ \tilde y_2=4.5+f(1,4.5)=11,\\ y_2=4.5+\frac{[f(1,4.5)+f(2,11)]}{2}=17.25;\\ \tilde y_3=14.75+f(2,14.75)=32.5,\\ y_3=14.75+\frac{[f(2,14.75)+f(3,32.5)]}{2}=41.875. \]

误差估计与收敛性

对常微分方程初值问题的误差估计，一般使用局部截断误差概念，这是因为如果从方法本身分析\(y_i-y_i^*\)的大小往往比较复杂。

• 局部截断误差：若\(y_{i+1}\)是在\(y_j=y(x_j),j\le i\)的假定下，由某数值方法得到的\(y(x_{i+1})\)的近似值，则称该数值方法的局部截断误差为\(R_{i+1}=y(x_{i+1})-y_{i+1}\)。
• 阶数：若某数值方法的局部截断误差为\(O(h^{p+1})\)，则称此方法的阶数为\(p\)，或称其为\(p\)阶方法。

分析一个方法的好坏，Taylor公式与全微分是常用的手段，例如对Euler方法，有

\[y(x_{i+1})=y(x_i+h)=y(x_i)+hy'(x_i)+\frac{h^2y''(x_i)}{2}+O(h^3),\quad \xi\in[x_i,x_{i}+h], \]

注意到\(y'(x_i)=f(x_i,y_i)\)精确成立，所以

\[R_{i+1}=y(x_{i+1})-[y(x_i)+hf(x_i,y_i)]=\frac{h^2y''(x_i)}{2}+O(h^3)=O(h^2), \]

当\(y''(x)\)在区间上有界时\(R_{i+1}=O(h^2)\)，故Euler方法是一个\(1\)阶方法。而对一次迭代的向后Euler方法，即

\[y(x_{i+1})=y(x_i)+hf(x_{i}+h,y_i+hf(x_i,y_i)), \]

因\(\displaystyle{\mathrm{d}f(x,y)=\frac{\partial f}{\partial x}\mathrm{d}x+\frac{\partial f}{\partial y}\mathrm{d}y}\)，特别在\((x_i,y_i)\)点处，

\[\mathrm{d}f(x_i,y_i)=\frac{\partial f(x_i,y_i)}{\partial x}\mathrm{d}x+\frac{\partial f(x_i,y_i)}{\partial y}\mathrm{d}y=\left(\frac{\partial f(x_i,y_i)}{\partial x}+\frac{\partial f(x_i,y_i)}{\partial y}y'(x_i) \right)\mathrm{d}x, \]

此式给出了一个由\(x\)的微小变化，先导致\(y\)变化，再导致\(f(x,y)\)微小变化的关系式，等同于\(x\)的一元函数，在这里就是\(y'(x)\)，所以\(\dfrac{\mathrm{d}f}{\mathrm{d}x}=y''(x)\)，于是

\[\begin{aligned} f(x_i+h,y_i+hf(x_i,y_i))&=f(x_i,y_i)+h\frac{\partial f(x_i,y_i)}{\partial x}+hf(x_i,y_i)\frac{\partial f(x_i,y_i)}{\partial y}+O(h^2)\\ &=y'(x_i)+h\frac{\mathrm{d}f(x_i,y_i)}{\mathrm{d}x}+O(h^2)\\ &=y'(x_i)+hy''(x_i)+O(h^2). \end{aligned} \]

这说明对向后Euler方法，有

\[R_{i+1}=y(x_{i+1})-y_{i+1}=-\frac{h^2y''(x_i)}{2}+O(h^3)=O(h^2), \]

即向后Euler方法也是\(1\)阶方法。但向前Euler方法和向后Euler方法的误差项里，\(h^2\)的系数正好相反，因此将它们作平均，得到的改进Euler方法就满足\(R_{i+1}=O(h^3)\)，说明改进Euler方法是一个二阶方法。

例：证明用下式求解初值问题

\[y_{i+1}=y(x_i)+hf\left(x_i+\frac{1}{2}h,y_i+\frac{1}{2}hf(x_i,y_i) \right), \]

该方法是一个\(2\)阶方法，即局部截断误差为\(O(h^3)\)。

解：有

\[\begin{aligned} &\quad f\left(x_i+\frac{1}{2}h,y_i+\frac{1}{2}hf(x_i,y_i) \right)\\ &=f(x_i,y_i)+\frac{h}{2}\left(\frac{\partial f(x_i,y_i)}{\partial x}+\frac{\partial f(x_i,y_i)}{\partial y}f(x_i,y_i) \right)+O(h^2)\\ &=y'(x_i)+\frac{h}{2}y''(x_i)+O(h^2), \end{aligned} \]

且

\[y(x_{i+1})=y(x_i)+hy'(x_i)+\frac{h^2}{2}y''(x_i)+O(h^3), \]

代入即得。

考虑方法的收敛性，指的是\(h\to 0\)时，数值解\(y_i\)是否收敛到方程准确解\(y(x_i)\)。

• 收敛：如果一个数值方法对任意的点\(x_i=x_0+ih\)，当\(h\to 0\)时都有\(y_i\to y(x_i)\)，则称该方法是收敛的。

现考虑Euler方法的收敛性，假定\(f(x,y)\)在区间上服从Lipschitz条件。在上面的分析中易得，如果\(y''(x)\)在区间上有界，则局部截断误差\(T_{i+1}=y(x_{i+1})-y_{i+1}\)满足\(|T_{i+1}|\le Ch^2\)，这里\(C=\displaystyle{\frac{1}{2}\max y''(x)}\)。记总体截断误差为\(e_{i}=y(x_i)-y_i\)，则

\[\begin{aligned} |e_{i+1}|&=|y(x_{i+1})-y_{i+1} |\\ &=|y(x_i)+hf(x_i,y(x_i))+T_{i+1}-[y_i+hf(x_i,y_i)]|\\ &\le|y(x_i)-y_i|+h|f(x_i,y(x_i))-f(x_i,y_i)|+|T_{i+1}|\\ &\le (1+hL)|e_i|+|T_{i+1}|,\quad \forall i,\\ |e_n|&\le |T_n|+(1+hL)|e_n|\\ &\le |T_n|+(1+hL)|T_{n-1}|+(1+hL)^2|e_{n-1}|\\ &\le |T_n|+(1+hL)|T_{n-1}|+(1+hL)^2|T_{n-2}|+(1+hL)^3|e_{n-2}|\\ &\le \cdots \\ &\le |T_n|+(1+hL)|T_{n-1}|+(1+hL)^2|T_{n-2}|+\cdots+(1+hL)^{n-1}|T_1|\\ &\le Ch^2\left[\sum_{k=0}^{n-1}(1+hL)^{k} \right]\\ &=\frac{(1+hL)^{n}-1}{1+hL-1}Ch^2\\ &\le\frac{Ch}{L}(e^{nhL}-1). \end{aligned} \]

这说明Euler法在各点上逐点收敛。

稳定性

稳定性指的是数值方法的稳定性，即误差的积累能否得到控制，以下绝对稳定的概念是用于刻画稳定性的一个标准。

• 绝对稳定：用一个数值方法求解微分方程\(y'=\lambda y\)，其中\(\lambda\)是一个复常数，对给定步长\(h>0\)，在计算\(y_i\)时引入了误差\(\rho_i\)。若这个误差在计算后面的\(y_{i+k},k=1,2,\cdots\)中引起的误差绝对值均不增加，就说这个数值方法对于步长\(h\)和复数\(\lambda\)是绝对稳定的。
• 绝对稳定域：一个绝对稳定的方法中，\(\lambda h\)的取值范围。

一般，为计算绝对稳定域，只需计算相邻两步的误差之比的绝对值（模），如果该比值\(\left|\dfrac{\rho_{i+1}}{\rho_i}\right|\le 1\)，就是绝对稳定的。如对Euler方法，有

\[y_{i+1}=y_{i}+h\lambda y_{i}, \]

若引入了误差\(\rho_{i}\)，就有

\[y_{i+1}+\rho_{i+1}=y_i+\rho_i+h\lambda(y_i+\rho_i),\\ \rho_{i+1}=(1+h\lambda)\rho_i,\\ \left|\frac{\rho_{i+1}}{\rho_i}\right|=|1+h\lambda|\le 1, \]

这说明绝对稳定域是以\((-1,0)\)为圆心，\(1\)为半径的一个圆。

例：求向后Euler方法

\[y_{i+1}=y_i+hf(x_{i+1},y_{i+1}) \]

的绝对稳定域。

解：有

\[\rho_{i+1}=\rho_i+h\lambda \rho_{i+1},\\ \left|\frac{\rho_{i+1}}{\rho_i} \right|=\left|\frac{1}{1-h\lambda }\right|\le 1, \]

得\(|1-h\lambda |\ge 1\)，这是以\((1,0)\)为半径，\(1\)为圆心的圆。

对一般的方差，要选择合适的步长\(h\)使得该方法稳定，可以取

\[\lambda\approx \frac{\partial f}{\partial y}. \]

R-K法

本部分仅针对计算方法。

由\(y(x_{i+1})-y(x_i)=y'(\xi)h\)可知，求解微分方程的核心，是对这一区间上的平均变化率\(y'(\xi)\)进行预估，如果预估方法好，则方法自然精度高。为估计出这个平均变化率，可通过求解\([x_i,x_{i+1}]\)内几个点的斜率近似值，加权平均得到平均变化率，具体地，有

\[y_{i+1}=y_i+h\left(\sum_{i=1}^{m}\alpha_iK_i \right),\\ K_1=f(x_i,y_i),\\ K_2=f(x_{i}+\lambda_2h,y_i+\mu_2h),\\ \cdots\\ K_m=f(x_i+\lambda_mh,y_i+\mu_mh). \]

这里的待估参数在诸\(\lambda_i,\mu_i\)，为使对平均变化率的预估较好，应尽可能使得求积公式的局部截断误差具有较高的阶数。具体做法仍是Taylor展开和全微分，这里略去过程，给出以下的经典R-K公式：

\[\left\{\begin{array}{} \displaystyle{y_{i+1}=y_i+\frac{h}{6}(K_1+2K_2+2K_3+K_4),}\\ K_1=f(x_i,y_i),\\ \displaystyle{K_2=f\left(x_i+\frac{h}{2},y_i+\frac{h}{2}K_1 \right),}\\ \displaystyle{K_3=f\left(x_i+\frac{h}{2},y_i+\frac{h}{2}K_2 \right),}\\ K_4=f(x_i+h,y_i+hK_3). \end{array}\right. \]

我们所指的R-K法，经常都是这个方法，它的局部截断误差为\(O(h^5)\)，是一个\(4\)阶方法。

求矩阵的特征值

幂法

幂法的想法是，由于特征向量线性无关，\(\mathbb{R}^{n}\)中的任一向量能通过特征向量线性表示，按特征值绝对值大小\(|\lambda_i|\)排列，它们对应的特征向量是\(v_i\)。现假定绝对值最大的特征值（称为主特征值）唯一，且仅有唯一的特征向量\(v_1\)，如果\(\displaystyle{x=\sum_{i=1}^{n}\alpha_iv_i}\)，则

\[A^{k}x=\sum_{i=1}^{n}\alpha_i\lambda_i^{k}v_i=\lambda_1^{k}\sum_{i=1}^{n}\alpha_i\left(\frac{\lambda_i}{\lambda_1} \right)^{k}v_i\to\lambda_1^{k}\alpha_1v_1. \]

即\(x^{(k)}=\lambda_1^{k}(\alpha_1v_1+\varepsilon_k)\)，经过充分的迭代后，此序列的极限为矩阵\(A\)的特征向量。实际计算时，应采用如下方法：

\[\lim_{k\to \infty}\frac{x^{(k+1)}_i}{x^{(k)}_i}=\lambda_1\lim_{k\to \infty}\frac{\alpha v_{1,i}+\varepsilon_{k+1,i}}{\alpha_1v_{1,i}+\varepsilon_{k,i}}=\lambda_1, \]

也就是对向量序列的某一分量连续作比，其极限就是矩阵的主特征值，且收敛速度为\(\gamma=\dfrac{\lambda_2}{\lambda_1}\)。在实际应用时，需按照以下的改进幂法进行，它比普通的幂法多增加了一个归一化的步骤，使得迭代序列的最大分量为\(1\)：

1. 首先取非零的初始向量\(x^{(0)}=y^{(0)}\ne 0\)。
2. 对任意\(k=1,2,\cdots\)，进行幂法迭代：\(x^{(k)}=Ay^{(k-1)}\)。
3. 归一化：\(y^{(k)}=\dfrac{x^{(k)}}{\max\{x^{(k)}_i\}}\)。

通过此改进幂法，得到以下结论：

\[\lim_{k\to \infty}y^{(k)}=\frac{v_1}{\max\{v_{1,i}\}},\\ \lambda_1=\lim_{k\to \infty}\max \{Ay^{(k)}_i\}. \]

为使此方法加速，可以采用原点平移法，即选取合适的\(p\)，对\(B=A-pI\)执行幂法，因\(B\)与\(A\)有相同的特征向量，且特征值之间差值为\(p\)。使用此幂法的前提是\(\lambda_1-p\)仍为矩阵\(B\)的主特征值，且收敛速度得到加快。

下讨论主特征值不唯一的情形。如果\(|\lambda_1|=|\lambda_2|=\cdots=|\lambda_r|\)，情况将有所不同。如\(\lambda_1=\cdots=\lambda_r\)，此时

\[x^{(k)}=\lambda_1^{k}\left(\sum_{i=1}^{r}\alpha_iv_i+\sum_{i=r+1}^{n}\alpha_i\left(\frac{\lambda_i}{\lambda_1} \right)^{k}v_n \right). \]

这说明幂法仍收敛，但向量极限\(\displaystyle{\sum_{i=1}^{r}\alpha_iv_i}\)将有赖于初值\(x^{(0)}\)的选择。又如\(\lambda_1=-\lambda_2\)，此时

\[x^{(k)}=\lambda_1^{k}\left(\alpha_1v_1+(-1)^{k}\alpha_2v_2+\sum_{i=3}^{n}\alpha_i\left(\frac{\lambda_r}{\lambda_1} \right)^{k}v_i \right). \]

可以看出，虽然幂法不收敛，但是其子序列\(\{x^{(2k)}\}\)和\(\{x^{(2k+1)}\}\)分别收敛，故仍可生效。但如果主特征值是复数，则情况将变得复杂，尽管仍可以应用，但是其子序列将不收敛。

反幂法

为求\(A\)按模最小的特征值，如果\(A\)没有零特征值，即可逆，那么\(A^{-1}\)的最大特征值的倒数就是\(A\)的特征值。不过，对\(A\)的求逆往往不直接求，因当\(n\)充分大时\(A^{-1}\)的计算比较复杂。

注意到构造迭代序列\(x^{(k+1)}=A^{-1}x^{(k)}\)，等价于求解线性方程组\(Ax^{(k+1)}=x^{(k)}\)，因此可以通过求解线性方程组来获得迭代序列，同时，每一步都进行归一化。具体操作如下：

1. 对\(A\)进行\(LU\)分解，取\(x^{(0)}=y^{(0)}\ne0\)。

2. 对\(K=1,2,\cdots\)，进行反幂法求解：

   \[Lz=y^{(k-1)},\\ Ux^{(k)}=z. \]

3. 归一化：令\(y^{(k)}=\dfrac{x^{(k)}}{\max\{x^{(k)}_i\}}\)。

得到的结论与幂法的相似，最后将特征值求倒数即可。

Jacobi方法

Jacobi方法适用于实对称矩阵特征值求算，且可以求出所有的特征值，它基于以下两个原则：实对称矩阵必可正交对角化，相似矩阵的特征值相等。因此，通过正交矩阵将实对称矩阵对角化，这样对角线上所有的元素都是矩阵的特征值。实际使用Jacobi方法执行正交旋转时，一般是采用二维旋转阵：

\[P(i,j,\theta)=\begin{bmatrix} 1 \\ & \ddots \\ && 1 \\ &&& \cos\theta&&&&\sin \theta \\ &&&& 1 \\ &&&&& \ddots \\ &&&&&& 1 \\ &&&-\sin \theta&&&& \cos\theta \\ &&&&&&&& 1 \\ &&&&&&&&& \ddots \\ &&&&&&&&&& \ddots \end{bmatrix} \]

这里\(\cos\theta\)和\(\sin \theta\)分列第\(i,j\)行与第\(i,j\)列的交点。确定了算法以后，需要明确两个问题：

• 每一步如何选择\(i,j\)。
• 如何选择旋转角\(\theta\)。

关于第一个问题，一般选择\(a_{i,j}\)为非对角线元素中的绝对值最大者，而关于第二个问题，可以进行以下的实验：设二阶矩阵\(A=(a_{i,j})\)，这里\(a_{1,2}=a_{2,1}\)，于是

\[PAP^{-1}=\begin{bmatrix} a_{1,1}\cos^2\theta+a_{2,2}\sin^2\theta+a_{1,2}\sin2\theta & \dfrac{a_{2,2}-a_{1,1}}{2}\sin2\theta+a_{1,2}\cos2\theta \\ *& a_{1,1}\sin ^2\theta+a_{2,2}\cos^2\theta-a_{1,2}\sin2\theta \end{bmatrix}. \]

为使它为对角阵，应有\(\cot2\theta=\dfrac{a_{1,1}-a_{2,2}}{2a_{1,2}}\)，特别当\(a_{1,1}=a_{2,2}\)时，有\(\theta=\dfrac{\pi}{4}\)。