本
文
摘
要
关于线性方程组的数值解法一般有两类：直接法和迭代法
1. 直接法：直接法就是经过有限步算术运算，可求得线性方程组精确解的方法（若计算过程中没有舍入误差）。常用于求解低阶稠密矩阵方程组及某些大型稀疏矩阵方程组（如大型带状方程组）。
2. 迭代法：迭代法就是用某种极限过程去逐步逼近线性方程组精确解的方法。优点：存储单元少、程序设计简单、原始系数矩阵在计算过程始终不变等；缺点：存在收敛性及收敛速度问题。常用于求解大型稀疏矩阵方程组（尤其是由微分方程离散后得到的大型方程组）。本文主要讲直接方法。
1 高斯消去法
高斯消去法也称为逐次消去法，是一个古老的求解线性方程组的方法，特别是由其改进而得到的选主元素消去法、三角分解法仍然是目前计算机上常用的有效方法。
高斯消去法的过程写成定理即为：
设 Ax=b\mathbf{Ax}=\mathbf{b} ，其中A∈Rn×n\mathbf{A}\in \mathbb{R}^{n\times n} 。
1. 如果 akk(k)≠0(k=1,2,⋯,n)a_{kk}^{(k)}\ne 0(k=1,2,\cdots,n) ，则可通过高斯消去法将 Ax=b\mathbf{Ax}=\mathbf{b} 约化为等价的三角形线性方程组（图1），且计算公式为：图1: 三角形线性方程组（1）消元计算(k=1,2,⋯,n−1)(k=1,2,\cdots,n-1)
mik=aik(k)akk(k),i=k+1,⋯,n,m_{ik}=\frac{a_{ik}^{(k)}}{a_{kk}^{(k)}},\qquad\qquad\qquad i=k+1,\cdots,n,
aij(k+1)=aij(k)−mikakj(k),i,j=k+1,⋯,n,a_{ij}^{(k+1)}=a_{ij}^{(k)}-m_{ik}a_{kj}^{(k)},\quad i,j=k+1,\cdots,n,
bi(k+1)=bi(k)−mikbk(k),i=k+1,⋯,n.b_i^{(k+1)}=b_i^{(k)}-m_{ik}b_k^{(k)},\quad i=k+1,\cdots,n.
（2）回代计算
xn=bn(n)ann(n),x_n=\frac{b_n^{(n)}}{a_{nn}^{(n)}},
xi=bi(i)−∑j=i+1naij(i)xjaii(i),i=n−1,⋯,2,1x_i=\frac{b_i^{(i)}-\sum\limits_{j=i+1}^n a_{ij}^{(i)}x_j}{a_{ii}^{(i)}},\quad i=n-1,\cdots,2,1
2. 如果 A\mathbf{A} 为非奇异矩阵，则可通过高斯消去法（及交换两行的初等变换）将方程组 Ax=b\mathbf{Ax}=\mathbf{b} 约化为方程组（图1）。上述消元和回代过程总的乘除法次数为 n33+n2−n3≈n33\frac{n^3}{3}+n^2-\frac n3\approx \frac{n^3}{3} ，加减法次数为 n33+n22−56n≈n33\frac{n^3}{3}+\frac{n^2}{2}-\frac56n\approx \frac{n^3}{3} 。
高斯消去法写成MATLAB代码为：
function[x,t] =GaussEli(A,b)% 高斯消去法 % A: 系数矩阵 % b: 载荷矩阵 % x: 解矩阵 % t: 时间 tic [row,col] = size(A); x = zeros(row,1); % 消元过程 for k = 1:row-1 if A(k,k) == 0 print(主对角元不能为0！); else for i = k+1:row m_ik = A(i,k)/A(k,k); for j = k+1:col A(i,j) = A(i,j) - m_ik*A(k,j); end b(i) = b(i) - m_ik*b(k); end end end % 回代过程 x(end) = b(end)/A(end,end); for k = row-1:-1:1 sum = 0; for j = k+1:col sum = sum + A(k,j)*x(j); end x(k) = (b(k)-sum)/A(k,k); end t = toc; end但是，由算法可知，高斯消去法对于某些简单的矩阵可能会失败，例如
图22 列主元消去法
由高斯消去法可知，在消元过程中可能出现 akk(k)=0a_{kk}^{(k)}=0 的情况，这时消去法将无法进行；而且即使满足主元素 akk(k)≠0a_{kk}^{(k)}\ne 0 但很小时，用其作除数，会导致其他元素数量级的严重增长和舍入误差的扩散，使得最后得到的计算结果不可靠。这也就是小主元带来的麻烦，所以还应该避免采用绝对值小的主元素 akk(k)a_{kk}^{(k)} 。
由此来看，对于一般矩阵来说，最好每一步选取系数矩阵（或消元后的低阶矩阵）中绝对值最大的元素作为主元素，以使得高斯消去法具有较好的数值稳定性。这就是全主元素消去法的思想，不足之处在于其在选主元时要花费较多机器时间，目前主要使用的是列主元消去法。
列主元消去法总体过程与高斯消去法基本一致，唯一不同之处在于在选主元之前增加了一步按列选主元，其过程直接写成MATLAB代码为：
function[x,t] =MainColGaussEli(A,b)% 列主元素高斯消去法 % A: 系数矩阵 % b: 载荷矩阵 % x: 解矩阵 % t: 时间 tic [row,col] = size(A); x = zeros(row,1); % 消元过程 for k = 1:row-1 % 按列选主元 max=abs(A(k,k)); %先假设对角元的元素均为最大值； r=k; %记录此时的行号； for i=k+1:row %比较同一列其他元素值与主对角元元素的大小；i是行号； if max<abs(A(i,k)) max=abs(A(i,k)); r=i; %更新同一列最大值元素所在的行号； end end % 换行 A([k,r],:) = A([r,k],:); b([k,r],:) = b([r,k],:); if A(k,k) == 0 print(主对角元不能为0！); else for i = k+1:row m_ik = A(i,k)/A(k,k); for j = k+1:col A(i,j) = A(i,j) - m_ik*A(k,j); end b(i) = b(i) - m_ik*b(k); end end end % 回代过程 x(end) = b(end)/A(end,end); for k = row-1:-1:1 sum = 0; for j = k+1:col sum = sum + A(k,j)*x(j); end x(k) = (b(k)-sum)/A(k,k); end t = toc; end3 矩阵三角分解法
高斯消去法有很多变形，有的是高斯消去法的改进、改写，有的是用于某一类特殊性质矩阵的高斯消去法的简化。
这里矩阵的三角分解一般用的是矩阵的LU分解，即将矩阵分解为一个下三角矩阵L与一个上三角矩阵U的乘积，当然前提条件是矩阵的顺序主子式不为 00 ，下面是矩阵的LU分解的MATLAB代码：
function[L,U,t] =LUDecomp(A)% 矩阵的 LU 分解 % A: 系数矩阵 % L: L 矩阵 % U: U 矩阵 % t: 时间 tic [row,col] = size(A); L = zeros(row,col); L(logical(eye(size(L)))) = 1; U = zeros(row,col); % 消元过程 for k = 1:row-1 if A(k,k) == 0 print(主对角元不能为0！); else for i = k+1:row m_ik = A(i,k)/A(k,k); L(i,k) = m_ik; for j = k+1:col A(i,j) = A(i,j) - m_ik*A(k,j); end end end end for i = 1:row for j = 1:col if i <= j U(i,j) = A(i,j); end end end t = toc; end输入矩阵
图3可得其LU分解得到的矩阵为：
L = 1 0 0 0 1 0 2 -1 1 U = 1 1 1 0 4 -1 0 0 -2一旦实现了矩阵的 A\mathbf{A} 的LU分解，那么求解 Ax=b\mathbf{Ax}=\mathbf{b} 的问题就等价于求解两个三角形方程组：
（1） Ly=b\mathbf{Ly}=\mathbf{b} ，求 y\mathbf{y} ；
（2） Ux=y\mathbf{Ux}=\mathbf{y} ，求 x\mathbf{x} 。
这里只展示不选主元的三角分解法。其计算公式为：
（1） u1i=a1i(i=1,2,⋯,n),li1=ai1u11,i=2,3,⋯,nu_{1i}=a_{1i}(i=1,2,\cdots,n),l_{i1}=\frac{a_{i1}}{u_{11}},\quad i=2,3,\cdots,n 。
计算 U\mathbf{U} 的第 rr 行， L\mathbf{L} 的第 rr 列元素 (r=2,,3,⋯,n)(r=2,,3,\cdots,n) ：
（2） uri=ari−∑k=1r−1lrkuki,i=r,r+1,⋯,nu_{ri}=a_{ri}-\sum\limits_{k=1}^{r-1}l_{rk}u_{ki},\quad i=r,r+1,\cdots,n ；
（3） 且；lir=air−∑k=1r−1likukrurr,i=r+1,⋯,n,且r≠n；l_{ir}=\frac{a_{ir}-\sum\limits_{k=1}^{r-1}l_{ik}u_{kr}}{u_{rr}},\quad i=r+1,\cdots,n,\text{且}r\ne n；
求解 Ly=b,Ux=y\mathbf{Ly}=\mathbf{b},\mathbf{Ux}=\mathbf{y} 的计算公式：
（4） y1=b1y_1=b_1 ，
yi=bi−∑k=1i−1likyk,i=2,3,⋯,ny_i=b_i-\sum\limits_{k=1}^{i-1}l_{ik}y_k,\quad i=2,3,\cdots,n ；
（5） xn=ynunnx_n=\frac{y_n}{u_{nn}} ，
。xi=yi−∑k=i+1nuikxkuii,i=n−1,n−2,⋯,1。x_i=\frac{y_i-\sum\limits_{k=i+1}^n u_{ik}x_k}{u_{ii}},\quad i=n-1,n-2,\cdots,1。
写成MATLAB代码则为：
function[x,t] =LUMethod(A,b)% 矩阵三角分解法 (这里看不到 LU 矩阵是因为其元素存放在原矩阵 A 中) % A: 系数矩阵 % b: 载荷矩阵 % x: 解矩阵 % t: 时间 tic [row,~] = size(A); x = zeros(row,1); y = zeros(row,1); for i = 2:row A(i,1) = A(i,1)/A(1,1); end % 计算 U 的第 r 行，L 的第 r 列元素 % 杜利特尔(Doolittle)分解 for r = 2:row for i = r:row sum = 0; for k = 1:r-1 sum = sum + A(r,k)*A(k,i); end A(r,i) = A(r,i) - sum; end for i = r+1:row if r < row sum = 0; for k = 1:r-1 sum = sum + A(i,k)*A(k,r); end A(i,r) = (A(i,r) - sum)/A(r,r); end end end % 求解 Ly=b y(1) = b(1); for i = 2:row sum = 0; for k = 1:i-1 sum = sum + A(i,k)*y(k); end y(i) = b(i) - sum; end % 求解 Ux=y x(end) = y(end)/A(end,end); for i = row-1:-1:1 sum = 0; for k = i+1:row sum = sum + A(i,k)*x(k); end x(i) = (y(i) - sum)/A(i,i); end t = toc; end除了以上几种方法，还有针对特定矩阵的平方根法（对称正定矩阵）、追赶法（三对角矩阵），其方法也都比较简单，呃……还是来讲一下吧。
4 平方根法
应用有限元法求解结构力学问题时，最后归结为求解线性方程组，系数矩阵大多具有对称正定性质。所谓平方根法，就是利用对称正定矩阵的三角分解而得到的求解对称正定方程组的一种有效方法，目前在计算机上广泛应用平方根法解此类方程组。
先给出两个三角分解定理
对称阵的三角分解定理 设 A\mathbf{A} 为 nn 阶对称矩阵，且 A\mathbf{A} 的所有顺序主子式均不为零，则 A\mathbf{A} 可唯一分解为 A=LDLT\mathbf{A}=\mathbf{LDL}^\mathrm{T} ，其中 L\mathbf{L} 为单位下三角矩阵，D\mathbf{D} 为对角矩阵。
对称正定矩阵的三角分解或楚列斯基(Cholesky)分解 如果 A\mathbf{A} 为 nn 阶对称正定矩阵，则存在一个实的非奇异下三角矩阵 L\mathbf{L} 使 A=LLT\mathbf{A}=\mathbf{LL}^\mathrm{T} ，当限定 L\mathbf{L} 的对角元素为正时，这种分解是唯一的。从而可得平方根法得计算过程为：
对于 j=1,2,⋯,nj=1,2,\cdots,n
（1） ljj=(ajj−∑k=1j−1ljk2)12l_{jj} = (a_{jj}-\sum\limits_{k=1}^{j-1}l_{jk}^2)^{\frac12} ；
（2）。lij=aij−∑k=1j−1likljkljj,i=j+1,⋯,n。l_{ij}=\frac{a_{ij}-\sum\limits_{k=1}^{j-1}l_{ik}l_{jk}}{l_{jj}},\quad i=j+1,\cdots,n。
求解 Ax=b\mathbf{Ax}=\mathbf{b} ，即求解两个三角形方程组：
，求，求1∘Ly=b，求y;2∘LTx=y，求x1^{\circ}\quad \mathbf{Ly}=\mathbf{b}，求 \mathbf{y};\qquad\qquad 2^{\circ}\quad\mathbf{L}^{\mathrm{T}}\mathbf{x}=\mathbf{y}，求\mathbf{x} 。
（3）yi=bi−∑k=1i−1likyklii,i=1,2,⋯,n.y_i=\frac{b_i-\sum\limits_{k=1}^{i-1}l_{ik}y_k}{l_{ii}},\quad i = 1,2,\cdots,n.
（4）xi=yi−∑k=i+1nlkixklii,i=n,n−1,⋯,1.x_i = \frac{y_i-\sum\limits_{k=i+1}^n l_{ki}x_k}{l_{ii}},\quad i=n,n-1,\cdots,1. [1]
其MATLAB代码为：
function[x,t] =SquareRoot(A,b)% 平方根法（L 矩阵元素存在原 A 矩阵中） % A: 系数矩阵 % b: 载荷矩阵 % x: 解矩阵 % t: 时间 tic [n,~] = size(A); x = zeros(n,1); y = zeros(n,1); for j = 1:n sum = 0; for k = 1:j-1 sum = sum + A(j,k)^2; end A(j,j) = (A(j,j) - sum)^0.5; for i = j+1:n sum = 0; for k = 1:j-1 sum = sum + A(i,k)*A(j,k); end A(i,j) = (A(i,j) - sum)/A(j,j); end end for i = 1:n sum = 0; for k = 1:i-1 sum = sum + A(i,k)*y(k); end y(i) = (b(i) - sum)/A(i,i); end for i = n:-1:1 sum = 0; for k = i+1:n sum = sum + A(k,i)*x(k); end x(i) = (y(i) - sum)/A(i,i); end t = toc; end从平方根法计算过程可以看出，在计算 L\mathbf{L} 的元素 liil_{ii} 时需要计用到开方运算，为避免开方，我们提出改进的平方根法：
其计算过程为：
d1=a11.d_1=a_{11}.
对于 i=2,3,\cdots,n.
（1） t_{ij}=a_{ij}-\sum\limits_{k=1}^{j-1}t_{ik}l_{jk},\quad j = 1,2,\cdots,i-1 ；
（2） l_{ij}=\frac{t_{ij}}{d_j},\quad j = 1,2,\cdots,i-1 ；
（3） d_i=a_{ii}-\sum\limits_{k=1}^{i-1}t_{ik}l_{ik} 。
求解 \mathbf{Ly}=\mathbf{b}，\mathbf{DL}^{\mathrm{T}}\mathbf{x}=\mathbf{y} 计算公式：
（4） y_1=b_1 ；
y_i=b_i-\sum\limits_{k=1}^{i-1}l_{ik}y_k,\quad i=2,3,\cdots,n 。
（5） x_n=\frac{y_n}{d_n} ；
x_i=\frac{y_i}{d_i}-\sum\limits_{k=i+1}^n l_{ki}x_k,\quad i = n-1,\cdots,2,1 。
写成MATLAB代码为：
function[x,t] =ImprovedSquareRoot(A,b)% 改进的平方根法（L 矩阵元素存在原 A 矩阵中） % A: 系数矩阵 % b: 载荷矩阵 % x: 解矩阵 % t: 时间 tic [n,m] = size(A); T = zeros(n,m); x = zeros(n,1); y = zeros(n,1); for i = 2:n for j = 1:i-1 sum = 0; for k = 1:j-1 sum = sum + T(i,k)*A(j,k); end T(i,j) = (A(i,j) - sum); end for j = 1:i-1 A(i,j) = T(i,j)/A(j,j); end sum = 0; for k = 1:i-1 sum = sum + T(i,k)*A(i,k); end A(i,i) = A(i,i) - sum; end y(1) = b(1); for i = 2:n sum = 0; for k = 1:i-1 sum = sum + A(i,k)*y(k); end y(i) = b(i) - sum; end x(n) = y(n)/A(n,n); for i = n-1:-1:1 sum = 0; for k = i+1:n sum = sum + A(k,i)*x(k); end x(i) = y(i)/A(i,i) - sum; end t = toc; end5 追赶法
在一些实际问题中，例如解常微分方程边值问题，解热传导方程以及船体数学放样中建立三次样条函数等，都会要求解系数矩阵为对角占优的三对角线方程组（图4）。
图4图4可简记为 \mathbf{Ax}=\mathbf{f} 。
我们利用矩阵的三角分解，得到 \mathbf{A}=\mathbf{LU} ，其中 \mathbf{L} 为下三角矩阵， \mathbf{U} 为单位上三角矩阵，即
图5求解 \mathbf{Ax}=\mathbf{f} 等价于解两个方程组：
1^{\circ}\quad \mathbf{Ly}=\mathbf{f}，求\mathbf{y}；\qquad\qquad 2^{\circ}\quad\mathbf{Ux}=\mathbf{y}，求\mathbf{x} 。
从而得到求解三对角线方程组的追赶法公式：
（1）计算 \{\beta_i\} 的递推公式
\beta_1=\frac{c_1}{b_1} ，
\beta_i=\frac{c_i}{b_i-a_i\beta_{i-1}},\quad i = 2,3,\cdots,n-1 ；
（2）解 \mathbf{Ly}=\mathbf{f}
y_1=\frac{f_1}{b_1} ，
y_i=\frac{f_i-a_iy_{i-1}}{b_i-a_i\beta_{i-1}},\quad i = 2,3,\cdots,n ；
（3）解 \mathbf{Ux}=\mathbf{y}
x_n=y_n ，
x_i=y_i-\beta_ix_{i+1},\quad i = n-1,n-2,\cdots,2,1 。
我们将计算系数 \beta_1\rightarrow\beta_2\rightarrow\cdots\rightarrow\beta_{n-1} 及 y_1\rightarrow y_2\rightarrow\cdots\rightarrow y_n 的过程称为追的过程，将计算方程组的解 x_n\rightarrow x_{n-1}\rightarrow\cdots\rightarrow x_1 的过程称为赶的过程。
其MATLAB代码为：
function[x,t] =ZhuiGan(a,b,c,f)% 追赶法（三对角矩阵） % a: 系数矩阵 -1 对角线元素 % b: 系数矩阵对角线元素 % c: 系数矩阵 1 对角线元素 % f: 载荷矩阵 % x: 解矩阵 % t: 时间 tic n = length(b); a = [0; a]; beta = zeros(n-1,1); y = zeros(n,1); x = zeros(n,1); beta(1) = c(1)/b(1); for i = 2:n-1 beta(i) = c(i)/(b(i)-a(i)*beta(i-1)); end y(1) = f(1)/b(1); for i = 2:n y(i) = (f(i)-a(i)*y(i-1))/(b(i)-a(i)*beta(i-1)); end x(n) = y(n); for i = n-1:-1:1 x(i) = y(i) - beta(i)*x(i+1); end t = toc; end6 总结
下面把几种方法放在一起，来看看这几种方法的运行效率，这里给的系数矩阵是一个 100 阶的三对角矩阵。输入下述代码：
clear n = 500; T1 = zeros(n,1); T2 = zeros(n,1); T3 = zeros(n,1); T4 = zeros(n,1); T5 = zeros(n,1); T6 = zeros(n,1); for k = 1:n a = -1*ones(99,1); b = 2*ones(100,1); c = -1*ones(99,1); A = diag(b) + diag(a,-1) + diag(c,1); f = zeros(99,1); f = [1; f]; [x1,t1] = GaussEli(A,f); [x2,t2] = MainColGaussEli(A,f); [x3,t3] = LUMethod(A,f); [x4,t4] = SquareRoot(A,f); [x5,t5] = ImprSquareRoot(A,f); [x6,t6] = ZhuiGan(a,b,c,f); T1(k) = t1; T2(k) = t2; T3(k) = t3; T4(k) = t4; T5(k) = t5; T6(k) = t6; end x = 1:n; plot(x,T1,r-,x,T2,g-,x,T3,b-,x,T4,c-,x,T5,m-,x,T6,k-) legend(GaussEli,MainColGaussEli,LUMethod,... SquareRoot,ImprSquareRoot,ZhuiGan) axis([0 n 0 0.002]) % 输出 fprintf( Method | Solve:x | Time \n) fprintf(================|========================|========\n) fprintf( GaussEli |); fprintf(%8.4f,x1); fprintf(|%8.4f,t1); fprintf(\n) fprintf(----------------|------------------------|--------\n) fprintf(MainColGaussEli |); fprintf(%8.4f,x2); fprintf(|%8.4f,t2); fprintf(\n) fprintf(----------------|------------------------|--------\n) fprintf( LUMethod |); fprintf(%8.4f,x3); fprintf(|%8.4f,t3); fprintf(\n) fprintf(----------------|------------------------|--------\n) fprintf( SquareRoot |); fprintf(%8.4f,x4); fprintf(|%8.4f,t4); fprintf(\n) fprintf(----------------|------------------------|--------\n) fprintf( ImprSquareRoot |); fprintf(%8.4f,x5); fprintf(|%8.4f,t5); fprintf(\n) fprintf(----------------|------------------------|--------\n) fprintf( ZhuiGan |); fprintf(%8.4f,x6); fprintf(|%8.4f,t6); fprintf(\n)运行得到数值结果为：
Method | Solve:x | Time ================|========================|======== GaussEli | 0.9901 ........ 0.0099| 0.0008 ----------------|------------------------|-------- MainColGaussEli | 0.9901 ........ 0.0099| 0.0014 ----------------|------------------------|-------- LUMethod | 0.9901 ........ 0.0099| 0.0010 ----------------|------------------------|-------- SquareRoot | 0.9901 ........ 0.0099| 0.0007 ----------------|------------------------|-------- ImprSquareRoot | 0.9901 ........ 0.0099| 0.0006 ----------------|------------------------|-------- ZhuiGan | 0.9901 ........ 0.0099| 0.0000 >>图像为：
图6: 运行时间对比从数值结果或图形上都可以很清楚地看到各种方法的效率高低，很直观地展现了追赶法对于解三对角矩阵的优势。
