MATLAB 求解线性方程组

• 前言
• 1 线性方程组的直接解法
• 1.1利用左除运算符的直接解法
• 1.2 LU分解法
• 1.2.1 LU分解法的基本思想
• 1.2.2 Matlab的LU分解函数
• 1.2.3 用LU分解法求解线性方程组
• 2 线性方程组的迭代解法
• 2.1 基本思想
• 2.2 雅可比（Jacobi）迭代法
• 2.3 高斯-赛德尔（Gauss-Seidel）迭代法
• 总结

前言

线性方程组计算机解法有直接法和迭代法两大类。

• 直接法
  用计算公式直接计算求出线性方程组解的方法。包括有高斯（Gauss）消元法、列主元消元法、全主元消元法、LU分解法（矩阵的三角分解法）、追赶法等。
• 迭代法
  用迭代公式通过迭代计算求出满足精度要求的线性方程组近似解的方法。包括有Jacobi迭代法、seidel迭代法（即Gauss-seidel迭代法）、sor法（超松弛迭代法）等。（迭代法是一种逐次逼近线性方程组解的方法）

1 线性方程组的直接解法

1.1利用左除运算符的直接解法

A x = b      ⇒ d e t ( A ) ≠ 0      x = A \ b Ax=b \;\;\xRightarrow{det(A) \neq 0} \;\;x=A\backslash b Ax=bdet(A)​=0 ​x=A\b
eg:
{ 2 x 1 + x 2 − 5 x 3 + x 4 = 13 x 1 − 5 x 2 + 7 x 4 = − 9 2 x 2 + x 3 − x 4 = 6 x 1 + 6 x 2 − x 3 − 4 x 4 = 0 \begin{cases} 2x_1+x_2-5x_3+x_4=13\\ x_1-5x_2+7x_4 = -9 \\ 2x_2+x_3-x_4=6\\x_1+6x_2-x_3-4x_4=0 \end{cases} ⎩⎪⎪⎪⎨⎪⎪⎪⎧​2x1​+x2​−5x3​+x4​=13x1​−5x2​+7x4​=−92x2​+x3​−x4​=6x1​+6x2​−x3​−4x4​=0​
demo：

A=[2 1 -5 1;1 -5 0 7;0 2 1 -1;1 6 -1 -4];
b=[13,-9,6,0]';
x=A\b

运行结果如下：

x =
 -66.555555555555543
  25.666666666666664
 -18.777777777777775
  26.555555555555550

1.2 LU分解法

1.2.1 LU分解法的基本思想

A x = b ⇒ A = L U L U x = b { L y = b U x = y Ax=b \quad \xRightarrow{A=LU}LUx=b \Longrightarrow \left\{\begin{matrix} Ly=b \\ Ux=y \end{matrix}\right. Ax=bA=LU ​LUx=b{Ly=bUx=y​

1.2.2 Matlab的LU分解函数

• [L,U]=lu(A):产生一个上三角矩阵U和一个变换形式的下三角矩阵L，使之满足A=LU。（！注意：矩阵A必须为方阵）
• [L,U,P]=lu(A):产生一个上三角阵U和一个下三角阵L以及一个置换矩阵P，使之满足PA=LU。（！注意：矩阵A也必须为方阵）

1.2.3 用LU分解法求解线性方程组

• A x = b L U x = b x = U \ ( L \ b ) Ax=b\longrightarrow LUx=b \longrightarrow x=U\backslash(L\backslash b) Ax=bLUx=bx=U\(L\b)
• A x = b P A x = P b L U x = P b x = U \ ( L \ p ∗ b ) Ax=b\longrightarrow PAx=Pb \longrightarrow LUx=Pb\longrightarrow x=U\backslash(L\backslash p*b) Ax=bPAx=PbLUx=Pbx=U\(L\p∗b)

eg:
{ 2 x 1 + x 2 − 5 x 3 + x 4 = 13 x 1 − 5 x 2 + 7 x 4 = − 9 2 x 2 + x 3 − x 4 = 6 x 1 + 6 x 2 − x 3 − 4 x 4 = 0 \begin{cases} 2x_1+x_2-5x_3+x_4=13\\ x_1-5x_2+7x_4 = -9 \\ 2x_2+x_3-x_4=6\\x_1+6x_2-x_3-4x_4=0 \end{cases} ⎩⎪⎪⎪⎨⎪⎪⎪⎧​2x1​+x2​−5x3​+x4​=13x1​−5x2​+7x4​=−92x2​+x3​−x4​=6x1​+6x2​−x3​−4x4​=0​
demo：

A=[2 1 -5 1;1 -5 0 7;0 2 1 -1;1 6 -1 -4];
b=[13,-9,6,0]';
[L,U]=lu(A);
x=U\(L\b)
% [L,U,P]=lu(A);
% x=U\(L\P*b)

运行结果如下：

x =
 -66.555555555555543
  25.666666666666664
 -18.777777777777775
  26.555555555555550

2 线性方程组的迭代解法

2.1 基本思想

将线性方程组 Ax=b 等价变形为 x=Bx+g,然后构造向量迭代公式:
x ( k + 1 ) = B x ( k ) + g k = 0 , 1 , 2... x^{(k+1)}=Bx^{(k)}+g \quad k=0,1,2... x(k+1)=Bx(k)+gk=0,1,2...
再给定一个初始向量 x ( 0 ) x^{(0)} x(0),代入迭代公式计算出迭代向量序列：
x ( 1 ) , x ( 2 ) , . . . x ( k ) , . . . x^{(1)},x^{(2)},...x^{(k)},... x(1),x(2),...x(k),...
eg:
{ 10 x 1 − x 2 = 9 − x 1 + 10 x 2 − 2 x 3 = 7 − 2 x 2 + 10 x 3 = 6 { x 1 = 10 x 2 − 2 x 3 − 7 x 2 = 10 x 1 − 9 x 3 = ( 6 + 2 x 2 ) / 10 \begin{cases} 10x_1-x_2=9\\ -x_1+10x_2-2x_3 = 7\\ -2x_2+10x_3=6\\ \end{cases} \quad \Longrightarrow \quad\begin{cases} x_1=10x_2-2x_3-7\\ x_2=10x_1-9\\ x_3=(6+2x_2)/10\\ \end{cases} ⎩⎪⎨⎪⎧​10x1​−x2​=9−x1​+10x2​−2x3​=7−2x2​+10x3​=6​⎩⎪⎨⎪⎧​x1​=10x2​−2x3​−7x2​=10x1​−9x3​=(6+2x2​)/10​
{ x 1 ( k + 1 ) = 10 x 2 ( k ) − 2 x 3 ( k ) − 7 x 2 ( k + 1 ) = 10 x 1 ( k ) − 9 x 3 ( k + 1 ) = ( 6 + 2 x 2 ( k ) ) / 10 \Longrightarrow \begin{cases} x_1^{(k+1)}=10x_2^{(k)}-2x_3^{(k)}-7\\ x_2^{(k+1)}=10x_1^{(k)}-9\\ x_3^{(k+1)}=(6+2x_2^{(k)})/10\\ \end{cases} ⎩⎪⎨⎪⎧​x1(k+1)​=10x2(k)​−2x3(k)​−7x2(k+1)​=10x1(k)​−9x3(k+1)​=(6+2x2(k)​)/10​

2.2 雅可比（Jacobi）迭代法

A x = b ⇒ ( A = D − L − U ) ( D − L − U ) x = b Ax=b \quad \xRightarrow{(A=D-L-U)} (D-L-U)x=b Ax=b(A=D−L−U) ​(D−L−U)x=b
D = [ a 11 a 22 ⋱ a n n ] L = − [ 0 a 21 0 ⋮ ⋱ ⋱ a n 1 ⋯ a n , n − 1 0 ] U = − [ 0 a 12 ⋯ a 1 n 0 ⋱ ⋮ ⋱ a n − 1 , n 0 ] D=\begin{bmatrix} a_{11} & \\ & a_{22}\\ & & \ddots & \\ &&& a_{nn} \end{bmatrix}\quad L=-\begin{bmatrix} 0 & \\ a_{21}&0\\ \vdots & \ddots & \ddots & \\a_{n1} &\cdots &a_{n,n-1}& 0 \end{bmatrix}\quad U=-\begin{bmatrix} 0 & a_{12}& \cdots & a_{1n} \\ & 0& \ddots &\vdots\\ & &\ddots &a_{n-1,n}& \\ &&& 0 \end{bmatrix} D=⎣⎢⎢⎡​a11​​a22​​⋱​ann​​⎦⎥⎥⎤​L=−⎣⎢⎢⎢⎡​0a21​⋮an1​​0⋱⋯​⋱an,n−1​​0​⎦⎥⎥⎥⎤​U=−⎣⎢⎢⎢⎡​0​a12​0​⋯⋱⋱​a1n​⋮an−1,n​0​​⎦⎥⎥⎥⎤​
求解公式为：
x = D − 1 ( L + U ) + D − 1 b x=D^{-1}(L+U)+D^{-1}b x=D−1(L+U)+D−1b
迭代公式为：
x ( k + 1 ) = D − 1 ( L + U ) x ( k ) + D − 1 b      ⇒ B = D − 1 ( L + U ) , g = D − 1 b      x ( k + 1 ) = B x ( k ) + g x^{(k+1)}=D^{-1}(L+U)x^{(k)}+D^{-1}b \;\;\xRightarrow{B=D^{-1}(L+U),g=D^{-1}b}\;\; x^{(k+1)}=Bx^{(k)} +g x(k+1)=D−1(L+U)x(k)+D−1bB=D−1(L+U),g=D−1b ​x(k+1)=Bx(k)+g
代码实现：
（雅可比迭代法的函数文件Jacobi.m）

function [y,n]=jacobi(A,b,x0,ep) %A:系数矩阵；b:常数矩阵；x0:迭代初值；ep:迭代精度；
D=diag(diag(A));	%取A的主对角线上的元素建立对角矩阵；
L=-tril(A,-1);		%取A的主对角线以下（不包括主对角线）的元素建立下三角矩阵；
U=-triu(A,1);		%取A的主对角线以上（不包括主对角线）的元素建立上三角矩阵；
B=D\(L+U);
f=D\b;
y=B*x0+g;
n=1;
while norm(y-x0)>=ep	%用范数判断误差是否满足精度要求
    x0=y;
    y=B*x0+g;
    n=n+1;
end

2.3 高斯-赛德尔（Gauss-Seidel）迭代法

D x ( k + 1 ) = ( L + U ) x ( k ) + b       D x ( k + 1 ) = L x k + 1 + U x ( k ) + b       ( D − L ) x ( k + 1 ) = U x ( k ) + b Dx^{(k+1)}=(L+U)x^{(k)}+b\;\Longrightarrow \;Dx^{(k+1)}=Lx^{k+1}+Ux^{(k)}+b \;\Longrightarrow \;(D-L)x^{(k+1)}=Ux^{(k)}+b Dx(k+1)=(L+U)x(k)+bDx(k+1)=Lxk+1+Ux(k)+b(D−L)x(k+1)=Ux(k)+b
x ( k + 1 ) = ( D − L ) − 1 U x ( k ) + ( D − L ) − 1 b    → B = ( D − L ) − 1 U , g = ( D − L ) − 1 b    x ( k + 1 ) = B x ( k ) + g x^{(k+1)}=(D-L)^{-1}Ux^{(k)}+(D-L)^{-1}b\;\xrightarrow{B=(D-L)^{-1}U,g=(D-L)^{-1}b}\;x^{(k+1)}=Bx^{(k)}+g x(k+1)=(D−L)−1Ux(k)+(D−L)−1bB=(D−L)−1U,g=(D−L)−1b ​x(k+1)=Bx(k)+g
代码实现：
（Gauss-Seidel迭代法的函数文件gauseidel.m）

function [y,n]=gauseidel(A,b,x0,ep)     %A：系数矩阵；b:常数矩阵；x0:迭代初值；ep：迭代精度；
D=diag(diag(A));
L=-tril(A,-1);
U=-triu(A,1);
B=(D-L)\U;
f=(D-L)\b;
y=B*x0+g;
n=1;
while norm(y-x0)>=ep	%用范数判断误差是否满足精度要求
    x0=y;
    y=B*x0+g;
    n=n+1;
end

eg1:
分别用雅可比迭代法和高斯-赛德尔迭代法求解线性方程组。设迭代初值为0，迭代精度为 1 0 − 6 10^{-6} 10−6。
[ 4 − 2 1 − 2 4 3 − 1    − 3 3 ] [ x 1 x 2 x 3 ] = [ 1 5 0 ] \begin{bmatrix} 4 \quad-2\quad1 \\ -2\quad4\quad3 \\ -1\;-3\quad3 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2\\x_3 \end{bmatrix}=\begin{bmatrix} 1 \\ 5\\0 \end{bmatrix} ⎣⎡​4−21−243−1−33​⎦⎤​⎣⎡​x1​x2​x3​​⎦⎤​=⎣⎡​150​⎦⎤​
demo:

clc,clear;
A=[4 -2 -1;-2 4 3;-1 -3 3];
b=[1 5 0]';
[x,n]=jacobi(A,b,[0 0 0]',1.0e-6)
[x1,n1]=gauseidel(A,b,[0 0 0]',1.0e-6)

运行结果如下：

x =

   0.970588483091126
   0.852941039929380
   1.176470992801838
n =

    35
x1 =

   0.970588105856975
   0.852941168882785
   1.176470537501776
n1 =

    16

eg2:
分别用雅可比迭代法和高斯-赛德尔迭代法求解下列线性方程组，看是否收敛。

[ 1 2    − 2 1 1 1 2 2 1 ] [ x 1 x 2 x 3 ] = [ 9 7 6 ] \begin{bmatrix} 1 \quad2\;-2 \\ 1\quad1\quad1\\ 2\quad2\quad1 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2\\x_3 \end{bmatrix}=\begin{bmatrix} 9 \\ 7\\6 \end{bmatrix} ⎣⎡​12−2111221​⎦⎤​⎣⎡​x1​x2​x3​​⎦⎤​=⎣⎡​976​⎦⎤​
demo:

clc,clear;
A1=[1 2 -2;1 1 1;2 2 1];
b1=[9;7;6];
[x11,n11]=jacobi(A1,b1,[0 0 0]',1.0e-6)
[x22,n22]=gauseidel(A1,b1,[0 0 0]',1.0e-6)

运行结果如下：

x11 =
   -27
    26
     8
n11 =
     4
x22 =
   NaN
   NaN
   NaN
n22 =
    1012

总结

• 直接法：以矩阵初等变换为基础，可以求得方程的精确解；占用的内存空间大、程序实现较复杂；一般适合求解低阶稠密线性方程组。
• 迭代法：从给定的初始值逐步逼近精确解的过程，求解过程占用存储空间小，程序设计简单；适用于求解大型稀疏矩阵线性方程组；要考虑算法的收敛性。

补充
迭代法收敛性判断：

• 迭代收敛定理:    x ( k + 1 ) = B x ( k ) + g    对 任 意 x ( 0 ) 都 收 敛           ρ ( B ) < 1 ( ρ ( B ) 为 矩 阵 B 的 谱 半 径 ) \;x^{(k+1)}=Bx^{(k)}+g\;对任意x^{(0)}都收敛\;\iff\;\rho(B)<1\quad\quad(\rho(B)为矩阵B的谱半径) x(k+1)=Bx(k)+g对任意x(0)都收敛ρ(B)<1(ρ(B)为矩阵B的谱半径)
• 迭代收敛判别条件：

如果有用，请记得给我点赞！（^ 0^）\(//∇//)\