数值计算day4-求解线性方程组(下)
上节课主要介绍了线性方程组的几种直接求解法，包括高斯消去法（主元消去）、高斯约当法(可以用来求解矩阵的逆)、LU分解法（高斯消去和Crout 方法），使用高斯消去法实现LU分解的原理可以通过一种拉格朗日形式来理解。本节课主要介绍线性方程组的几种迭代求解法，包括Jacobi法和高斯塞德尔法。
1. 使用MATLAB内置函数实现LU分解以及矩阵求逆
LU分解
MATLAB使用lu命令对矩阵\(A\)进行LU分解，分解前对矩阵\(A\)进行行变换，以满足主元消去的条件：
\(P A= LU\) \(A=(P^{-1}L)U\)
>> help lu
lu - lu 矩阵分解
此 MATLAB 函数 返回矩阵 Y，其中包含严格下三角矩阵 L（即不带其单位对角线）和上三角矩阵 U 作为子矩阵。也即，如果 [L,U,P] =
lu(A)，Y = U+L-eye(size(A))。不会返回置换矩阵 P。
Y = lu(A)
[L,U] = lu(A)
[L,U,P] = lu(A)
[L,U,P,Q] = lu(A)
[L,U,P,Q,R] = lu(A)
[...] = lu(A,'vector')
[...] = lu(A,thresh)
[...] = lu(A,thresh,'vector')
另请参阅 cond, det, ilu, inv, qr, rref
lu 的参考页
名为 lu 的其他函数
>> A = [1 2 3;4 5 6;7 8 9]
A =
1 2 3
4 5 6
7 8 9
>> format long
>> det(A)
ans =
-9.516197353929915e-16
>> [l u p] = lu(A)
l =
1.000000000000000 0 0
0.142857142857143 1.000000000000000 0
0.571428571428571 0.500000000000000 1.000000000000000
u =
7.000000000000000 8.000000000000000 9.000000000000000
0 0.857142857142857 1.714285714285714
0 0 -0.000000000000000
p =
0 0 1
1 0 0
0 1 0
>> p*A
ans =
7 8 9
1 2 3
4 5 6
>> l*u
ans =
7 8 9
1 2 3
4 5 6
>> [l_1,u_1] = lu(A);
>> inv(p)*l == l_1
ans =
3×3 logical 数组
1 1 1
1 1 1
1 1 1
>> [l_2,u_2,p_2] = lu(p*A)
l_2 =
1.000000000000000 0 0
0.142857142857143 1.000000000000000 0
0.571428571428571 0.500000000000000 1.000000000000000
u_2 =
7.000000000000000 8.000000000000000 9.000000000000000
0 0.857142857142857 1.714285714285714
0 0 -0.000000000000000
p_2 =
1 0 0
0 1 0
0 0 1
矩阵求逆
MATLAB使用inv命令求矩阵的逆，但当矩阵的数量级很小时，inv求解会带来较大的舍入误差（Round-off error）:
>> A = 10^(-4)*[1 2 3;4 5 6;7 8 9]
A =
1.0e-03 *
0.100000000000000 0.200000000000000 0.300000000000000
0.400000000000000 0.500000000000000 0.600000000000000
0.700000000000000 0.800000000000000 0.900000000000000
>> inv(A)
警告: 矩阵接近奇异值，或者缩放错误。结果可能不准确。RCOND = 1.317607e-17。
ans =
1.0e+19 *
-0.527049830677416 1.054099661354831 -0.527049830677416
1.054099661354831 -2.108199322709663 1.054099661354832
-0.527049830677415 1.054099661354831 -0.527049830677416
>> A*inv(A)
警告: 矩阵接近奇异值，或者缩放错误。结果可能不准确。RCOND = 1.317607e-17。
ans =
1.250000000000000 0 0
0.500000000000000 1.000000000000000 -0.500000000000000
1.000000000000000 0 0
2. Jacobi法求解线性方程组
Jacobi法是一种同步迭代求解法，考虑如下线性方程构成的系统：$$\begin{cases}a_{11}x_1+\cdots+a_{1n}x_n = b_1\a_{21}x_1+\cdots+a_{2n}x_n = b_2\ \vdots \ a_{n1}x_1+\cdots+a_{nn}x_n = b_n\end{cases}$$ 简单运算，可推导出：
- \(B_n=\frac{B_{n}-b_nB_{n-1}}{d_n-b_na_{n-1}}\)
- 使用回代法求解：\(x_n=B_n\), \(x_i=B_i-a_ix_{i+1},i=n-1,n-2,...,1\)
MATLAB实现：
function x = Tridiagonal(A,B)
% The function solve a tridiagonal system of linear equations [a][x]=[b]
% using Thomas algorithm.
% Input variables:
% A The matrix of coefficients.
% B A column vector of constants.
% Output variable:
% x A colum vector with the solution.
[nR, nC] = size(A);
for i = 1:nR
d(i) = A(i,i);
end
for i = 1:nR-1
ad(i) = A(i,i+1);
end
for i = 2:nR
bd(i) = A(i,i-1);
end
ad(1) = ad(1)/d(1);
B(1) = B(1)/d(1);
for i = 2:nR-1
ad(i) = ad(i)/(d(i)-bd(i)*ad(i-1));
B(i)=(B(i)-bd(i)*B(i-1))/(d(i)-bd(i)*ad(i-1));
end
B(nR)=(B(nR)-bd(nR)*B(nR-1))/(d(nR)-bd(nR)*ad(nR-1));
x(nR,1) = B(nR);
for i = nR-1:-1:1
x(i,1) = B(i)-ad(i)*x(i+1);
end
5. 范数与条件数
5.1 向量范数
一般来说，向量范数需要满足一下三个条件：
- 正定性 $$|x|\geq 0,\forall x\in R^{n}, |x|=0\text{当且仅当}x=0$$
- 齐次性 $$|kx|=k|x|,\forall k\in R, x\in R^{n}$$
- 三角不等式 $$|x+y|\leq |x|+|y|,\forall x,y\in R^{n} $$
几种常见的向量范数：
- 1范数： $$|x|1=\sum^n|x_i|$$
- 2范数： $$|x|_2=\sqrt{x2_1+...+x2_n}$$
- \(\infty\)范数： $$|x|{\infty}=\max|x_i|$$
5.2 矩阵\(p\)范数
矩阵范数要比向量范数多一个条件：$$|AB|\leq |A||B|$$
矩阵范数中的\(p\)范数均是由向量范数诱导而来：$$|A|_p=max\frac{|Ax|_p}{|x|p}=\max|Ax|_p$$
- 1范数 $$|A|1=\max|Ax|1=\max|Ae_i|_1\text{（列绝对值和中的最大值）}$$ 证明：对任意\(\|x\|_1=1\)，注意$$x=k_1e_1+...+k_ne_n$$ 因此 $$|k_1|+...|k_n|=1$$ 根据向量范数的三角不等式，可得 $$|Ax|_1=|k_1Ae_1+...+k_nAe_n|_1\leq |k_1||Ae_1|_1+...|k_n||Ae_n|1\leq (|k_1|+...+|k_n|)\max|Ae_i|1=\max|Ae_i|_1$$ 即$$|A|1\leq \max|Ae_i|_1 $$ 当\(x=e_i\)时，能够取等号，因此$$|A|1= max|Ae_i|_1 $$
- 2范数 $$|A|2=\max{|x|2=1}|Ax|2=\sqrt{\lambda(A^TA)}$$ 证明：假设\(A^TA\)的特征值为\(\lambda_1,...,\lambda_n\geq 0\), 对应的单位特征向量为\(x^((1)),...,x^{(n)}\text{（该组特征向量为\)R^{n}\(空间的一组单位正交基有待证明）}\) 对\(\|x\|_2=1,x=k_1x^{(1)}+...+k_{n}x^{(n)}\)，有$$\sqrt{k2_1+...+k2_n}=1$$ 因此$$|Ax|2=\sqrt{xTATAx}=\sqrt{(k_1x{(1)}+...+k_{n}x)T(\lambda_1k_1x+...+\lambda_nkx{(n)})}=\sqrt{\lambda_1k2_1+...+\lambda_nk^2_n}\leq \sqrt{\lambda(A^TA)}$$ 当\(x\)取最大特征值对应的单位特征向量时，等号成立，因此\(\|A\|_2=\sqrt{\lambda_{max}(A^TA)}\)
- \(\infty\)范数 $$|A|{\infty}=\max{|x|{\infty}=1}|Ax|=\max_{1\leq i\leq n}|A^Te_i|1\text{（行绝对值和中的最大值）}$$ 证明：$$x=k_1e_1+...k_ne_n$$ 当\(\|x\|_{\infty}=1\)时，有\(\max_{1\leq i\leq n}|k_i|=1\)，因此 $$|Ax|=|k_1Ae_1+...+k_nAe_n|{\infty}\leq |Ae_1+...+Ae_n|=|A^Te_i|_1\text{（行绝对值和中的最大值）}$$ 当\(x^T=[1 ... 1]\)时，等号成立
5.3 条件数与相对误差
对线性系统\(Ax=b\)，假设其真实解为\(x_E\)，数值解为\(x_{NS}\)，则真实误差为\(TrueError = x_E-x_{NS}\)，一般真实解无法给出，因此，真实误差难以计算，但系统残差（residual）则可以通过下式给出：$$R=Ax_E-Ax_{NS}=b-Ax_{NS}$$
例：$$\begin{cases}1.02x_1+0.98x_2=2\0.98x_1+1.02x_2\end{cases}$$ 该系统的真实解为$$x_E=\begin{bmatrix}1\1\end{bmatrix}$$
- 若数值解为$$x_{NS}=\begin{bmatrix}1.02\1.02\end{bmatrix}$$ 真实误差为$$TrueError=-\begin{bmatrix}0.02\0.02\end{bmatrix}$$ 系统残差为 $$R=\begin{bmatrix}1.02 & 0.98\0.98 & 1.02\end{bmatrix}\begin{bmatrix}-0.02\-0.02\end{bmatrix}=\begin{bmatrix}-0.04\-0.04\end{bmatrix}$$ 系统残差与真实误差的数量级时一致的
- 若数值解为$$x_{NS}=\begin{bmatrix}2\0\end{bmatrix}$$ 真实误差为$$TrueError=\begin{bmatrix}-1\1\end{bmatrix}$$ 系统残差为 $$R=\begin{bmatrix}1.02 & 0.98\0.98 & 1.02\end{bmatrix}\begin{bmatrix}-1\1\end{bmatrix}=\begin{bmatrix}-0.04\0.04\end{bmatrix}$$ 此时，系统残差的数量级与真实误差的数量级差别很大，因此无法总是用残差来估计真实误差
一个好的想法是用范数来衡量估计误差，\(TrueError = x_E-x_{NS}\)， \(R=A(x_E-x_{NS})=A*TrueError\) ，注意 \(TrueError=A^{-1}R\)，根据范数的三角不等式，有 $$|TrueError|\leq |A^{-1}||R|$$ 进一步考虑相对误差和相对残差$$\frac{|TrueError|}{|x_E|},\ \frac{|R|}{|b|}$$ 有 $$\frac{|TrueError|}{|x_E|}\leq |A^{-1}|\frac{|b|}{|x_E|}\frac{|R|}{|b|}$$ 注意$$Ax_E=b,\ |b|\leq |A||x_E|\rightarrow |x_E|\geq \frac{|b|}{|A|}\rightarrow \frac{1}{|x_E|}\leq \frac{|A|}{|b|}$$ 因此有 $$\frac{|TrueError|}{|x_E|}\leq |A^{-1}||A|\frac{|R|}{|b|}$$ 另一方面，$$|R|\leq |A||TrueError|\rightarrow |TrueError|\geq \frac{|R|}{|A|}$$ 因此 $$\frac{|TrueError|}{|x_E|}\geq \frac{1}{|x_E|}\frac{|R|}{|A|}\rightarrow \frac{|TrueError|}{|x_E|}\geq \frac{1}{|A|}\frac{|b|}{|x_E|}\frac{|R|}{|b|}$$ 同时$$x_E=A^{-1}b\rightarrow |x_E|\leq |A^{-1}||b|\rightarrow \frac{|b|}{|x_E|}\geq \frac{1}{|A^{-1}|}$$ 因此有 $$\frac{|TrueError|}{|x_E|}\geq \frac{1}{|A|}\frac{|b|}{|x_E|}\frac{|R|}{|b|}\geq \frac{1}{|A^{-1}||A|}\frac{|R|}{|b|}$$ 综上，有 $$\frac{1}{|A^{-1}||A|}\frac{|R|}{|b|}\leq \frac{|TrueError|}{|x_E|}\leq |A^{-1}||A|\frac{|R|}{|b|}$$ \(\|A\|\|A^{-1}\|\)就称作条件数(Condition number)。
例：考虑第3节的线性方程组 $$9x_1-2x_2+3x_3+2x_4 = 54.5$$ $$2x_1+8x_2-2x_3+3x_4 = -14$$ $$-3x_1+2x_2+11x_3-4x_4 = 12.5$$ $$-2x_1+3x_2+2x_3+10x_4 = -21$$ 其真实解为 $$x_E=\begin{bmatrix}5 \ -2 \2.5\-1\end{bmatrix}$$ 而使用第3节的高斯赛德尔法迭代6次后，得到的数值解为 $$x_{NS}=\begin{bmatrix}4.98805\ -1.99511\ 2.49806 \ -1.00347 \end{bmatrix}$$ 因此真实误差为 $$TrueError=x_E-x_{NS}=\begin{bmatrix}0.0119\-0.0049\0.0019\0.0035\end{bmatrix}$$ 系统残差为 $$R=\begin{bmatrix}0.13009\ -0.00869\-0.03817\0.00001\end{bmatrix}$$ 现在以向量\(\infty\)为例，检查一下上述不等式的真实性：$$|x_E|{\infty}=5,\ |TrueError|=0.0119,\ |R|{\infty} = 0.13009,\ |b|=54.5, |A|=20,\ |A^{-1}|_{\infty}=0.19019$$ 注意 $$|A||A^{-1}|=3.8038\rightarrow \frac{1}{|A^{-1}||A|}\frac{|R|}{|b|}=0.00062749 \leq \frac{|TrueError|}{|x_E|}=0.00238 \leq |A^{-1}||A|\frac{|R|}{|b|}=0.009079$$
病态系统（ill-conditioned）：病态系统是指系统参数有很小的变化时都会引起解很大的变化的一类系统，且可以证明当条件数远大于\(1\)时，系统是一个病态的。
6. 总结
本节课主要介绍了求解线性方程组的两种迭代数值求解法，包括Jacobi法和高斯塞德尔法，其中Jacobi法是同步更新，而高斯赛德尔法是异步更新。对于特殊的三对角矩阵，可以采用Thomas算法求解，相比高斯消去、LU分解法更节省内存和时间。最后，介绍了向量范数和矩阵范数的概念，并引入条件数来确定线性系统数值解的相对误差范围。
