研一课程笔记-数值分析

绪论与预备知识

概论

衡量数值算法的几个标准

衡量算法的优劣有两个标准,一是要有可靠的理论基础,包括收敛性、数值稳定性和算法精度等;二是要有较好的计算复杂度。

  • 数值稳定性

    • 数值稳定性即算法对舍入误差的敏感性。舍入误差对计算结果的精确性影响小的算法,具有较好的数值稳定性;反之,算法的数值稳定性差。
  • 算法精度

    • 数值算法的精度衡量了近似解对精确解的近似程度。算法的精度越高并不一定导致近似解的误差越小。精度的具体定义依赖于所考虑的问题类型。在偏微分方程的数值算法中我们一般考虑收敛精度,它指的是算法的收敛速率。
  • 计算复杂度

    • 计算复杂度指的是通过数值计算解决问题的困难程度。最常见的是时间复杂度和空间复杂度
  • 收敛性

    • 算法的收敛性一般指数值解逼近精确解的性质,而对于不同的问题,其定义也是有区别的。例如在迭代算法中,收敛性指的是近似解能够随着迭代过程趋向于精确解,它通常与迭代矩阵的谱半径相关。在偏微分方程数值算法中,收敛性指的是随网格尺寸等参数变化的数值解的收敛性质,它与算法的收敛精度密切相关。

误差理论

误差按照它们的来源可分为以下四种

  1. 模型误差:反映实际问题有关量之间关系的计算公式,即数学模型,通常只是近似的。由此产生的数学模型的解与实际问题的解之间的误差称为模型误差。

  2. 观测误差:数学模型中包含的某些参数(如时间、长度、电压等等)往往通过观测而获得。由观测得到的数据与实际的数据之间是有误差的。这种误差称为观测误差。

  3. 截断误差求解数学模型所用的数值计算方法如果是一种近似的方法,那么只能得到数学模型的近似解,由此产生的误差称为截继误差或方法误差。例如,由Taylor (泰勒)公式,函数 f(x)f(x) 可表示为

f(x)=f(0)+f(0)x+f(0)2!x2++f(n)(0)n!xn+f(n+1)(θx)(n+1)!xn+1,(0<θ<1)f(x)=f(0)+f'(0)x+\frac{f''(0)}{2!}x^2+\cdots+\frac{f^{(n)}(0)}{n!}x^n+\frac{f^{(n+1)}(\theta x)}{(n+1)!}x^{n+1},(0<\theta<1)

近似即:(此近似公式的误差就是截断误差)

f(x)f(0)+f(0)x+f(0)2!x2++f(n)(0)n!xnf(x)\approx f(0)+f'(0)x+\frac{f''(0)}{2!}x^2+\cdots+\frac{f^{(n)}(0)}{n!}x^n

  1. 舍入误差:由于计算机的字长有限,参加运算的数据以及运算结果在计算机上存放会产生误差

数值计算的误差度量

对于准确值xx、近似值 aa 和误差限 ϵ\epsilon

  • e=xae=x-a ,称 ee 为近似值 aa绝对误差,简称误差

  • 如果 e|e| 的一个上界已知,记为 ϵ\epsilon ,即 eϵ|e|\leq \epsilon 则称 ϵ\epsilon 为近似值 aa绝对误差限或绝对误差界,简称误差限或误差界

  • er=ex=xaxe_r=\frac ex=\frac{x-a}xere_r 为近似值 aa 的相对误差。由于 xx 未知,实际上总把 ea\frac e a作为 aa 的相对误差,并且也记为er=ea=xaae_r=\frac ea=\frac{x-a}a。一般用百分比表示

  • er|e_r| 的上界,即 ϵr=ϵa\epsilon_r=\frac{\epsilon}{|a|} 称为近似值 aa相对误差限或相对误差界。

  • 有效数字:由准确值经过四舍五入得到的近似值(经过对于某一绝对误差限舍入后),从它的末位数字到第一位非零数字都是有效数字

    • 严格来说:如果 aa 的绝对误差限是它的某一位的半个单位,并且从该位到它的第一位非零数字共有 nn 位,则称用 aa 近似 xx 时具有 nn 位有效数字。

PS. aϵxa+ϵa-\epsilon\leq x\leq a+\epsilon,即 x=a±ϵ.x=a\pm\epsilon.

函数值的误差估计

给定多元函数A=f(x1,x2,,xn)A=f\left(x_1,x_2,\cdots,x_n\right) ,且 x1,x2,,xnx_1^*,x_2^*,\cdots,x_n^*x1,x2,,xnx_1,x_2,\cdots,x_n 近似值,于是可求 AA 的近似值 A=f(x1,,xn)A^*=f\left(x_1^*,\cdots,x_n^*\right)

绝对误差

在点 x=(x1,x2,,xn)x=\left(x_1,x_2,\cdots,x_n\right) 进行泰勒展开,由

f(x1,x2,,xn)=f(a1,a2,,an)+i=1nfxi(xiai)+12!i,j=1n2fxixj(xiai)(xjaj)+f(x_{1}, x_{2}, \ldots, x_{n}) = f(a_{1}, a_{2}, \ldots, a_{n}) + \sum_{i=1}^{n} \frac{\partial f}{\partial x_{i}} (x_{i} - a_{i}) + \frac{1}{2!} \sum_{i,j=1}^{n} \frac{\partial^{2} f}{\partial x_{i} \partial x_{j}} (x_{i} - a_{i})(x_{j} - a_{j}) + \cdots

其中, fxi\frac{\partial f}{\partial x_{i}} 表示函数 ff 对变量 xix_{i} 的偏导数,2fxixj\frac{\partial^{2} f}{\partial x_{i} \partial x_{j}} 表示函数ff 对变量 xix_{i}xjx_{j} 的混合偏导数,依此类推。所以:

AA=f(x1,,xn)f(x1,x2,,xn)j=1nf(x)xj(xjxj)\begin{aligned}A^{*}-A&=f\left(x_{1}^{*},\cdots,x_{n}^{*}\right)-f\left(x_{1},x_{2},\cdots,x_{n}\right)\\&\approx\sum_{j=1}^n\frac{\partial f(x)}{\partial x_j}\left(x_j^*-x_j\right)\end{aligned}

PS.这个等式也可f以用中值定理直接得出

e(xi)=xjxje(x_i)=x_j^*-x_j ,也即:

e(A)j=1nf(x)xje(xj),e(A)j=1nf(x)xje(xj)e(A)\approx\sum_{j=1}^n\frac{\partial f(x)}{\partial x_j}e\left(x_j\right),|e(A)|\leqslant\sum_{j=1}^n\left|\frac{\partial f(x)}{\partial x_j}\right||e\left(x_j\right)|

PS.

相对误差

er(A)=e(A)Aj=1nxjf(x)f(x)xjer(xj)e_r(A)=\frac{e(A)}{A}\approx\sum_{j=1}^n\frac{x_j}{f(x)}\frac{\partial f(x)}{\partial x_j}e_r\left(x_j\right)

  • er(xj)=e(xj)xje_r (x_j) = \frac{e(x_j)}{x_j} ,即 e(xj)=er(xj)xje(x_j)= e_r (x_j)\cdot {x_j}
  • 因子 xjf(x)f(x)xj\frac{x_j}{f(x)}\frac{\partial f(x)}{\partial x_j} 反映 xjx_j 相对误差 er(xj)e_r (x_j) 对相对误差 er(A)e_r(A) 影响的程度

函数值误差的代数运算

  • f(x,y)=xyf(x,y)=xy
    • e(xy)ye(x)+xe(y)e(xy)\approx ye(x)+xe(y)
    • er(xy)er(x)+er(y)e_r(xy)\approx e_r(x)+e_r(y)
  • f(x,y)=x/yf(x,y)=x/y
    • e(x/y)1ye(x)xy2e(y)e(x/y)\approx\frac{1}{y}e(x)-\frac{x}{y^{2}}e(y)
    • er(x/y)er(x)er(y)e_{r}(x/y)\approx e_{r}(x)-e_{r}(y)
  • f(x,y)=x±yf(x,y)=x\pm y
    • e(x±y)e(x)±e(y)e(x\pm y)\approx e(x)\pm e(y)
    • er(x±y)xx±yer(x)±yx±yer(y)e_{r}(x\pm y)\approx\frac{x}{x\pm y}e_{r}(x)\pm\frac{y}{x\pm y}e_{r}(y)x±y0x\pm y\neq0
  • f(x)=xf(x)=\sqrt{x}
    • er(x)12er(x)e_r(\sqrt{x})\approx\frac{1}{2}e_r(x)
  • f(x)=xnf(x)=x^{n}
    • er(xn)ner(x)e_{r} (x^{n})\approx ne_{r}(x)

几点注意:

  • 对于 er(xy)xxyer(x)yxyer(y)e_{r}(x- y)\approx\frac{x}{x- y}e_{r}(x)-\frac{y}{x- y}e_{r}(y) ,如果 xxyy 非常接近,会出现很大计算误差
  • 同样的,对于 e(x/y)1ye(x)xy2e(y)e(x/y)\approx\frac{1}{y}e(x)-\frac{x}{y^{2}}e(y) ,若 y0y\to 0,同理

算法在数值计算中应注意的几个问题

  • 要有数值稳定性。即能控制舍入误差的传播
  • 两数相加要防止较小的数加不到较大的数中所引起的严重后果
  • 避免两个相近的近似值相减,以免严重损失有效数字
  • 除法运算中,要尽量避免除数的绝对值远小于被除数的绝对值

范数

向量范数

定义向量大小的量,又称为向量的模。是定义在 RnR^n 上的实值函数 \|\cdot\|,满足:

  • x0\|x\|\geq 0 且仅 x=0,x=0x=0,\|x\|=0
  • kx=kx,kR\|kx\|=|k|\cdot\|x\|,k\in R
  • x+yx+y\|x+y\|\leq\|x\|+\|y\|

证明是向量范数,只需要满足上面三个条件

三种向量范数(pp-范数)

RnR^n 中的任一向量 x=(x1,x2,,xn)Tx=(x_{1},x_{2},\cdots,x_{n})^{\mathrm{T}}

x1=i=1nxix2=i=1nxi2x=max1inxj\|x\|_{1}=\sum_{i=1}^{n}|x_{i}|\\\|x\|_{2}=\sqrt{\sum_{i=1}^{n}x_{i}^{2}}\\\|x\|_{\infty}=\max_{1\leqslant i\leqslant n}|x_{j}|

范数等价性定理

α,β\|\cdot\|_\alpha,\|\cdot\|_\betaRnR^n 上的任意两种向量范数,则存在与向量 xx 无关的常数 mmM(0<m<M)M(0 < m < M), 使下列关系成立

mxαxβMxα,xRnm\|x\|_\alpha\leqslant\|x\|_\beta\leqslant M\|x\|_\alpha,\quad\forall x\in\mathbb{R}^n

此外:

摘自 张绍飞,赵迪. 矩阵论教程[M]. 机械工业出版社,2012. p98

矩阵范数

用于定义矩阵“大小”的量,又称为矩阵的模,定义在 Rn×nR^{n\times n} 上的实值函数 \|\cdot\|,满足对任意矩阵 A,BA,B

  • A0\|A\|\geq 0 且仅 A=O,A=0A=O,\|A\|=0
  • kA=kA,kR\|kA\|=|k|\cdot\|A\|,k\in \mathbb R
  • A+BA+B\|A+B\|\leq\|A\|+\|B\|
  • ABAB\|AB\|\leq\|A\|\cdot\|B\|
矩阵范数相容性

对于给定的向量范数 \|\cdot\| 和矩阵范数 \|\cdot\| , 如果对任一个 xRnx\in \mathbb R^n 和任一个 ARn×nA\in \mathbb R^{n\times n}, 满足 AxAx\|Ax\|\leqslant\|A\|\|x\| 则称所给的矩阵范数与向量范数是相容的

当在同一个问题中需要同时使用矩阵范数和向量范数时,这两种范数应当是相容的

矩阵的算子范数定理

设在 RnR^n 中给定了一种向量范数,对任一矩阵 ARn×nA\in\mathbb{R}^{n\times n}, 令:

A=maxx=1Ax(1.3.1)\|A\|=\max_{\|x\|=1}\|Ax\| \tag{1.3.1}

则由式1.3.1定义的 \|\cdot\|是一种矩阵范数,并且它与所给定的向量范数相容

称式1.3.1所定义的矩阵范数为从属于所给定向量范数的矩阵范数,又称为矩阵的算子范数

maxx=1Ax\max_{\|x\|=1}\|Ax\| 表示矩阵 AA 作用在单位范数向量 xx 上所得到的结果的最大范数。具体来说:

  • 对于给定的矩阵 AA,我们考虑所有单位范数向量 xx,即满足 x=1\|x\|=1 的向量。
  • 我们将矩阵 AA 作用在每一个单位范数向量 xx 上,得到向量 AxAx
  • 然后计算每个 AxAx 的范数 Ax\|Ax\|
  • 最终,我们找到使得 Ax\|Ax\| 最大的单位范数向量 xx,即求解 maxx=1Ax\max_{\|x\|=1}\|Ax\|

证明为什么是一种矩阵范数:

设给定的向量范数为 p\|\cdot\|_p,则从属于向量范数 p\|\cdot\|_p 的矩阵范数仍记为 p\|\cdot\|_p,即

Ap=maxxp=1Axp\|A\|_p=\max_{\|x\|_p=1}\|Ax\|_p

矩阵 pp-范数

A=[aij]Rn×nA=[a_{ij}]\in\mathbb{R}^{n\times n},则

A1=max1jni=1naijA2=λmax(ATA)A=max1inj=1naij\begin{gathered} \|A\|_{1}=\max_{1\leqslant j\leqslant n}\sum_{i=1}^{n}|a_{ij}| \\ \|A\|_{2}=\sqrt{\lambda_{\max}\left(A^{\mathrm{T}}A\right)} \\ \|A\|_{\infty}=\max_{1\leqslant i\leqslant n}\sum_{j=1}^{n}|a_{ij}| \end{gathered}

其中 λmax(ATA)\lambda_{\mathrm{max}} (A^{\mathrm{T}}A)表示矩阵 ATAA^{\mathrm{T}}A 的最大特征值。

矩阵范数 A1,A2,A\|A\|_{1},\|A\|_{2},\|A\|_{\infty} 又分别称为矩阵的列范数、谱范数和行范数

佛罗贝尼乌斯范数

另外还有一种常用的矩阵范数,就是Frobenius(佛罗贝尼乌斯) 范数,又称为Euclid 范数,

AF=i,j=1naij2\|A\|_\mathrm{F}=\sqrt{\sum_{i,j=1}^na_{ij}^2}

其中 A=[aij]Rn×nA=[a_{ij}]\in\mathbb{R}^{n\times n},又可以记作 AE\|A\|_\mathrm{E}

可以证明,F\|\cdot\|_\mathrm{F} 与向量范数 2\|\cdot\|_{2} 相容,即

Ax2AFx2,ARn×n,xRn\|Ax\|_2\leqslant\|A\|_\text{F}\|x\|_2,\quad A\in\mathbb{R}^{n\times n},x\in\mathbb{R}^n

F\|\cdot \|_\mathrm{F} 不从属于任何向量范数,即不是算子范数

算子范数的估计

单位矩阵 II 的任何一种算子范数都有

I=maxx=1Ix=1\|I\|=\max\limits_{\|x\|=1}\|Ix\|=1


设矩阵 ARn×nA \in R^{n\times n} 的某种范数 A<1\left\|A\right\|<1 ,则 I±AI\pm A 为非奇异矩阵,并且当该种范数为算子范数时,下式成立

(I±A)111A\left\|(I\pm A)^{-1}\right\|\leqslant\frac1{1-\|A\|}

I±AI\pm A 奇异,即 (I±A)x=0(I\pm A)\cdot x=0,展开并移项,即 x=Axx=\mp Ax,取范数,即 x=AxAx\|x\|=\|Ax\| \leq \|A\|\cdot \|x\|

A1\|A\|\geq1,矛盾

当然这里还涉及到一个点,即矩阵范数相容性,即满足 AxAx\|Ax\|\leqslant\|A\|\|x\|

此外也可以通过这种方式来证明:

ei=(010)e_i=\left(\begin{array}{c}0\\1\\0\end{array}\right)

则:

xei=[0x100x200x30]x\cdot e_i=\begin{bmatrix} 0 & x_1 & 0 \\ 0 & x_2 & 0 \\ 0 & x_3 & 0 \end{bmatrix}

取范数,即 x=xei\|x\|=\|x\cdot e_i\|,相当于我们构建了一个向量范数与矩阵范数的等式

那么 Ax=AxeiAxei=Ax\|Ax\|=\|Ax\cdot e_i\|\leq\|A\|\cdot\|x\cdot e_i\|=\|A\|\cdot\|x\|

(IA)(IA)1=I(I-A)(I-A)^{-1}=I ,则 (IA)1A(IA)1=I(I-A)^{-1}-A\cdot (I-A)^{-1}=I,即 (IA)1=A(IA)1+I(I-A)^{-1}=A\cdot (I-A)^{-1}+I

两边同取算子范数:(IA)1=A(IA)1+IA(IA)1+I\|(I-A)^{-1}\|=\|A\cdot (I-A)^{-1}\|+\|I\|\leq\|A\|\cdot\|(I-A)^{-1}\|+\|I\|

移项,即 (1A)(IA)1I=1(1-\|A\|)\|(I-A)^{-1}\|\leq\|I\|=1,而 A<1\|A\|<1,故上式得证

PS.一个方阵是奇异的,如果它的行列式为零。非奇异矩阵则是行列式非零的方阵。

  1. 逆矩阵:奇异矩阵是不可逆的,因为逆矩阵的存在要求矩阵的行列式不为零。因此,奇异矩阵没有逆矩阵。

  2. :奇异矩阵的秩小于其阶数。一个 n×nn \times n 的奇异矩阵的秩最多为 n1n-1

  3. 零空间:奇异矩阵的零空间(也称为核)不仅包含零向量,还包含其他非零向量。这是因为奇异矩阵将某些向量映射到零向量。

  4. 特征值:奇异矩阵的特征值中至少有一个为零。这是因为特征值是行列式与矩阵迹的根,而奇异矩阵的行列式为零。

  5. 解的存在性:对于线性方程组 Ax=bAx = b,其中 AA 是奇异矩阵,解可能不存在,或者有无穷多个解。这取决于向量 bb 是否位于矩阵 AA 的列空间中。

线性方程组的解法

主要讨论如下形式的线性方程组的求解问题.

{a11x1+a12x2++a1nxn=b1,a21x1+a22x2++a2nxn=b2,an1x1+an2x2++annxn=bn,(2.1.1)\left.\left\{\begin{array}{c}a_{11}x_1+a_{12}x_2+\cdots+a_{1n}x_n=b_1,\\a_{21}x_1+a_{22}x_2+\cdots+a_{2n}x_n=b_2,\\\vdots\\a_{n1}x_1+a_{n2}x_2+\cdots+a_{nn}x_n=b_n,\end{array}\right.\right.\tag{2.1.1}

也即:

Ax=bAx=b

其中

A=(a11a1nan1ann)Rn×n,x=(x1xn)Rn,b=(b1bn)Rn.\left.A=\left(\begin{array}{ccc}a_{11}&\cdots&a_{1n}\\\vdots&&\vdots\\a_{n1}&\cdots&a_{nn}\end{array}\right.\right)\in\mathbb{R}^{n\times n},\quad x=\left(\begin{array}{c}x_1\\\vdots\\x_n\end{array}\right)\in\mathbb{R}^n,\quad b=\left(\begin{array}{c}b_1\\\vdots\\b_n\end{array}\right)\in\mathbb{R}^n.

设系数矩阵 AA 非奇异,即 det(A)0det(A) \neq 0,则方程组有唯一解向量 xx.

Gauss 消去法

由消元和回代两部分组成

消元即是对增广矩阵 [A,b][A,b] 做初等行变换,经过 n1n-1 次消元得到

(A(n),b(n))=(a11(1)a1n(1)b1(1)a22(2)a2n(2)b2(2)Oann(n)bn(n))\left.\left(A^{(n)},b^{(n)}\right)=\left(\begin{array}{ccccc}a_{11}^{(1)}&\cdots&\cdots&a_{1n}^{(1)}&b_{1}^{(1)}\\&a_{22}^{(2)}&\cdots&a_{2n}^{(2)}&b_{2}^{(2)}\\&&\ddots&\vdots&\vdots\\O&&&a_{nn}^{(n)}&b_{n}^{(n)}\end{array}\right.\right)

即线性方程组:

{a11(1)x1+a12(1)x2++a1n(1)xn=b1(1),a22(2)x2++a2n(2)xn=b2(2),ann(n)xn=bn(n).\begin{aligned}\begin{cases}a_{11}^{(1)}x_1+a_{12}^{(1)}x_2+\cdots+a_{1n}^{(1)}x_n=b_1^{(1)},\\a_{22}^{(2)}x_2+\cdots+a_{2n}^{(2)}x_n=b_2^{(2)},\\\vdots\\a_{nn}^{(n)}x_n=b_n^{(n)}.\end{cases}\end{aligned}

通过逐步回代,依次求出 xn1,xn2,,x1x_{n-1},x_{n-2},\ldots,x_1

顺序Gauss 消去法

消元过程的第 kk 步,对矩阵需做 (nk)2(n-k)^2 次乘法运算及 (nk)(n - k) 次除法运算,对右端向量需作 (nk)(n - k) 次乘法运算,所以消去过程总的乘除法运算工作量为

k=1n1(nk)2+k=1n1(nk)+k=1n1(nk)=13n3+12n256n.\sum_{k=1}^{n-1}(n-k)^2+\sum_{k=1}^{n-1}(n-k)+\sum_{k=1}^{n-1}(n-k)=\frac13n^3+\frac12n^2-\frac56n.

回代过程中,计算每个 xkx_k 需作 (nk+1)(n - k + 1) 次乘除法运算,其工作量为

k=1n(nk+1)=12n(n+1).\sum_{k=1}^n(n-k+1)=\frac12n(n+1).

算法运行条件

要使得顺序消去法可运行,必须使得前 n1n - 1 个主元素 akk(k)(k=12n1)a^{(k)}_{kk} (k = 1,2, \dots,n - 1) 均不为零,其充要条件是系数矩阵 AA 的前 n1n-1 个顺序主子式不为零,即

Dk=a11(1)a1k(1)ak1(1)akk(1)0,k=1,2,,n1.\left.D_k=\left|\begin{array}{ccc}a_{11}^{(1)}&\cdots&a_{1k}^{(1)}\\\vdots&&\vdots\\a_{k1}^{(1)}&\cdots&a_{kk}^{(1)}\end{array}\right.\right|\neq0,\quad k=1,2,\cdots,n-1.

为什么运行条件是主元素不为零充要是系数主子式不为零,这个其实也可以用LU分解的思路来思考

当主元素 akk(k)a^{(k)}_{kk} 均不为零时,我们可以考虑使用LU分解的思路来证明对应系数矩阵的顺序主子式不为零。

假设我们有一个 n×nn \times n 的方阵 AA,其顺序主子式都不为零,即 Ak0|A_k| \neq 0,其中 AkA_k 表示 AA 的第 kk 阶顺序主子式。

我们知道LU分解可以写为 A=LUA = LU,其中 LL 是单位下三角矩阵,UU 是上三角矩阵。我们可以将 AA 分解为:

A=[a1100a21a220an1an2ann]=[100l2110ln1ln21][u11u12u1n0u22u2n00unn]A = \begin{bmatrix} a_{11} & 0 & \cdots & 0 \\ a_{21} & a_{22} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ a_{n1} & a_{n2} & \cdots & a_{nn} \end{bmatrix} = \begin{bmatrix} 1 & 0 & \cdots & 0 \\ l_{21} & 1 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ l_{n1} & l_{n2} & \cdots & 1 \end{bmatrix} \begin{bmatrix} u_{11} & u_{12} & \cdots & u_{1n} \\ 0 & u_{22} & \cdots & u_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & u_{nn} \end{bmatrix}

我们可以看到,对于LU分解中的 UU 矩阵,其对角线上的元素 ukku_{kk} 对应于原矩阵 AA 的主元素 akk(k)a^{(k)}_{kk}。由于主元素均不为零,因此LU分解中的 UU 矩阵的对角线上的元素 ukku_{kk} 也均不为零。

现在我们来看LU分解中的 LL 矩阵。由LU分解的性质可知,LL 的对角线上的元素均为1,而且 LL 的第 kk 行第 jj 列的元素 lkjl_{kj} 是原矩阵 AA 的第 jj 阶顺序主子式与 UU 的第 kk 阶顺序主子式的商。由于 AA 的顺序主子式都不为零,因此 LL 的元素 lkjl_{kj} 也不为零。

一个问题

前面我们提到,除法运算中,要尽量避免除数的绝对值远小于被除数的绝对值

如果对角线元素 aiia_{ii} 很小(导致 mik|m_{ik}| 很大),就会导致很大的误差

列主元素Gauss 消去法

kk 次消元之前,通过行变换将 aik(k)(i=kk+1,n)a^{(k)}_{ik} (i =k,k + 1,\dots ,n) 中绝对值最大的元素交换到第 kk 行的主对角线位置上,然后进行消元计算,即:

其他过程同顺序法,算法要求各个列的主元素不为零

对于列主元素交换,相当于多乘了一个初等行变换矩阵,所以最终还是可以化为一个LU分解的形式

比如乘以

P=[010100001]P = \begin{bmatrix} 0 & 1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \end{bmatrix}

交换一二行

LU分解及其与Gauss消去的关系

考试不考,仅供拓展思路

对于 Ax=bAx=bAA 分解为 A=LUA=LU,其中 LL 是下三角矩阵, UU 是上三角矩阵,则可将原方程化为:

Ly=b,Ux=yLy=b,\,\,Ux=y

先由 Ly=bLy=b 解出 yy,再由 Ux=yUx=y 解出 xx

而Gauss消去也可以看作一个LU分解的过程,即:

L(n1)L(n2)L1A=A(n)L^{(n-1)}L^{(n-2)}\cdots L^1 \cdot A=A^{(n)}

其中 L1L^1 对应于将第一列非主对角线元素都化为零的过程,即

L1=[1000l1100l20100ln1001]L^1=\begin{bmatrix} 1 & 0 & 0 & \cdots &0\\ l_1 & 1 & 0 & \cdots &0\\ l_2 & 0 & 1 & \cdots &0\\ \vdots & \vdots & \vdots &&0\\ l_{n-1} & 0 & 0 & \cdots &1\\ \end{bmatrix}

其中 lil_i 对应 mi1-m_{i1} 。以此类推

不妨设 L=(L(n1)L(n2)L1)1,U=A(n)L=(L^{(n-1)}L^{(n-2)}\cdots L^1)^{-1},U=A^{(n)},显然,LL 是下三角矩阵,AA 为上三角矩阵,即为 A=LUA=LU

矩阵条件数与病态方程组

AAΔARn×n\Delta A\in\mathbb{R}^n\times nbbΔbRn\Delta b\in\mathbb{R}^nAA非奇异,b0b\neq0xx 是方程组 Ax=bAx=b 的解向量。若 ΔA<1A1\|\Delta A\|<\frac1{\|{A}^{-1}\|} ,则有

  • 方程组

    (A+ΔA)(x+Δx)=b+Δb(A+\Delta A)(x+\Delta x)=b+\Delta b

    有唯一解 x+Δxx+\Delta x

  • 下列估计式成立:

    ΔxxAA11A1ΔA(ΔAA+Δbb)\frac{\parallel\Delta x\parallel}{\parallel x\parallel}\leqslant\frac{\parallel A\parallel\parallel A^{-1}\parallel}{1-\parallel A^{-1}\parallel\parallel\Delta A\parallel}(\frac{\parallel\Delta A\parallel}{\parallel A\parallel}+\frac{\parallel\Delta b\parallel}{\parallel b\parallel})

对于第一点,我们需要证明 A+ΔAA+\Delta A 非奇异,证明关键其实是该定理:

设矩阵 ARn×nA \in R^{n\times n} 的某种范数 A<1\left\|A\right\|<1 ,则 I±AI\pm A 为非奇异矩阵,并且当该种范数为算子范数时,下式成立

(I±A)111A\left\|(I\pm A)^{-1}\right\|\leqslant\frac1{1-\|A\|}

A+ΔA=A(I+A1ΔA)A+\Delta A=A(I+A^{-1}\Delta A),其中 AA 非奇异,即证 I+A1ΔAI+A^{-1}\Delta A 非奇异,由上定理,即证 A1ΔA<1\|A^{-1}\Delta A\|<1

A1ΔAA1ΔA<1\parallel A^{-1}\Delta A\parallel\leqslant\parallel A^{-1}\parallel\parallel\Delta A\parallel<1,故得证

对于第二点,对原方程组进行展开,并代入 Ax=bAx=b

Ax+AΔx+ΔAx+ΔAΔx=b+ΔbAΔx+ΔAx+ΔAΔx=ΔbΔx=A1ΔbA1ΔAxA1ΔAΔxAx+A\cdot \Delta x+\Delta A \cdot x+\Delta A\cdot \Delta x=b+\Delta b\\ A\cdot \Delta x+\Delta A \cdot x+\Delta A\cdot \Delta x=\Delta b\\ \Delta{x}={A}^{-1}\Delta{b}-{A}^{-1}\Delta{A}{x}-{A}^{-1}\Delta{A}\Delta{x}

两边同时取范数,由矩阵范数相容性,即 b=AxAx\|b\|=\|Ax\|\leqslant\|A\|\|x\|

Δx=A1Δb+A1ΔAx+A1ΔAΔxA1Δb+A1ΔAx+A1ΔAΔx\begin{aligned} \|\Delta x\| &= \|A^{-1}\Delta b\|+\|A^{-1}\Delta A\cdot x\|+\|A^{-1}\Delta{A}\Delta{x}\|\\ &\leq \parallel A^{-1}\parallel\parallel\Delta b\parallel+\parallel A^{-1}\parallel\parallel\Delta A\parallel\parallel x\parallel+\parallel A^{-1}\parallel\parallel\Delta A\parallel\parallel\Delta x\parallel \end{aligned}

(1A1ΔA)ΔxxAA1Δbb+A1ΔA(1-\parallel A^{-1}\parallel\parallel\Delta A\parallel)\frac{\parallel\Delta x\parallel}{\parallel x\parallel}\leqslant\parallel A\parallel\parallel A^{-1}\parallel\frac{\parallel\Delta b\parallel}{\parallel b\parallel}+\parallel A^{-1}\parallel\parallel\Delta A\parallel

由第一点 A1ΔA<1\parallel A^{-1}\parallel\parallel\Delta A\parallel<1,上式得证

一个细节就是这里用到的都是算子范数,使得 矩阵范数与向量范数相容,从而满足范数相容性

对于上面定理第二点的式子,我们可以改写成:

ΔxxAA11AA1ΔAA(ΔAA+Δbb)\frac{\parallel\Delta x\parallel}{\parallel x\parallel}\leqslant\frac{\parallel A\parallel\parallel A^{-1}\parallel}{1-\parallel A\parallel\parallel A^{-1}\parallel\frac{\parallel\Delta A\parallel}{\parallel A\parallel}}(\frac{\parallel\Delta A\parallel}{\parallel A\parallel}+\frac{\parallel\Delta b\parallel}{\parallel b\parallel})

那么 AA1\parallel A\parallel\parallel A^{-1}\| 越小,A,bA,b 的误差对解向量的相对误差影响就越小,反之则大。

条件数与病态

关于病态方程组求解

对于病态方程组的求解:

  • 用更高精度的运算
  • 平衡方法:AA 元素数量级差别很大时进行行平衡或列平衡
  • 残差矫正法:病态但不特别严重时的一种缓解病态方法

迭代法

迭代法是目前求解大规模稀疏线性方程组的主要方法之一.

迭代法的一般形式及其收敛性

对于线性方程组,将系数矩阵分解为两个矩阵的差,即 A=NPA=N-P ,代入 Ax=bAx=b,即:

Nx=Px+bx=N1Px+N1bNx=Px+b\\x=N^{-1}Px+N^{-1}b

G=N1P,d=N1bG=N^{-1}P,\quad d=N^{-1}b 代入,得 x=Gx+dx=Gx+d ,取 x(0)Rn×nx^{(0)}\in\mathbb{R}^{n\times n} 作为初始近似解,使用递推式

x(k+1)=Gx(k)+d(k=0,1,)x^{(k+1)}=Gx^{(k)}+d\quad(k=0,1,\cdots)

产生一个向量序列,当 kk 足够大时, xkx^k 可以看作原方程的近似解

序列收敛性

只有在迭代法收敛的情况下,用它所产生的向量序列 {x(k)}\{x^{(k)}\} 中的向量作为方程组的近似解才有意义,而且,k 越大,x(k)x^{(k)} 作为方程组的解就越精确。关于向量序列收敛定义如下:

设有向量序列

x(k)=(x1(k),x2(k),,xn(k))(k=0,1,)x^{(k)}=\left({x_1}^{(k)},{x_2}^{(k)},\cdots,{x_n}^{(k)}\right)^\top\quad(k=0,1,\cdots)

如果存在常向量x=(x1,x2,,xn)x^*=(x_1^*,x_2^*,\cdots,x_n^*)^\top,使得

limkxi(k)=xi(i=1,2,,n)\lim\limits_{k\to\infty}x_i^{(k)}=x_i^*\quad(i=1,2,\cdots,n)

则称向量序列{x(k)}\{x^(k)\}收敛于常向量xx^*,记为

limkx(k)=x\lim\limits_{k\to\infty}x^{(k)}=x^*

相关定理

设有向量序列{x(k)}\{x^{(k)}\}和常向量 xx^*,如果对某种范数,有

limkx(k)x=0\lim\limits_{k\to\infty}\left\|x^{(k)}-x^*\right\|=0

则必有

limkx(k)=x\lim\limits_{k\to\infty}x^{(k)}=x^*

由向量范数等价性,原条件即:

limkmax1inxi(k)xi=0\lim_{k\to\infty}\max_{1\leqslant i\leqslant n}\left|x_i^{(k)}-x_i^*\right|=0

一个后面用到的等式

x(k)x=i=k(x(i)x(i+1))x^{(k)}-x^*=\sum_{i=k}^\infty\left(x^{(i)}-x^{(i+1)}\right)

要证明这个公式,我们从Jacobi迭代法的定义出发。Jacobi迭代法用于求解线性方程组 Ax=bAx = b,迭代公式为:

x(k+1)=D1(b(L+U)x(k))x^{(k+1)} = D^{-1}(b - (L + U)x^{(k)})

其中 A=D+L+UA = D + L + UDDAA 的对角矩阵,LLAA 的严格下三角矩阵,UUAA 的严格上三角矩阵。

我们要证明:x(k)x=i=k(x(i)x(i+1))x^{(k)} - x^* = \sum_{i=k}^\infty (x^{(i)} - x^{(i+1)})

首先,我们知道对于任何迭代步 kk,有以下关系:

x(k)x=(x(k)x(k+1))+(x(k+1)x)x^{(k)} - x^* = (x^{(k)} - x^{(k+1)}) + (x^{(k+1)} - x^*)

利用这一点,我们可以递推:

x(k+1)x=(x(k+1)x(k+2))+(x(k+2)x)x^{(k+1)} - x^* = (x^{(k+1)} - x^{(k+2)}) + (x^{(k+2)} - x^*)

这样继续下去,我们得到:

x(k)x=(x(k)x(k+1))+(x(k+1)x(k+2))+(x(k+2)x(k+3))+x^{(k)} - x^* = (x^{(k)} - x^{(k+1)}) + (x^{(k+1)} - x^{(k+2)}) + (x^{(k+2)} - x^{(k+3)}) + \cdots

在极限的意义下:

x(k)x=i=k(x(i)x(i+1))x^{(k)} - x^* = \sum_{i=k}^\infty (x^{(i)} - x^{(i+1)})

这个公式成立的关键在于迭代过程的收敛性,即 x(i)xx^{(i)} \to x^*ii \to \infty。因此,尾项 (x(i)x)(x^{(i)} - x^*) 在极限下趋于零。

谱半径与收敛性

n×nn\times n 矩阵 GG 的特征值是 λ1,λ2,,λn\lambda_1,\lambda_2,\cdots,\lambda_n ,称

ρ(G)=max1inλi\rho(G)=\max_{1\leqslant i\leqslant n}|\lambda_i|

为矩阵 GG 的谱半径.

谱半径与收敛性的定理
  • 定理:任意一个矩阵的谱半径与其范数是有关系的,即任一矩阵AA的谱半径不大于AA的任一与某一向量范数相容的矩阵范数,即

ρ(A)A.\rho(A)\leqslant\|A\|.

按照谱半径的定义 ρ(A)=max1inλi\rho(A)=\max_{1\leq i\leq n}|\lambda_i|, 其中,λi(i=1,2,,n)\lambda_i(i=1,2,\cdots,n)nn 阶方阵 AA 的特征值

xx 为对应于 λ\lambdaAA 的特征向量,则有 λx=Ax\lambda x=Ax,两边取范数可得

λxAx.|\lambda|\|x\|\leqslant\|A\|\|x\|.

因为 xx 为非零向量,x0\|x\|\neq0,故有λA.|\lambda|\leqslant\|A\|.

上述这个不等式对AA的任何特征值均成立,即可得 ρ(A)A.\rho(A)\leqslant\|A\|.

  • 对上定理,利用矩阵的Jordan 标准型得:设 AA 是任意 nn 阶矩阵, 则 AAkk 次幂 Ak0A^k\to 0 ( 当 k+k \to +\infty) 的充要条件为谱半径 ρ(A)<1\rho(A)<1

Jordan标准型是一种特殊形式的矩阵,通常用于简化矩阵的分析和计算特征值、特征向量等问题。对于一个给定的方阵,其Jordan标准型可以通过相似变换得到。以下是关于Jordan标准型的一些重要概念:

  1. 定义:对于一个 n×nn \times n 的方阵 AA,存在一个可逆矩阵 PP,使得 P1APP^{-1}AP 是Jordan标准型。Jordan标准型是一个块对角矩阵,其中每个块称为Jordan块。
  2. Jordan块:Jordan块是一个形式为 λI+N\lambda I + N 的矩阵,其中 λ\lambda 是矩阵的特征值,II 是单位矩阵,NN 是上三角矩阵,其对角线元素为0或者都是1。Jordan块的大小取决于特征值的代数重数和几何重数。
  3. 特征值和特征向量:Jordan标准型将矩阵的特征值和特征向量整合在一起,使得特征向量更易于计算和理解。每个Jordan块对应一个特征值,而特征向量则可以从Jordan块中推导出来。
  • 定理:对任意的向量 dd, 迭代法收敛的充分必要条件是 ρ(G)<1\rho(G)<1

对迭代形式 x(k)=Gx(k1)+dx^{(k)}=Gx^{(k-1)}+d 和最终解 x=Gx+dx^*=Gx^*+d 做差并递推:

x(k)x=G(x(k1)x)==Gk(x(0)x)x^{(k)}-x^*=G\left(x^{(k-1)}-x^*\right)=\cdots=G^{k}\left(x^{(0)}-x^*\right)

若收敛,则 x(k)x0x^{(k)}-x^*\to 0 ,即 Gk(x(0)x)0G^{k}\left(x^{(0)}-x^*\right)\to 0,即 Gk0,kG^{k}\to 0,k\to \infty,由上引理得 ρ(G)<1\rho(G)<1

同理可证充分性

收敛性判断

用迭代矩阵的谱半径判断迭代公式是否收敛往往不容易,可以使用如下定理

如果矩阵 GG 的某种范数G<1\|G\|<1,则

  • 方程组 x=Gx+dx=Gx+d 的解 xx^* 存在且唯一;

  • 对于迭代公式 x(k+1)=Gx(k)+dx^{(k+1)}=Gx^{(k)}+d

    limkx(k)=x,x(0)Rn\lim\limits_{k\to\infty}x^{(k)}=x^*,\quad\forall x^{(0)}\in\mathbb{R}^n

并且下列两式成立

x(k)xGk1Gx(1)x(0)(2.3.6)\left\|x^{(k)}-x^*\right\|\leqslant\frac{\|G\|^k}{1-\|G\|}\left\|x^{(1)}-x^{(0)}\right\| \tag{2.3.6}

x(k)xG1Gx(k)x(k1)(2.3.7)\left\|x^{(k)}-x^*\right\|\leqslant\frac{\|G\|}{1-\|G\|}\left\|x^{(k)}-x^{(k-1)}\right\|\tag{2.3.7}

对于第一点,因 G<1\|G\| < 1 , 根据定理算子范数的估计,可知矩阵 IGI-G非奇异, 其中 II 是单位矩阵

对于第二点,由 G<1\|G\| < 1 ,则 ρ(G)<1\rho(G)<1 ,即迭代法收敛,则 limkx(k)=x\lim\limits_{k\to\infty}x^{(k)}=x^*

则对 x(k)x=i=k(x(i)x(i+1))x^{(k)}-x^*=\sum_{i=k}^\infty\left(x^{(i)}-x^{(i+1)}\right) 两边取范数,由x+yx+y\|x+y\|\leq\|x\|+\|y\| 得:

x(k)xi=kx(i)x(i+1)=i=kGi(x(0)x(1))i=kGix(0)x(1)=Gk1Gx(1)x(0),\begin{aligned}\left\|x^{(k)}-x^*\right\|&\leqslant\sum_{i=k}^\infty\left\|x^{(i)}-x^{(i+1)}\right\|=\sum_{i=k}^\infty\left\|G^i\left(x^{(0)}-x^{(1)}\right)\right\|\\&\leqslant\sum_{i=k}^\infty\|G\|^i\left\|x^{(0)}-x^{(1)}\right\|\\&=\frac{\|G\|^k}{1-\|G\|}\left\|x^{(1)}-x^{(0)}\right\|,\end{aligned}

PS. x(k)x=i=k(x(i)x(i+1))x^{(k)}-x^*=\sum_{i=k}^\infty\left(x^{(i)}-x^{(i+1)}\right) 这个等式上面证过

对于2.3.7

x(k)x=G(x(k1)x)=G(x(k1)x(k))+G(x(k)x).x^{(k)}-x^*=G\left(x^{(k-1)}-x^*\right)=G\left(x^{(k-1)}-x^{(k)}\right)+G\left(x^{(k)}-x^*\right).

则:

x(k)xGx(k1)x(k)+Gx(k)x,x(k)xG1Gx(k)x(k1).\begin{aligned}\left\|x^{(k)}-x^*\right\|&\leqslant\|G\|\left\|x^{(k-1)}-x^{(k)}\right\|+\|G\|\left\|x^{(k)}-x^*\right\|,\\&\left\|x^{(k)}-x^*\right\|\leqslant\frac{\|G\|}{1-\|G\|}\left\|x^{(k)}-x^{(k-1)}\right\|.\end{aligned}

Jacobi 迭代法

若系数矩阵 A=(aij)Rn×nA=(a_{ij})\in\mathbb{R}^{n\times n} 满足条件 aii0(i=1,2,,n)a_{ii}\neq0(i=1,2,\cdots,n),则将 AA 分解为 A=D+L+UA=D+L+U

其中 DD 为对角矩阵,LL 为严格下三角矩阵,UU 为严格上三角矩阵

相当于 D=diag(aii)D=diag(a_{ii}) 且显然 D1D^{-1} 存在

对于一般形式,取 N=D,P=(L+U)N=D,P=-(L+U) ,得迭代式:

x(k+1)=D1(L+U)x(k)+D1b(k=0,1,)x^{(k+1)}=-D^{-1}(L+U)x^{(k)}+D^{-1}b\quad(k=0,1,\cdots)

其中 x(0)Rnx^{(0)}\in\mathbb{R}^n 任取,迭代矩阵:

GJ=D1(L+U)G_J=-D^{-1}(L+U)

此外,d=D1bd=D^{-1}b

D1=diag(aii1)D^{-1}=diag(a_{ii}^{-1}),分量形式为:

xi(k+1)=(j=1naijxj(k)+bi)/aii(i=1,2,,n;k=0,1,)xi(k+1)=j=1naijaiixj(k)+biaii(i=1,2,,n;k=0,1,)x_i^{(k+1)}=(-\sum_{j=1}^na_{ij}x_j^{(k)}+b_i)/a_{ii}\quad(i=1,2,\ldots,n;k=0,1,\ldots)\\ x_i^{(k+1)}=-\sum_{j=1}^n\frac{a_{ij}}{a_{ii}}x_j^{(k)}+\frac{b_i}{a_{ii}}\quad(i=1,2,\ldots,n;k=0,1,\ldots)

  • 由定理:对任意的向量 dd, 迭代法收敛的充分必要条件是 ρ(G)<1\rho(G)<1,得: Jacobi 迭代法收敛的充分必要条件是 ρ(GJ)<1\rho(G_J)<1

  • 由定理:矩阵 GG 的某种范数G<1\|G\|<1,则方程组 x=Gx+dx=Gx+d 的解 xx^* 存在且唯一,得 :如果 GJ<1\|G_{J}\|<1, 则 Jacobi 迭代法收敛.

引理:

若矩阵 ARn×nA\in\mathbb{R}^{n\times n} 是主对角线按行 (或按列) 严格占优阵,则 AA 是非奇异矩阵.

主对角线按行 (或按列) 严格占优阵,即满足:

aii>j=1jinaij,i=1,2,,n.|a_{ii}|>\sum_{j=1\atop j\neq i}^n|a_{ij}|\:,\quad i=1,2,\cdots,n.

主对角线元素大于该行非对角线元素绝对值之和

如果系数矩阵AA是主对角线按行(或按列) 严格占优阵,则用Jacobi迭代法求解必收敛.

AA对称正定矩阵,则雅可比迭代收敛的充要条件是 2DA2D - A 也是对称正定矩阵.

Gauss-Seidel迭代法

在迭代法一般形式中,取N=D+L,P=UN=D+L,P=-U,形成以下的迭代公式

x(k+1)=(D+L)1Ux(k)+(D+L)1b(k=0,1,)x^{(k+1)}=-(D+L)^{-1}Ux^{(k)}+(D+L)^{-1}b\quad(k=0,1,\cdotp\cdotp\cdotp)

其中:

GG=(D+L)1Ud=(D+L)1bG_G=-(D+L)^{-1}U\\ d=(D+L)^{-1}b

使用时为了避免求 D+LD+L 的逆矩阵,可以通过这种方式改写:

(D+L)x(k+1)=Ux(k)+bDx(k+1)=Lx(k+1)Ux(k)+bx(k+1)=D1Lx(k+1)D1Ux(k)+D1b\begin{aligned}&(D+L)x^{(k+1)}=-Ux^{(k)}+b\\&Dx^{(k+1)}=-Lx^{(k+1)}-Ux^{(k)}+b\\&x^{(k+1)}=-D^{-1}Lx^{(k+1)}-D^{-1}Ux^{(k)}+D^{-1}b\end{aligned}