科学计算/计算方法备考纲要
科学计算/计算方法 备考纲要
声明:本备考纲要以科学计算为主,如有计算方法的相关内容将额外标注。
误差
关于误差的分类,需掌握绝对误差与相对误差。
-
绝对误差:若用\(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\)的近似时,有
事实上,可将相应的误差看作相应的微分,即\(\mathrm{d}x\approx x^*-x\),\(\mathrm{d}y\approx y^*-y\),于是
此时相对误差为\(\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多项式,需要用到差商。一般地,有
且重结点差商与导数之间具有如下的关系:下设\(k\)为结点重数,则
例题:设\(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)\),则多项式插值的误差为
证明即对任何\(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)\),则插值多项式的误差为
类似构造\(\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\)的点。
始终记住,求导次数等于条件数。
三次样条插值
本部分只针对计算方法,因科学计算对样条插值的计算不作要求。
三次样条插值需要增加两个边界条件,常见的有:
- 端点处一阶导数值:\(S'(x_0)\)和\(S'(x_n)\)已知。
- 端点处二阶导数值:\(S''(x_0)\)和\(S''(x_n)\)已知。
- 周期条件:\(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)\)是不高于三次的多项式,故其二阶导数不超过一次,可构造插值基函数为
称\(h_i=x_i-x_{i-1}\)为步长。利用基函数思想对此函数二次积分(自行体会系数的配凑),得到
这样,每一个方程都被一个\(M_i\)唯一确定,为求\(M_i\),需用一阶导数条件,即\(S_i'(x_{i})=S_{i+1}'(x_i)\),所以
利用它们相等,\(i=1,2,\cdots,n-1\),可得
为求解,需将边界条件予以应用。最简单的边界条件是给出\(S''(x_0)\)和\(S''(x_n)\)时的边界,此时\(M_0\)和\(M_n\)直接给出;而如果给出\(S'(x_0)\)和\(S'(x_n)\),则要对\(S_1(x)\)和\(S_n(x)\)求导,得到
数据拟合
拟合即使用最小二乘法计算,由于正交多项式不作要求,因此只考察基础最小二乘法,至多扩展至多项式拟合。回归方程为\(y=(a_0,a_1,\cdots,a_k)(1,x,\cdots,x^k)'\),考察残差平方和为
对\(k+1\)个参数分别求偏导,就得到法方程组为
由此可以解出最优统计量。特别当回归方程为直线时,有
这是回归分析的基础结果。
例:给定以下超定方程组
\[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. \]
数值微分与数值积分
数值微分
我们所用的数值微分基于插值多项式构造,即用插值多项式的微商来替代函数的微商,且一般用等距节点构造的插值多项式求微分,常用的有两点式和三点式。对于两点式,有
事实上就是两点之间的斜率。对于三点式,有
对其误差估计,可以直接由插值多项式的误差估计推得,有
注意\(f^{(n+1)}(\xi)\)的导数一般无法求得,但若在插值结点处,则由\(\omega_n(x_i)=0\)可以忽略这一项,此时
对两点式,有
对三点式,则可以类似推导,得到其截断误差为\(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]\)上,不同的求积公式表现为
这里\(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\)分别为
一般只会用到至多四个结点的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\),所以
但当\(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]\)上不变号,所以可以用积分中值定理,得到
而对抛物线公式则有一个小技巧,因为它具有\(4\)阶代数精度,因此可以构造一个\(4\)次多项式来作为插值多项式,不妨选择\(x=\dfrac{1}{2}\)处的重结点,这样做的好处是\(x(x-\dfrac{1}{2})^2(x-1)\)会在\([0,1]\)上不变号,从而可以使用积分中值定理,有
需要注意,如果是在\([a,b]\)上作积分,误差估计需要乘上\((b-a)^{n+2}\),这里\(n\)是方法的代数精度。通过积分中值定理的做法,只能处理梯形公式与抛物线公式的误差估计,对于更多结点数的Cotes系数,则需要更高级的方法。
复合求积与逐次分半
复合求积指的是在小区间上作积分,常常在小区间上使用梯形求积公式或者抛物线求积公式。这里给出复合梯形、抛物线求积公式的误差估计(\(n\)均代表复合求积的复合区间数):
为了保证计算速度的同时保证计算精度,往往不用事前估计来确定需要使用的结点数,而是使用逐次分半法。我们知道\(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_{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=S_n\),然而利用同样的外推手段,可以获得
这就是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\)个方程精确成立:
然而,由于此时的\(\{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\)只有第\(i\)列有元素,主对角线上元素为\(1\),是一个多重倍加阵,由倍加的性质显然
即反向倍加即可。并且我们有
显然\(\det(A^{(1)})=\det(A^{(n)})\)。
为使计算精度增大,常使用列主元消去法,它基于除数不能太小的想法,每次将当前主元位置中,系数最大的那个方程移到当前行。如果当前主元是\(x_i\),\(A^{(i)}\)中第\(i\)列系数最大的是\(a_{k,i}\),则需将第\(k\)行与第\(i\)行互换。在矩阵的操作中,相当于在\(M_i\)处理之前,先用对换矩阵\(P_i\)处理一遍,也就是
三角分解
三角分解基于如下的想法:如果线性方程组的系数矩阵是上三角或者下三角的,则可以分别通过前代与后代,递推地给出方程的解。将一个矩阵分成两个三角矩阵的乘积,就能通过简单的代入求解方程,而\(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\)分解,设
那么计算过程中,我们有
随着\(k\)增大,逐层先\(u_{k,i}\)后\(l_{i,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分解,有
分别计算,就有
类推计算,可得
为避免开方运算,可以使用\(LDL'\)分解,即额外用一个对角矩阵来存储一个需要开方的结果。
误差分析
先给出向量范数的定义:
- 向量范数:满足以下三个性质的映射称为向量范数。
- 正定性:\(\forall x\in\mathbb{R}^{n}\),\(\|x\|\ge 0\),等号成立当期仅当\(x=0\)。
- 齐次性:\(\forall x\in\mathbb{R}^{n}\),\(\lambda\in\mathbb{R}\),成立\(\|\lambda x\|=\lambda\|x\|\)。
- 三角不等式:\(\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\)的最大、最小特征值。可以验证有
- \(\forall A\),\(\mathrm{Cond}(A)\ge 1\)。
- 单位阵,置换阵,正交阵的谱条件数\(\mathrm{Cond}_2=1\)。
- 对任意非零实数\(\alpha\),有\(\mathrm{Cond}(A)=\mathrm{Cond}(\alpha A)\)。
对于某些线性方程组\(Ax=b\),对\(A\)和\(b\)进行微小扰动将极大影响解的精确度,事实上这主要由\(A\)影响。若\(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迭代法相对最为直观,对方程组
先对角元素归一化,即第\(i\)行除以\(a_{i,i}\),得到
此时的系数矩阵出现了一个单位阵\(I\),由于\(Ix=x\),可以构造如此的迭代序列:
从上面的操作中,我们可以提取\(A\)的对角部分为\(D\),这样,第一步我们用\(D^{-1}\)左乘\(Ax+b\)的两边,第二步从\(D^{-1}A\)中分离了一个\(I\),并将其移到等式的另一边,因此Jacobi迭代的步骤是
同时,可以写出其分量形式为
欲理解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)}\),这时候,分量形式就是
如果将\(I-D^{-1}A\),具体地,令\(A=D-L-U\),\(I-D^{-1}A=D^{-1}(L+U)\),就将原先的迭代矩阵拆成两个严格三角矩阵\(D^{-1}L\)和\(D^{-1}U\),即
就可以将分量迭代写成向量形式,即
所以\(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\),有
可以看作\(x^{(k+1)}+\Delta x\)。如果在误差前面加上一个松弛因子\(\omega\),记\(x^{(k+1)}=x^{(k)}+\omega\Delta x\),就有
用分量形式将更直观,有
这是上一个分量与新算出分量的加权平均,且同样具有Gauss-Seidel法的节省空间优点。可以计算其迭代矩阵,有
迭代法的收敛
迭代法的迭代序列不一定收敛,能否收敛主要由其迭代矩阵决定。
- 迭代收敛定理:对任何初始向量\(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\),有
欲使向量序列\(\{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迭代法,有一些收敛定理。
- 对\(Ax=b\),若系数矩阵\(A\)具有严格对角优势,则Jacobi迭代法和Gauss-Seidel迭代法收敛。
- 对\(Ax=b\),若\(A\)对称正定,则Gauss-Seidel迭代法收敛。
- 为使松弛法收敛,必要条件是\(|1-\omega|<1\),即\(\omega\in(0,1)\)。
- 对\(Ax=b\),若\(A\)具有严格对角优势,则松弛法收敛的充分条件是\(\omega\in (0,1]\)。
- 对\(Ax=b\),若\(A\)是对角线元素均为正实数,则\(\omega\in(0,2)\)时松弛法收敛的充要条件是\(A\)正定。
最后,由于松弛法的松弛因子是可选的,因此存在一个最优松弛因子使松弛法收敛最快。在\(A\)正定时,定理保证了\(\omega\in(0,2)\)时最优松弛因子是
这里\(\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\)使
则称该解序列\(\{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). \]
二分法较为特殊,它能够保证每次二分都缩小绝对误差限到原先的一半,因而按照收敛速度的预期,它是线性收敛的;但是二分法不能保证上述极限是收敛的,即
可能不存在极限,甚至可能是无界数列。不过,我们仍称其为线性收敛的方法。
加速迭代
松弛加速法:对一个特定的收敛迭代\(x=\varphi(x)\),两边加上\(\lambda x\)将其改写为
为使迭代方程更快,应取\(\lambda=-\varphi'(\xi)\),用\(x_k\)替代,只要\(\varphi'(x_k)\ne -1\),就令
Aitken加速法避免了松弛法要计算导数的缺点,步骤如下:
-
根据\(x_n\)求得\(y_n=\varphi(x_n)\)和\(z_n=\varphi(y_{n})\)。
-
估计误差,令\(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}, \] -
由\(\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级数改写而成:
Newton法在单根附近是二阶收敛的,在重根附近一阶收敛,如果想使得Newton迭代法在\(m\)重根附近也具有二阶收敛速度,需将迭代方程改写为
非线性方程组牛顿迭代法
对\(k\)个未知数\(X=(x_1,x_2,\cdots,x_k)'\)以及\(k\)个方程\(F(X)=(u_1(X),\cdots,u_k(X))’\),记
类似推广得到
这就是非线性方程组的牛顿迭代法。
常微分方程初值问题
一步解法
考虑常微分方程初值问题:
Euler法即节选向量场中的部分点,用一条折线来替代函数在该点附近的性状,折线的终点即下一步计算的起点,即
注意到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_{i+1}\),因为求\(y_{i+1}\)的过程中用到了\(y_{i+1}\),因此应将其视作一个关于\(y_{i+1}\)方程,使用迭代法求解,初值自然地使用向前差分的结果,即
迭代法收敛的条件是\(\{y_{i-1}^{(k)}\}\)是一个Cauchy序列,不妨假设\(f(x,y)\)满足Lipschitz条件(这是ode初值问题收敛的充分条件),从而
这说明只要\(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})]}\)。只迭代一次,就是
称改进Euler方法的公式为梯形公式,是基于下式:
其中\((*)\)部分是使用梯形公式对积分进行估计的结果。
例:用改进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)=f(x_i,y_i)\)精确成立,所以
当\(y''(x)\)在区间上有界时\(R_{i+1}=O(h^2)\),故Euler方法是一个\(1\)阶方法。而对一次迭代的向后Euler方法,即
因\(\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)\)点处,
此式给出了一个由\(x\)的微小变化,先导致\(y\)变化,再导致\(f(x,y)\)微小变化的关系式,等同于\(x\)的一元函数,在这里就是\(y'(x)\),所以\(\dfrac{\mathrm{d}f}{\mathrm{d}x}=y''(x)\),于是
这说明对向后Euler方法,有
即向后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\),则
这说明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方法,有
若引入了误差\(\rho_{i}\),就有
这说明绝对稳定域是以\((-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\)使得该方法稳定,可以取
R-K法
本部分仅针对计算方法。
由\(y(x_{i+1})-y(x_i)=y'(\xi)h\)可知,求解微分方程的核心,是对这一区间上的平均变化率\(y'(\xi)\)进行预估,如果预估方法好,则方法自然精度高。为估计出这个平均变化率,可通过求解\([x_i,x_{i+1}]\)内几个点的斜率近似值,加权平均得到平均变化率,具体地,有
这里的待估参数在诸\(\lambda_i,\mu_i\),为使对平均变化率的预估较好,应尽可能使得求积公式的局部截断误差具有较高的阶数。具体做法仍是Taylor展开和全微分,这里略去过程,给出以下的经典R-K公式:
我们所指的R-K法,经常都是这个方法,它的局部截断误差为\(O(h^5)\),是一个\(4\)阶方法。
求矩阵的特征值
幂法
幂法的想法是,由于特征向量线性无关,\(\mathbb{R}^{n}\)中的任一向量能通过特征向量线性表示,按特征值绝对值大小\(|\lambda_i|\)排列,它们对应的特征向量是\(v_i\)。现假定绝对值最大的特征值(称为主特征值)唯一,且仅有唯一的特征向量\(v_1\),如果\(\displaystyle{x=\sum_{i=1}^{n}\alpha_iv_i}\),则
即\(x^{(k)}=\lambda_1^{k}(\alpha_1v_1+\varepsilon_k)\),经过充分的迭代后,此序列的极限为矩阵\(A\)的特征向量。实际计算时,应采用如下方法:
也就是对向量序列的某一分量连续作比,其极限就是矩阵的主特征值,且收敛速度为\(\gamma=\dfrac{\lambda_2}{\lambda_1}\)。在实际应用时,需按照以下的改进幂法进行,它比普通的幂法多增加了一个归一化的步骤,使得迭代序列的最大分量为\(1\):
- 首先取非零的初始向量\(x^{(0)}=y^{(0)}\ne 0\)。
- 对任意\(k=1,2,\cdots\),进行幂法迭代:\(x^{(k)}=Ay^{(k-1)}\)。
- 归一化:\(y^{(k)}=\dfrac{x^{(k)}}{\max\{x^{(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\),此时
这说明幂法仍收敛,但向量极限\(\displaystyle{\sum_{i=1}^{r}\alpha_iv_i}\)将有赖于初值\(x^{(0)}\)的选择。又如\(\lambda_1=-\lambda_2\),此时
可以看出,虽然幂法不收敛,但是其子序列\(\{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)}\),因此可以通过求解线性方程组来获得迭代序列,同时,每一步都进行归一化。具体操作如下:
-
对\(A\)进行\(LU\)分解,取\(x^{(0)}=y^{(0)}\ne0\)。
-
对\(K=1,2,\cdots\),进行反幂法求解:
\[Lz=y^{(k-1)},\\ Ux^{(k)}=z. \] -
归一化:令\(y^{(k)}=\dfrac{x^{(k)}}{\max\{x^{(k)}_i\}}\)。
得到的结论与幂法的相似,最后将特征值求倒数即可。
Jacobi方法
Jacobi方法适用于实对称矩阵特征值求算,且可以求出所有的特征值,它基于以下两个原则:实对称矩阵必可正交对角化,相似矩阵的特征值相等。因此,通过正交矩阵将实对称矩阵对角化,这样对角线上所有的元素都是矩阵的特征值。实际使用Jacobi方法执行正交旋转时,一般是采用二维旋转阵:
这里\(\cos\theta\)和\(\sin \theta\)分列第\(i,j\)行与第\(i,j\)列的交点。确定了算法以后,需要明确两个问题:
- 每一步如何选择\(i,j\)。
- 如何选择旋转角\(\theta\)。
关于第一个问题,一般选择\(a_{i,j}\)为非对角线元素中的绝对值最大者,而关于第二个问题,可以进行以下的实验:设二阶矩阵\(A=(a_{i,j})\),这里\(a_{1,2}=a_{2,1}\),于是
为使它为对角阵,应有\(\cot2\theta=\dfrac{a_{1,1}-a_{2,2}}{2a_{1,2}}\),特别当\(a_{1,1}=a_{2,2}\)时,有\(\theta=\dfrac{\pi}{4}\)。