求解线性方程组

0. 引言
1. 直接消元法
- 1.1 高斯消元法
- - 1.1.1 顺序高斯消去法
  - 1.1.2 主元素高斯消去法
  - 1.1.3 高斯-约当(Gauss-Jordan)消去法
- 1.2 三角分解方法
- 1.3 解对角方程的追赶法
- 1.4 误差分析
2. 迭代法
- 2.1 雅克比(Jacobi)迭代法
- 2.2 高斯-赛德尔迭代法
- 2.3 超松弛迭代法
- 2.4 迭代方法的收敛性
小结
参考文献

0. 引言

在很多工程求解问题中，大部分的数值计算都可以归结为求解线性方程组的问题，所以求解线性方程组是数值计算中最为重要的一部分。其中还有很多其他数值计算中，其中的关键步骤都涉及到线性方程组的求解。由克莱姆法则可知，线性方程组Ax=bAx=bAx=b中，系数矩阵的行列式det⁡(A)≠0\det(A)\neq 0det(A)=0，则线性方程组具有唯一解。克莱姆法则在计算高维度情况下的矩阵运算中计算量很大，不适合求解一般情况下的线性方程组。常见的求解线性方程组的方法有两种，即直接方法和迭代方法。直接方法是指通过有限的步骤计算直接求解出方程组的精确解的方法。迭代的方法值得是用极限的过程去逐步逼近线性方程组精确解的方法。本节将实现这几种求解的方法。

1. 直接消元法

直接消元的方法是一种便捷的求解线性方程组的方法。高斯消去方法是最为古老的一种线性方程组求解的方法。通过改进高斯消去法得到更加便捷的主元素消去法、矩阵的三角分解(即Doolittle分解方法)，是用于求解低阶稠密方程组的有效方法。特点是，通过消元将一般的线性方程组求解问题转化为上三角方程组求解问题。

1.1 高斯消元法

1.1.1 顺序高斯消去法

一般地，设nnn阶线性方程组为
{a11(0)x1+a12(0)x2+...+a1n(0)xn=b1(0)a21(0)x1+a22(0)x2+...+a2n(0)xn=b2(0):an1(0)x1+an2(0)x2+...+ann(0)xn=bn(0)\begin{cases} a_{11}^{(0)}x_{1}+a_{12}^{(0)}x_{2}+...+a_{1n}^{(0)}x_{n}=b_{1}^{(0)}\\ a_{21}^{(0)}x_{1}+a_{22}^{(0)}x_{2}+...+a_{2n}^{(0)}x_{n}=b_{2}^{(0)}\\ :\\ a_{n1}^{(0)}x_{1}+a_{n2}^{(0)}x_{2}+...+a_{nn}^{(0)}x_{n}=b_{n}^{(0)}\\ \end{cases}⎩⎪⎪⎪⎨⎪⎪⎪⎧a11(0)x1+a12(0)x2+...+a1n(0)xn=b1(0)a21(0)x1+a22(0)x2+...+a2n(0)xn=b2(0):an1(0)x1+an2(0)x2+...+ann(0)xn=bn(0)
写成增广矩阵的方式为
(A(0)∣b(0))=[a11(0)a12(0)...a1n(0)b1(0)a21(0)a22(0)...a2n(0)b2(0):::::an1(0)an2(0)...ann(0)bn(0)](A^{(0)}|b^{(0)})=\left[\begin{array}{cccc} a_{11}^{(0)}&a_{12}^{(0)}&...&a_{1n}^{(0)}&b_{1}^{(0)}\\ a_{21}^{(0)}&a_{22}^{(0)}&...&a_{2n}^{(0)}&b_{2}^{(0)}\\ :&:&:&:&:\\ a_{n1}^{(0)}&a_{n2}^{(0)}&...&a_{nn}^{(0)}&b_{n}^{(0)}\\ \end{array}\right](A(0)∣b(0))=⎣⎢⎢⎢⎡a11(0)a21(0):an1(0)a12(0)a22(0):an2(0)......:...a1n(0)a2n(0):ann(0)b1(0)b2(0):bn(0)⎦⎥⎥⎥⎤
设li1=ai1(0)/a11(0),(i=2,3,...,n)l_{i1}=a_{i1}^{(0)}/a_{11}^{(0)},(i=2,3,...,n)li1=ai1(0)/a11(0),(i=2,3,...,n)，将第一行的−li1-l_{i1}−li1倍加到第iii行上，那么矩阵变为
(A(1)∣b(1))[a11(0)a12(0)...a1n(0)b1(0)0a22(1)...a2n(1)b2(1):::::an1(1)an2(1)...ann(1)bn(1)](A^{(1)}|b^{(1)})\left[\begin{array}{cccc} a_{11}^{(0)}&a_{12}^{(0)}&...&a_{1n}^{(0)}&b_{1}^{(0)}\\ 0&a_{22}^{(1)}&...&a_{2n}^{(1)}&b_{2}^{(1)}\\ :&:&:&:&:\\ a_{n1}^{(1)}&a_{n2}^{(1)}&...&a_{nn}^{(1)}&b_{n}^{(1)}\\ \end{array}\right](A(1)∣b(1))⎣⎢⎢⎢⎡a11(0)0:an1(1)a12(0)a22(1):an2(1)......:...a1n(0)a2n(1):ann(1)b1(0)b2(1):bn(1)⎦⎥⎥⎥⎤
其中，式子中
aij(1)=aij(0)−li1a1j(0),(i=2,3,...,n;j=1,2,...,n)a_{ij}^{(1)}=a_{ij}^{(0)}-l_{i1}a_{1j}^{(0)},(i=2,3,...,n;j=1,2,...,n)aij(1)=aij(0)−li1a1j(0),(i=2,3,...,n;j=1,2,...,n)

bi(1)=bi(0)−li1b1(0),(i=2,3,...,n)b_{i}^{(1)}=b_{i}^{(0)}-l_{i1}b_{1}^{(0)},(i=2,3,...,n)bi(1)=bi(0)−li1b1(0),(i=2,3,...,n)

第kkk步骤的时候，设主元素akk(k−1)≠0a_{kk}^{(k-1)}\neq 0akk(k−1)=0时，设lik=aik(k−1)/akk(k−1)l_{ik}=a_{ik}^{(k-1)}/a_{kk}^{(k-1)}lik=aik(k−1)/akk(k−1)，并且将(A(k−1)∣b(k−1))(A^{(k-1)|b^{(k-1)}})(A(k−1)∣b(k−1))的第kkk行的−lik-l_{ik}−lik倍加到第iii行(i=k+1,k+2,...,n)(i=k+1,k+2,...,n)(i=k+1,k+2,...,n)得到：
(A(k)∣b(k))=[a11(0)a12(0)...a1k(0)a1(k+1)(0)...a1n(0)b1(0)::::::akk(k−1)ak(k+1)(0)...akn(k−1)bk(k−1)a(k+1)(k+1)(k)...a(k+1)n(k)bn(k)::::an(k+1)(k)...ann(k)bn(k)](A^{(k)}|b^{(k)})=\left[\begin{array}{cccc} a_{11}^{(0)}&a_{12}^{(0)}&...&a_{1k}^{(0)}&a_{1(k+1)}^{(0)}&...&a_{1n}^{(0)}&b_{1}^{(0)}\\ &&:&:&:&:&:&:\\ &&&a_{kk}^{(k-1)}&a_{k(k+1)}^{(0)}&...&a_{kn}^{(k-1)}&b_{k}^{(k-1)}\\ &&&&a_{(k+1)(k+1)}^{(k)}&...&a_{(k+1)n}^{(k)}&b_{n}^{(k)}\\ &&&&:&:&:&:\\ &&&&a_{n(k+1)}^{(k)}&...&a_{nn}^{(k)}&b_{n}^{(k)}\\ \end{array}\right](A(k)∣b(k))=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡a11(0)a12(0)...:a1k(0):akk(k−1)a1(k+1)(0):ak(k+1)(0)a(k+1)(k+1)(k):an(k+1)(k)...:......:...a1n(0):akn(k−1)a(k+1)n(k):ann(k)b1(0):bk(k−1)bn(k):bn(k)⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤

式中
aij(k)=aij(k−1)−likakj(k−1),(i,j=k+1,k+2,...,n)a_{ij}^{(k)}=a_{ij}^{(k-1)}-l_{ik}a_{kj}^{(k-1)},(i,j=k+1,k+2,...,n)aij(k)=aij(k−1)−likakj(k−1),(i,j=k+1,k+2,...,n)

bi(k)=bi(k−1)−likbk(k−1),(i=k+1,k+2,...,n)b_{i}^{(k)}=b_{i}^{(k-1)}-l_{ik}b_{k}^{(k-1)},(i=k+1,k+2,...,n)bi(k)=bi(k−1)−likbk(k−1),(i=k+1,k+2,...,n)

所以经过n−1n-1n−1步消元之后，可以得到原方程组的同解三角形方程组A(n−1)x=b(n−1)A^{(n-1)}x=b^{(n-1)}A(n−1)x=b(n−1)，即
[a11(0)a12(0)...a1n(0)a22(1)...a2n(1):...:ann(n−1)][x1x2:xn]=[b1(n−1)b2(n−1):bn(n−1)]\left[\begin{array}{cccc} a_{11}^{(0)}&a_{12}^{(0)}&...&a_{1n}^{(0)}\\ &a_{22}^{(1)}&...&a_{2n}^{(1)}\\ &:&...&:\\ &&&a_{nn}^{(n-1)}\\ \end{array}\right] \left[\begin{array}{cccc} x_{1}\\ x_{2}\\ :\\ x_{n}\\ \end{array}\right]= \left[\begin{array}{cccc} b_{1}^{(n-1)}\\ b_{2}^{(n-1)}\\ :\\ b_{n}^{(n-1)}\\ \end{array}\right]⎣⎢⎢⎢⎡a11(0)a12(0)a22(1):.........a1n(0)a2n(1):ann(n−1)⎦⎥⎥⎥⎤⎣⎢⎢⎡x1x2:xn⎦⎥⎥⎤=⎣⎢⎢⎢⎡b1(n−1)b2(n−1):bn(n−1)⎦⎥⎥⎥⎤

综上所述，整个方程的求解过程如下所示
对k=1,2,...,n−1k=1,2,...,n-1k=1,2,...,n−1逐次计算
{lik=aik(k−1)/akk(k−1),(i=k+1,k+2,...,n)aij(k)=aij(k−1)−likakj(k−1),(i,j=k+1,k+2,...,n)bi(k)=bi(k−1)−likbi(k−1),(i=k+1,...,n)\begin{cases} l_{ik}=a_{ik}^{(k-1)}/a_{kk}^{(k-1)}&,(i=k+1,k+2,...,n)\\ a_{ij}^{(k)}=a_{ij}^{(k-1)}-l_{ik}a_{kj}^{(k-1)}&,(i,j=k+1,k+2,...,n)\\ b_{i}^{(k)}=b_{i}^{(k-1)}-l_{ik}b_{i}^{(k-1)}&,(i=k+1,...,n) \end{cases}⎩⎪⎨⎪⎧lik=aik(k−1)/akk(k−1)aij(k)=aij(k−1)−likakj(k−1)bi(k)=bi(k−1)−likbi(k−1),(i=k+1,k+2,...,n),(i,j=k+1,k+2,...,n),(i=k+1,...,n)

逐步回代求得原方程组的解，回代过程如下所示：
{xn=bn(n−1)/ann(n−1)xk=(bk(k−1)−∑j=k+1nakj(k−1)xj)/akk(k−1),(k=n−1,n−2,...,1)\begin{cases} x_{n}=b_{n}^{(n-1)}/a_{nn}^{(n-1)}\\ x_{k}=(b_{k}^{(k-1)}-\sum\limits_{j=k+1}^{n}a_{kj}^{(k-1)}x_{j})/a_{kk}^{(k-1)}&,(k=n-1,n-2,...,1)\\ \end{cases}⎩⎨⎧xn=bn(n−1)/ann(n−1)xk=(bk(k−1)−j=k+1∑nakj(k−1)xj)/akk(k−1),(k=n−1,n−2,...,1)

一般地，如果nnn阶矩阵AAA的顺序主子式均不为零，那么高斯消去法中求解线性方程组Ax=bAx=bAx=b的主元素满足akk(k−1)≠0a_{kk}^{(k-1)}\neq 0akk(k−1)=0
代码实现如下所示：

def GaussianSeq(Mat,Bais):if (Mat.shape[0]!=Mat.shape[1]) or (Mat.shape[0]!=Bais.shape[0]):raise ValueError("size of matrix doesn't match!")A = Mat.copy()b = Bais.copy()n = A.shape[0]# 消元过程for k in range(n):for i in range(k+1,n):l_ik = A[i,k]/A[k,k]for j in range(k,n):A[i, j] = A[i, j] - l_ik * A[k, j]b[i] = b[i] - l_ik*b[k]# 回代过程x = np.zeros((n,))x[n-1] = b[n-1]/A[n-1,n-1]for k in range(n-2,-1,-1):x[k] = (b[k]-np.sum(A[k,k+1:]*x[k+1:]))/A[k,k]return x

1.1.2 主元素高斯消去法

主元素高斯消去法主要解决的是舍入误差的扩大和传播而提出来的。事实上，从消元过程中可以看出来，当∣lik∣=∣aik(k−1)/akk(k−1)∣>1|l_{ik}|=|{a_{ik}^{(k-1)}}/{a_{kk}^{(k-1)}}|>1∣lik∣=∣aik(k−1)/akk(k−1)∣>1时候，若akj(k∗−1)a_{kj}^{(k*-1)}akj(k∗−1)有舍入误差ek−1e_{k-1}ek−1，则经过计算aij(k)=aij(k−1)−likakj(k−1)a_{ij}^{(k)}=a_{ij}^{(k-1)}-l_{ik}a_{kj}^{(k-1)}aij(k)=aij(k−1)−likakj(k−1)，akja_{kj}akj会有比ek−1e_{k-1}ek−1更大的舍入误差ek=∣lik∣ek−1e_{k}=|l_{ik}|e_{k-1}ek=∣lik∣ek−1。解决的方法是用交换的方法或者交换未知数次序的方法，选择绝对值尽可能大的系数作为第kkk次消元过程中的主元素akk(k−1)a_{kk}^{(k-1)}akk(k−1)。所以基本的算法过程如下所示：
输入： 矩阵系数AAA和bbb

记(A∣b)(A|b)(A∣b)为An×(n+1)A_{n\times{(n+1)}}An×(n+1)

for k=1,2,...,n−1\text{ for }k=1,2,...,n-1 for k=1,2,...,n−1

消元过程：按照第kkk列选定主元素，确定rrr使得∣ark∣=max⁡k≤i≤k+1∣aik∣|a_{rk}|=\max\limits_{k\leq i \leq k+1}|a_{ik}|∣ark∣=k≤i≤k+1max∣aik∣

判断∣ark∣|a_{rk}|∣ark∣是否等于0，如果是，则系数矩阵为奇异矩阵，退出计算过程；如果不是，则
- 如果r>kr>kr>k，则
  - 对于j=k+1,...,n+1j=k+1,...,n+1j=k+1,...,n+1交换r,kr,kr,k两行元素：arj↔akja_{rj}\leftrightarrow a_{kj}arj↔akj
- 对于i=k+1,k+2,...i=k+1,k+2,...i=k+1,k+2,...，计算
  - aik=aik/akka_{ik}=a_{ik}/a_{kk}aik=aik/akk
    - 对于j=k+1,k+2,...,n,n+1j=k+1,k+2,...,n,n+1j=k+1,k+2,...,n,n+1，计算
      - aij=aij−aikakja_{ij}=a_{ij}-a_{ik}a_{kj}aij=aij−aikakj
        end for \text{ end for } end for

回代过程：ann+1=ann+1/anna_{nn+1}=a_{nn+1}/a_{nn}ann+1=ann+1/ann
for k=1,2,...,n−1\text{ for }k=1,2,...,n-1 for k=1,2,...,n−1

计算akn+1=(akn+1−∑j=k+1nakjajn+1)/akka_{kn+1}=(a_{kn+1}-\sum\limits_{j=k+1}^{n}a_{kj}a_{jn+1})/a_{kk}akn+1=(akn+1−j=k+1∑nakjajn+1)/akk

end for \text{ end for } end for
输出： 未知数方程的解x=(a1n+1,a2n+1,...,ann+1)x=(a_{1n+1},a_{2n+1},...,a_{nn+1})x=(a1n+1,a2n+1,...,ann+1)
用numpy实现过程为：

def GaussianMain(Mat,Bais):if (Mat.shape[0]!=Mat.shape[1]) or (Mat.shape[0]!=Bais.shape[0]):raise ValueError("size of matrix doesn't match!")A = Mat.copy()b = Bais.copy()n = A.shape[0]# 消元过程for k in range(n):# 选定主元素r = kfor j in range(k,n):if np.abs(A[j,k])>np.abs(A[r,k]):r = j# 交换行A[[r, k], :] = A[[k, r], :]b[r],b[k] = b[k],b[r]if np.abs(A[k,k])< 1e-150:raise np.linalg.LinAlgErrorfor i in range(k+1,n):l_ik = A[i,k]/A[k,k]for j in range(k,n):A[i, j] = A[i, j] - l_ik * A[k, j]b[i] = b[i] - l_ik*b[k]# 回代过程x = np.zeros((n,))x[n-1] = b[n-1]/A[n-1,n-1]for k in range(n-2,-1,-1):x[k] = (b[k]-np.sum(A[k,k+1:]*x[k+1:]))/A[k,k]return x

1.1.3 高斯-约当(Gauss-Jordan)消去法

高斯消去法是将系数矩阵AAA变为上三角矩阵。高斯约当消去法是将对角线上方和下方的元素都消去，使得系数矩阵变为单位矩阵，这样就不需要回代过程。
所以计算过程变为：
{akj(k)=akj(k−1)/akk(k−1),(j=k,k+1,...,n)aij(k)=aij(k−1)−aik(k−1)×akj(k),(i=1,2,...,n,i≠k,j=k,k+1,...,n)bk(k)=bk(k−1)/akk(k−1)bk(k)=bk(k−1)−aik(k−1)×bk(k)\begin{cases} a_{kj}^{(k)}=a_{kj}^{(k-1)}/a_{kk}^{(k-1)}&,(j=k,k+1,...,n)\\ a_{ij}^{(k)}=a_{ij}^{(k-1)}-a_{ik}^{(k-1)}\times a_{kj}^{(k)}&,(i=1,2,...,n,i\neq k,j=k,k+1,...,n)\\ b_{k}^{(k)}=b_{k}^{(k-1)}/a_{kk}^{(k-1)}\\ b_{k}^{(k)}=b_{k}^{(k-1)}-a_{ik}^{(k-1)}\times b_{k}^{(k)}\\ \end{cases}⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧akj(k)=akj(k−1)/akk(k−1)aij(k)=aij(k−1)−aik(k−1)×akj(k)bk(k)=bk(k−1)/akk(k−1)bk(k)=bk(k−1)−aik(k−1)×bk(k),(j=k,k+1,...,n),(i=1,2,...,n,i=k,j=k,k+1,...,n)

用numpy实现为：

def Gauss_Jordan(Mat,Bais):if (Mat.shape[0] != Mat.shape[1]) or (Mat.shape[0] != Bais.shape[0]):raise ValueError("size of matrix doesn't match!")A = Mat.copy()b = Bais.copy()n = A.shape[0]# 消元过程for k in range(n):# 选定主元素r = kfor j in range(k, n):if np.abs(A[j, k]) > np.abs(A[r, k]):r = j# 交换行A[[r, k], :] = A[[k, r], :]b[r], b[k] = b[k], b[r]if np.abs(A[k, k]) < 1e-150:raise np.linalg.LinAlgErrortmp = A[k, k]for j in range(k,n):A[k,j] = A[k,j]/tmpb[k] = b[k]/tmpfor i in range(n):if i!=k:tmp = A[i,k]for j in range(k,n):A[i,j] = A[i,j] - tmp *A[k,j]b[i] = b[i] - tmp*b[k]return b

1.2 三角分解方法

这种算法也称作是Doolittle分解法，或者是LU分解。一般地，任意一个矩阵不一定有LU分解。有以下定理：
定理一：对任意矩阵A∈Rn×nA\in{R^{n\times n}}A∈Rn×n，若AAA的顺序主子式均不为零，则AAA具有唯一的LU分解。
定理二：若矩阵A∈Rn×nA\in{R^{n\times n}}A∈Rn×n非奇异，并且其LULULU分解存在，则AAA的LULULU分解是唯一的。
这样，如果对于线性方程组Ax=bAx=bAx=b中系数矩阵A=LUA=LUA=LU，那么方程可以写为(LU)x=b(LU)x=b(LU)x=b。所以设
{Ux=yLy=b\begin{cases} Ux=y\\ Ly=b\\ \end{cases}{Ux=yLy=b

所以进一步我们改进这种方程的求解方法。
设
A=[a11a12...a1na21a22...a2n::::an1an2...ann]=[l11l21l22::ln1ln2...lnn][u11u12...u1nu22...u2n::unn]A= \left [\begin{array}{cccc} a_{11}&a_{12}&...&a_{1n}\\ a_{21}&a_{22}&...&a_{2n}\\ :&:&:&:\\ a_{n1}&a_{n2}&...&a_{nn}\\ \end{array}\right] =\left [\begin{array}{cccc} l_{11}\\ l_{21}&l_{22}\\ :&:\\ l_{n1}&l_{n2}&...&l_{nn} \end{array}\right] \left [\begin{array}{cccc} u_{11}&u_{12}&...&u_{1n}\\ &u_{22}&...&u_{2n}\\ &&:&:\\ &&&u_{nn}\\ \end{array}\right] A=⎣⎢⎢⎡a11a21:an1a12a22:an2......:...a1na2n:ann⎦⎥⎥⎤=⎣⎢⎢⎡l11l21:ln1l22:ln2...lnn⎦⎥⎥⎤⎣⎢⎢⎡u11u12u22......:u1nu2n:unn⎦⎥⎥⎤
我们在这里设lii=1,(i=1,2,..,n)l_{ii}=1,(i=1,2,..,n)lii=1,(i=1,2,..,n)
这样可以得到，对于k=1,2,...,nk=1,2,...,nk=1,2,...,n
ukj=akj−∑r=1k−1lkrurj,(j=k,k+1,...,n)lik=(aik−∑r=1k−1lirurk)/ukk,(i=k+1,k+2,...,n)u_{kj}=a_{kj}-\sum_{r=1}^{k-1}l_{kr}u_{rj},(j=k,k+1,...,n)\\ l_{ik}=(a_{ik}-\sum_{r=1}^{k-1}l_{ir}u_{rk})/u_{kk},(i=k+1,k+2,...,n) ukj=akj−r=1∑k−1lkrurj,(j=k,k+1,...,n)lik=(aik−r=1∑k−1lirurk)/ukk,(i=k+1,k+2,...,n)
回代过程为：
{y1=b1yk=bk−∑j=1k−1lkjyi,(k=2,3,...,n)\begin{cases} y_{1}=b_{1}\\ y_{k}=b_{k}-\sum\limits_{j=1}^{k-1}l_{kj}y_{i}&,(k=2,3,...,n) \end{cases}⎩⎪⎨⎪⎧y1=b1yk=bk−j=1∑k−1lkjyi,(k=2,3,...,n)

{xn=yn/unnxk=(yk−∑j=k+1nukjxj)/ukk,(k=n−1,n−2,...,1)\begin{cases} x_{n}=y_{n}/u_{nn}\\ x_{k}=(y_{k}-\sum\limits_{j=k+1}^{n}u_{kj}x_{j})/u_{kk}&,(k=n-1,n-2,...,1) \end{cases}⎩⎨⎧xn=yn/unnxk=(yk−j=k+1∑nukjxj)/ukk,(k=n−1,n−2,...,1)
用numpy实现为：

def Doolittle_LU(A):if A.shape[0]!=A.shape[1]:raise ValueError("size of matrix doesn't match!")L_U_mat = np.zeros(A.shape)n = A.shape[0]for k in range(n):for j in range(k,n):if k==0:L_U_mat[k,j] = A[k,j]else:L_U_mat[k,j] = A[k,j] - np.sum(L_U_mat[k,:k]*L_U_mat[:k,j])for i in range(k+1,n):if k==0:L_U_mat[i,k] = A[i,k]/L_U_mat[k,k]else:L_U_mat[i,k] = (A[i,k] -np.sum(L_U_mat[i,:k]*L_U_mat[:k,k]))/L_U_mat[k,k]#L = np.tril(L_U_mat,k=-1)#U = np.triu(L_U_mat,k=0)#L = L + np.diag([1.0]*n)return L_U_mat
def Doolittle_solve(A,b):if A.shape[0]!=b.shape[0]:raise ValueError("size of matrix doesn't match!")L_U_mat = Doolittle_LU(A)n = L_U_mat.shape[0]y = np.zeros((n,))for k in range(n):if k==0:y[k] = b[k]else:y[k] = b[k] - np.sum(L_U_mat[k,:k]*y[:k])x = np.zeros((n,))for k in range(n-1,-1,-1):if k == n-1:x[k] = y[k]/L_U_mat[k,k]else:x[k] = (y[k] - np.sum(L_U_mat[k,k+1:]*x[k+1:]))/L_U_mat[k,k]return x

一般地，对于对角分解方法中A=LUA=LUA=LU当LLL的对角线元素为111的时候，称为Doolittle分解；当UUU的对角线元素为111的时候，称为Crout分解。

1.3 解对角方程的追赶法

有一种特殊的方程被称作是系数矩阵对角占优的线性方程组，设Ax=fAx=fAx=f，式子中
A=[b1c1a2b2c2:::an−1bn−1cn−1anbn],f=[f1f2:fn]A=\left [\begin{array}{cccc} b_{1}&c_{1}\\ a_{2}&b_{2}&c_{2}\\ &&:&:&:\\ &&&a_{n-1}&b_{n-1}&c_{n-1}\\ &&&&a_{n}&b_{n} \end{array}\right ], f=\left [\begin{array}{cccc} f_{1}\\ f_{2}\\ :\\ f_{n}\\ \end{array}\right ]A=⎣⎢⎢⎢⎢⎡b1a2c1b2c2::an−1:bn−1ancn−1bn⎦⎥⎥⎥⎥⎤,f=⎣⎢⎢⎡f1f2:fn⎦⎥⎥⎤
并且满足
{∣b1∣>∣c1∣>0∣bi∣≥∣ai∣+∣ci,∣(aici≠0)∣bn∣>∣an∣>0\begin{cases} |b_{1}|>|c_{1}|>0\\ |b_{i}|\geq |a_{i}|+|c_{i}&,|(a_{i}c_{i}\neq 0)\\ |b_{n}|>|a_{n}|>0 \end{cases}⎩⎪⎨⎪⎧∣b1∣>∣c1∣>0∣bi∣≥∣ai∣+∣ci∣bn∣>∣an∣>0,∣(aici=0)

我们对矩阵AAA进行Crout分解，即
A=LU=[l1a2l2a3l3::anln][1u11u21u3::1un]A=LU=\left [\begin{array}{cccc} l_{1}\\ a_{2}&l_{2}\\ &a_{3}&l_{3}\\ &&:&:\\ &&&a_{n}&l_{n} \end{array}\right ] \left [\begin{array}{cccc} 1&u_{1}\\ &1&u_{2}\\ &&1&u_{3}\\ &&&:&:\\ &&&&1&u_{n} \end{array}\right ]A=LU=⎣⎢⎢⎢⎢⎡l1a2l2a3l3::anln⎦⎥⎥⎥⎥⎤⎣⎢⎢⎢⎢⎡1u11u21u3::1un⎦⎥⎥⎥⎥⎤

这样解对角占优的方程过程如下所示：
{l1=b1ui=ci/lili+1=bi+1−ai+1ui,(i=2,...,n−1)\begin{cases} l_{1}=b_{1}\\ u_{i}=c_{i}/l_{i}\\ l_{i+1}=b_{i+1}-a_{i+1}u_{i}&,(i=2,...,n-1)\\ \end{cases}⎩⎪⎨⎪⎧l1=b1ui=ci/lili+1=bi+1−ai+1ui,(i=2,...,n−1)

回代过程为
{y1=f1/l1yi=(fi−aiyi−1)/li,(i=2,...,n)\begin{cases} y_{1}=f_{1}/l_{1}\\ y_{i}=(f_{i}-a_{i}y_{i-1})/l_{i}&,(i=2,...,n)\\ \end{cases}{y1=f1/l1yi=(fi−aiyi−1)/li,(i=2,...,n)
{xn=ynxi=yi−uixi+1,(i=n−1,...,1)\begin{cases} x_{n}=y_{n}\\ x_{i}=y_{i}-u_{i}x_{i+1}&,(i=n-1,...,1) \end{cases}{xn=ynxi=yi−uixi+1,(i=n−1,...,1)

1.4 误差分析

有一种线性方程组被称作是病态方程组，这类方程组有这样的一个特点。设线性方程为Ax=bAx=bAx=b，那么若AAA或者bbb有微小的扰动∣∣δA∣∣||\delta A||∣∣δA∣∣和∣∣δb∣∣||\delta b||∣∣δb∣∣时候，它的解xxx扰动很大∣∣δx∣∣||\delta x||∣∣δx∣∣，这样的方程组被称为是病态方程组，相应的矩阵称为病态矩阵。反之，如若∣∣δx∣∣||\delta x||∣∣δx∣∣比很小，称方程为良态方程。一般称为∣∣δx∣∣∣∣x∣∣\frac{||\delta x||}{||x||}∣∣x∣∣∣∣δx∣∣、∣∣δA∣∣∣∣A∣∣\frac{||\delta A||}{||A||}∣∣A∣∣∣∣δA∣∣和∣∣δb∣∣∣∣b∣∣\frac{||\delta b||}{||b||}∣∣b∣∣∣∣δb∣∣为相对扰动。线性方程有扰动的时候，方程则变为
(A+δA)(x+δx)=b+δb(A+\delta A)(x+\delta x)=b+\delta b(A+δA)(x+δx)=b+δb

(1) δA=0,δb≠0,b≠0\delta A=0,\delta b\neq 0,b\neq 0δA=0,δb=0,b=0时候，则方程化为：
A(x+δx)=b+δbA(x+\delta x)=b+\delta bA(x+δx)=b+δb

由于Ax=bAx=bAx=b则会有
Aδx=δbA\delta x=\delta bAδx=δb

由范数的定义可以得到
∣∣δx∣∣≤∣∣A−1∣∣⋅∣∣δb∣∣,∣∣b∣∣≤∣∣A∣∣⋅∣∣x∣∣||\delta x||\leq||A^{-1}||\cdot||\delta b||,||b||\leq||A||\cdot||x||∣∣δx∣∣≤∣∣A−1∣∣⋅∣∣δb∣∣,∣∣b∣∣≤∣∣A∣∣⋅∣∣x∣∣

所以得到
∣∣δx∣∣⋅∣∣b∣∣≤∣∣A∣∣⋅∣∣A−1∣∣⋅∣∣x∣∣⋅∣∣δb∣∣||\delta x||\cdot||b||\leq||A||\cdot||A^{-1}||\cdot||x||\cdot||\delta b||∣∣δx∣∣⋅∣∣b∣∣≤∣∣A∣∣⋅∣∣A−1∣∣⋅∣∣x∣∣⋅∣∣δb∣∣

从而有
∣∣δx∣∣∣∣x∣∣≤∣∣A∣∣⋅∣∣A−1∣∣⋅∣∣δb∣∣∣∣b∣∣\frac{||\delta x||}{||x||}\leq||A||\cdot||A^{-1}||\cdot{\frac{||\delta b||}{||b||}}∣∣x∣∣∣∣δx∣∣≤∣∣A∣∣⋅∣∣A−1∣∣⋅∣∣b∣∣∣∣δb∣∣

(2) δA≠0,δb=0,A≠0\delta A\neq0,\delta b=0,A\neq0δA=0,δb=0,A=0时候，则方程化为
(A+δA)(x+δx)=b(A+\delta A)(x+\delta x)=b(A+δA)(x+δx)=b

所以得到
δA(x+δx)+Aδx=0,δx=−AδA(x+δx)\delta A(x+\delta x)+A\delta x=0,\delta x=-A\delta A(x+\delta x)δA(x+δx)+Aδx=0,δx=−AδA(x+δx)

于是又有
∣∣δx∣∣≤∣∣A−1∣∣⋅∣∣δA∣∣⋅∣∣x+δx∣∣||\delta x||\leq||A^{-1}||\cdot||\delta A||\cdot||x+\delta x||∣∣δx∣∣≤∣∣A−1∣∣⋅∣∣δA∣∣⋅∣∣x+δx∣∣

从而有
∣∣δx∣∣∣∣x+δx∣∣≤∣∣A∣∣⋅∣∣A−1∣∣⋅∣∣δA∣∣∣∣A∣∣\frac{||\delta x||}{||x+\delta x||}\leq||A||\cdot||A^{-1}||\cdot{\frac{||\delta A||}{||A||}}∣∣x+δx∣∣∣∣δx∣∣≤∣∣A∣∣⋅∣∣A−1∣∣⋅∣∣A∣∣∣∣δA∣∣

定义条件数cond(A)=∣∣A∣∣⋅∣∣A−1∣∣\text{cond}(A)=||A||\cdot||A^{-1}||cond(A)=∣∣A∣∣⋅∣∣A−1∣∣。一般情况下，当系数矩阵的条件数较大的时候为病态方程，较小的时候为良态方程。
解决病态方程求解方法有以下的方法：
①采用高精度的运算方法计算；
②采用预处理的方法改善系数矩阵的条件数，例如可以寻找非奇异的对角矩阵P,QP,QP,Q将方程Ax=bAx=bAx=b化为
{PAQy=Pby=Q−1x\begin{cases} PAQy=Pb\\ y=Q^{-1}x\\ \end{cases}{PAQy=Pby=Q−1x

③当矩阵AAA的元素数量级较大的时候，对于系数矩阵AAA的行或者列乘以适当的比例因子，使得矩阵AAA所有的行列按照∞\infty∞范数大体上有相同的长度，使得系数矩阵AAA均衡，使得AAA条件数得到基本的改善。

2. 迭代法

一般情况下，对于低维度的系数矩阵求解可以直接使用高斯消去方法求解线性方程组。但是对于高维度的线性方程矩阵，求解显得力不从心。迭代方法是一种线性方程组的近似求解方法。所以这里提出了几种矩阵求解的迭代方法。

2.1 雅克比(Jacobi)迭代法

设线性方程组Ax=bAx=bAx=b，将系数矩阵AAA进行分解处理：A=D−L−UA=D-L-UA=D−L−U，其中矩阵D,L,UD,L,UD,L,U如下所示
D=[a11a22a33:ann]D=\left[ \begin{array}{cccc} a_{11}\\ &a_{22}\\ &&a_{33}\\ &&&:\\ &&&&a_{nn} \end{array}\right]D=⎣⎢⎢⎢⎢⎡a11a22a33:ann⎦⎥⎥⎥⎥⎤

−L=[0a210a31a320a41a42a43:::::an1an2an3...an(n−1)0]-L=\left[ \begin{array}{cccc} 0\\ a_{21}&0\\ a_{31}&a_{32}&0\\ a_{41}&a_{42}&a_{43}&:\\ :&:&:&:\\ a_{n1}&a_{n2}&a_{n3}&...&a_{n(n-1)}&0\\ \end{array}\right]−L=⎣⎢⎢⎢⎢⎢⎢⎡0a21a31a41:an10a32a42:an20a43:an3::...an(n−1)0⎦⎥⎥⎥⎥⎥⎥⎤

−U=[0a12a13...a1n0a23:a2n:::0a(n−1)n0]-U=\left[ \begin{array}{cccc} 0&a_{12}&a_{13}&...&a_{1n}\\ &0&a_{23}&:&a_{2n}\\ &&:&:&:\\ &&&0&a_{(n-1)n}\\ &&&&0\\ \end{array}\right]−U=⎣⎢⎢⎢⎢⎡0a120a13a23:...::0a1na2n:a(n−1)n0⎦⎥⎥⎥⎥⎤

则线性方程组化为(D−L−U)x=b(D-L-U)x=b(D−L−U)x=b，从而得到Dx=(L+U)x+bDx=(L+U)x+bDx=(L+U)x+b，于是得到
x=D−1(L+U)x+D−1b=(I−D−1A)x+D−1bx=D^{-1}(L+U)x+D^{-1}b=(I-D^{-1}A)x+D^{-1}bx=D−1(L+U)x+D−1b=(I−D−1A)x+D−1b

设BJ=I−D−1A,fJ=D−1bB_{J}=I-D^{-1}A,f_{J}=D^{-1}bBJ=I−D−1A,fJ=D−1b，可见通过以下的表达式进行递推可以实现方程的求解：
x(k+1)=BJx(k)+fJx^{(k+1)}=B_{J}x^{(k)}+f_{J}x(k+1)=BJx(k)+fJ

设线性方程组为∑j=1naijxj=bi,(i=1,2,...,n)\sum\limits_{j=1}^{n}a_{ij}x_{j}=b_{i},(i=1,2,...,n)j=1∑naijxj=bi,(i=1,2,...,n)，假设系数矩阵AAA是非奇异矩阵，并且aii≠0,(i=1,2,...,n)a_{ii}\neq 0,(i=1,2,...,n)aii=0,(i=1,2,...,n),得到的等价形式：
xi=1aii(bi−∑j=1,j≠inaijxj)x_{i}=\frac{1}{a_{ii}}(b_{i}-\sum\limits_{j=1,j\neq i}^{n}a_{ij}x_{j})xi=aii1(bi−j=1,j=i∑naijxj)

取初始向量x(0)=(x1(0),x2(0),...,xn(0))Tx^{(0)}=(x_{1}^{(0)},x_{2}^{(0)},...,x_{n}^{(0)})^{T}x(0)=(x1(0),x2(0),...,xn(0))T，使用迭代的方法可以得到迭代公式：
xi(k+1)=1aii(−∑j=1,j≠inaijxj(k)+bi)x_{i}^{(k+1)}=\frac{1}{a_{ii}}(-\sum\limits_{j=1,j\neq i}^{n}a_{ij}x_{j}^{(k)}+b_{i})xi(k+1)=aii1(−j=1,j=i∑naijxj(k)+bi)

这就是雅克比迭代法。

2.2 高斯-赛德尔迭代法

雅克比迭代法是一种同步迭代的方法，它计算过程中不能够充分利用计算过程中的变量。所以对于雅克比迭代方法提出了一种改进的迭代方法，如果每一次计算出一个新的分量便立即用它取代对应的旧的变量进行迭代，可能收敛更快，故而这种方法称作高斯-赛德尔迭代方法，迭代的公式为：
xi(k+1)=1aii(−∑j=1i−1aijxj(k+1)−∑j=inaijxj(k)+bi)x_{i}^{(k+1)}=\frac{1}{a_{ii}}(-\sum\limits_{j=1}^{i-1}a_{ij}x_{j}^{(k+1)}-\sum\limits_{j=i}^{n}a_{ij}x_{j}^{(k)}+b_{i})xi(k+1)=aii1(−j=1∑i−1aijxj(k+1)−j=i∑naijxj(k)+bi)

依然将系数矩阵进行分解A=D−L−UA=D-L-UA=D−L−U，将方程变为Dx=Lx+Ux+bDx=Lx+Ux+bDx=Lx+Ux+b。将新的分量代替旧分量，得到
Dx(k+1)=Lx(k+1)+Ux(k)+bDx^{(k+1)}=Lx^{(k+1)}+Ux^{(k)}+bDx(k+1)=Lx(k+1)+Ux(k)+b
即
(D−L)x(k+1)=Ux(k)+b(D-L)x^{(k+1)}=Ux^{(k)}+b(D−L)x(k+1)=Ux(k)+b

所以得到迭代表达式x(k+1)=(D−L)−1Ux(k)+(D−L)−1bx^{(k+1)}=(D-L)^{-1}Ux^{(k)}+(D-L)^{-1}bx(k+1)=(D−L)−1Ux(k)+(D−L)−1b，设系数BG−S=(D−L)−1UB_{G-S}=(D-L)^{-1}UBG−S=(D−L)−1U，fG−S=(D−L)−1bf_{G-S}=(D-L)^{-1}bfG−S=(D−L)−1b，所以得到
x(k+1)=BG−Sx(k)+fG−Sx^{(k+1)}=B_{G-S}x^{(k)}+f_{G-S}x(k+1)=BG−Sx(k)+fG−S

2.3 超松弛迭代法

逐次超松弛迭代法可以看作带参数的ω\omegaω的高斯-赛德尔迭代方法，是一种修正或加速，是求解大型稀疏矩阵的有效方法之一。
在高斯-赛德尔迭代方法中，由公式
xi(k+1)=1aii(−∑j=1i−1aijxj(k+1)−∑j=inaijxj(k)+bi)x_{i}^{(k+1)}=\frac{1}{a_{ii}}(-\sum\limits_{j=1}^{i-1}a_{ij}x_{j}^{(k+1)}-\sum\limits_{j=i}^{n}a_{ij}x_{j}^{(k)}+b_{i})xi(k+1)=aii1(−j=1∑i−1aijxj(k+1)−j=i∑naijxj(k)+bi)

得到
xi(k+1)=(1−ω)xi(k)+ωxi(k+1)=xi(k)+ω(xi(k+1)−xi(k))x_{i}^{(k+1)}=(1-\omega)x_{i}^{(k)}+\omega x_{i}^{(k+1)}=x_{i}^{(k)}+\omega(x_{i}^{(k+1)}-x_{i}^{(k)})xi(k+1)=(1−ω)xi(k)+ωxi(k+1)=xi(k)+ω(xi(k+1)−xi(k))

其中ω\omegaω称作是松弛因子。
现在有两种构造超松弛迭代算法，即基于G-S迭代方法和Jacobi迭代方法：
A. 基于G-S迭代方法
{xi(k+1)=xi(k)+ΔxiΔxi=ωri(k+1)ri(k+1)=1aii(−∑j=1i−1aijxj(k+1)−∑j=inaijxj(k)+bi)\begin{cases} x_{i}^{(k+1)}=x_{i}^{(k)}+\Delta x_{i}\\ \Delta x_{i}=\omega r_{i}^{(k+1)}\\ r_{i}^{(k+1)}=\frac{1}{a_{ii}}(-\sum\limits_{j=1}^{i-1}a_{ij}x_{j}^{(k+1)}-\sum\limits_{j=i}^{n}a_{ij}x_{j}^{(k)}+b_{i})\\ \end{cases}⎩⎪⎪⎪⎨⎪⎪⎪⎧xi(k+1)=xi(k)+ΔxiΔxi=ωri(k+1)ri(k+1)=aii1(−j=1∑i−1aijxj(k+1)−j=i∑naijxj(k)+bi)

显然ω=1\omega=1ω=1的时候，上式就是G-S迭代公式。
B. 基于Jacobi迭代方法
{xi(k+1)=xi(k)+ΔxiΔxi=ωri(k+1)ri(k+1)=1aii(−∑j=1naijxj(k)+bi)\begin{cases} x_{i}^{(k+1)}=x_{i}^{(k)}+\Delta x_{i}\\ \Delta x_{i}=\omega r_{i}^{(k+1)}\\ r_{i}^{(k+1)}=\frac{1}{a_{ii}}(-\sum\limits_{j=1}^{n}a_{ij}x_{j}^{(k)}+b_{i})\\ \end{cases}⎩⎪⎪⎪⎨⎪⎪⎪⎧xi(k+1)=xi(k)+ΔxiΔxi=ωri(k+1)ri(k+1)=aii1(−j=1∑naijxj(k)+bi)

一般地，ω\omegaω的取值决定着不同的松弛迭代方法：0<ω<10<\omega<10<ω<1时候称作低松弛方法；ω=1\omega=1ω=1时候为G-S方法或者是Jacobi方法；1<ω<21<\omega<21<ω<2时候是(渐次)超松弛方法。

2.4 迭代方法的收敛性

一般地，有迭代收敛方法的基本定理：假设有迭代方程组x=Bx+fx=Bx+fx=Bx+f，对于任意的初始向量x(0)x^{(0)}x(0)，迭代公式x(k+1)=Bx(k)+fx^{(k+1)}=Bx^{(k)}+fx(k+1)=Bx(k)+f收敛的充要条件是迭代矩阵BBB的谱半径ρ(B)<1\rho(B)<1ρ(B)<1。
设Ax=bAx=bAx=b，如果AAA是严格对角占优，则Jacobi迭代法、G-S迭代法均收敛。

对于迭代方法和直接求解的方法求解过程，以上具体代码可以参见笔者的github,上面有详细的计算过程。

小结

本小结详细介绍线性方程组的求解方法，在实际应用中应当合理利用这些方法来求解线性方程组来解决实际问题。

参考文献

[1] https://wenku.baidu.com/view/e1bca4b63968011ca2009117.html
[2] 计算方法，电子工业出版社

计算机求解线性方程组相关推荐

【原创】开源Math.NET基础数学类库使用(06)直接求解线性方程组
阅读目录前言 1.数值分析与线性方程 2.Math.NET解线性方程源码分析 3.Math.NET求解线性方程的实例 4.资源本博客所有文章分类的总目录:[总目录]本 ...
c语言解决方程的论文,c语言编程求解线性方程组论文1.doc
本科专业学年论文题目:线性方程组求解方法比较姓名付春雨专业计算机科学与技术专业班级 09级本科(1)班指导教师完成日期:2011 年 12月 8日题目:线性方 ...
用直接分解法求方程组的C语言程序,c语言编程求解线性方程组论文
计算机编程求解线性方程组第一章绪论在自然科学.工程技术.经济和医学各领域中产生的许多实际问题都可以通过数学语言描述为数学问题,也就是说,由实际问题建立数学模型,然后应用各种数学方法和技巧来求解 ...
数值计算大作业：Jacobi与Gauss -Seidel迭代求解线性方程组（Matlab实现）
作为研究生的入门课,数值计算的大作业算是所有研究生开学的重要编程作业. 我把Jacobi与Gauss -Seidel迭代求解线性方程组的数值计算作业在MATLAB中编程实现.具体的程序详细标注后放在文 ...
雅克比（Jacobi）迭代法求解线性方程组
长博文不利于翻阅,于是又将Jacobi迭代法单独出来了. 这篇博文把高斯-赛德尔迭代法和雅克比迭代法都放到一起了,个人觉得看着有点累.(迭代法求解线性方程组),不过还是要看的,因为它引出了迭代法. 进 ...
迭代法求解线性方程组的收敛问题总结
本讲之前,先将高斯-赛德尔迭代法和雅克比迭代法以及迭代法求解线性方程组贴出来,毕竟收敛问题研究的是迭代方法的收敛问题. 进入主题: 判断迭代法收敛的办法: 1.首先根据方程组的系数矩阵A的特点判断: ...
Gauss-Seidel迭代求解线性方程组
高斯-赛德尔迭代法考试比较多,所以考虑再三,还是单独提取出来独立一篇,方便查阅,突出重点. 首先举例引入: 通过手动求解下面的线性方程组得到精确解: 再用高斯-赛德尔迭代法求解比较: 本人拙见,将每一 ...
Guass列主元消去法求解线性方程组
上篇博文讲到:Guass消去法求解线性方程组,也提到了此方法求解线性方程组存在的问题,因此有如下: 基本思想: 看不懂没关系,直接看下面的例题,会有更直观的理解,然后再回到这里继续看. 选主元的步骤( ...
Guass消去法求解线性方程组
基本思想: 具体消元过程: 回代解方程组: 高斯顺序消去法求解线性方程组的计算公式: 存在的问题: 这种方法不是很实用,但是通常讲求解线性方程组都会先提到这种方法,这是基础,还是需要看看.

计算机求解线性方程组