上节课主要介绍了非线性方程的几种数值解法,其中包括交叉法(二分法、线性插值法)和开放法(牛顿法、割线法、固定点法)。本节课主要介绍线性方程组的数值求解方法,主要分为直接法和迭代法两类。直接法包括高斯消去法(Gauss elimination)、高斯约当法(Gauss-Jordan)以及LU分解法,迭代法包括Jacobi法和高斯塞德尔法(Gauss-Seidel) 。
1. 高斯消去法(Gauss Elimination Method)
高斯消去法是通过将线性方程组转化为上三角的形式,然后一步一步回代(back substitution)求解的方法。
例:$$ egin{cases}x_1+2x_2=5
3x_1+4x_2=6end{cases} $$
- 系统法:$$ egin{cases}x_1+2x_2=5
3x_1+4x_2=6end{cases} ightarrow egin{cases}x_1+2x_2=5
0-2x_2=-9end{cases} ightarrow egin{cases}x_1+2x_2=5
x_2=4.5end{cases} ightarrow egin{cases}x_1=-4
x_2=4.5end{cases}$$ - 系统法的增广矩阵形式:将方程写为矩阵的形式 $$egin{bmatrix}1 & 2 3 & 4end{bmatrix} egin{bmatrix}x_1 x_2 end{bmatrix}=egin{bmatrix}56end{bmatrix}$$ 系统法的求解过程可以写成如下增广矩阵的变换形式 $$egin{bmatrix}1 & 2 & 5 3 & 4 & 6end{bmatrix} ightarrow egin{bmatrix}1 & 2 & 5 & -2 & -9end{bmatrix} ightarrow egin{bmatrix}1 & 2 & 5 & 1 & 4.5end{bmatrix} ightarrow egin{bmatrix}1 & 0 & -4 & 1 & 4.5end{bmatrix} $$
例:
matrix:
类似的,(n imes n) 形式的方程组: $$ egin{cases}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{cases}$$
可以写作如下矩阵形式:$$A=egin{bmatrix}
a_{11} & a_{12} & cdots & a_{1n}
a_{21} & a_{22} & cdots & a_{2n}
vdots
a_{n1} & a_{n2} & cdots & a_{nn}
end{bmatrix}, x=egin{bmatrix}x_1 cdots x_nend{bmatrix}, b=egin{bmatrix}b_1 cdots b_nend{bmatrix}, Ax=b$$
高斯消去法的核心:
- 方程组写为矩阵形式
- 矩阵行变换:将一行中所有元素乘上一个标量,减去(加上)另外一行,该操作不会改变方程的解
- 逐行变换,直至将矩阵转换为便于求解的形式
一般而言,如下三种矩阵形式方便求解:
- 对角形式: (egin{bmatrix}1 & & 0\ & ddots & \0 & & 1 end{bmatrix}), (x_n=b_n)
- 上三角形式: (egin{bmatrix} a_{11} & a_{12} & cdots & a_{1n}\ & a_{22} & cdots & a_{2n}\ & & ddots&vdots\ & & & a_{nn}\ end{bmatrix}),使用回代法求解(back substitution): $$a_{nn}x_n=b_n ightarrow x_n = frac{b_n}{a_{nn}}$$ $$a_{n-1,n-1}x_{n-1}+ a_{n-1,n}x_n=b_{n-1} ightarrow x_{n-1}=frac{b_{n-1}-a_{n-1,n}x_n}{a_{n-1,n-1}} $$ 重复下去 $$x_i=frac{b_i-sum^n_{j=i+1}a_{ij}x_j}{a_{ii}}$$
- 下三角形式 :(egin{bmatrix} a_{11} & & &\ a_{21} & a_{22} & &\ vdots & vdots & ddots&\ a_{n1}& a_{n2}& cdots & a_{nn}\ end{bmatrix}) 使用前向替换法(forward substitution): $$a_{11}x_1=b_1 ightarrow x_1 = frac{b_1}{a_{11}}$$ $$a_{2,1}x_{1}+ a_{2,2}x_2=b_{2} ightarrow x_{2}=frac{b_{2}-a_{2,1}x_1}{a_{2,2}}$$ 重复下去 $$x_i=frac{b_i-sum^{i-1}{j=1}a{ij}x_j}{a_{ii}}$$
Matlab实现(求解(Ax=b), 使用左除(x=Aackslash b)):
>> A = [1 2;3 4]
A =
1 2
3 4
>> b = [5 6]'
b =
5
6
>> x = A
x =
-4.0000
4.5000
>> format long
>> x = A
x =
-3.999999999999999
4.499999999999999
MATLAB实现(高斯消去法):
function x=Gauss(a,b)
ab=[a,b];
[R,C]=size(ab);
for j=1:R-1
for i=j+1:R
ab(i,j:C)=ab(i,j:C)-ab(i,j)/ab(j,j)*ab(j,j:C);
end
end
x=zeros(R,1);
x(R)=ab(R,C)/ab(R,R);
for i=R-1:-1:1
x(i)=(ab(i,C)-ab(i,i+1:R)*x(i+1:R))/ab(i,i);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
>> A = [1 2;3 4]
A =
1 2
3 4
>> b = [5 6]'
b =
5
6
>> format long
>> Gauss(A,b)
ans =
-4.000000000000000
4.500000000000000
MATLAB实现(回代法与前向替换法)
function y = BackwardSub(a,b)
% The function solves a system of linear equations ax=b
% where a is upper triangular by using backward substitution.
% Input variables:
% a The matrix of coefficients.
% b A column vector of constants.
% Output variable:
% y A colum vector with the solution.
n = length(b);
y(n,1) = b(n)/a(n,n);
for i = n-1:-1:1
y(i,1)=(b(i)-a(i,i+1:n)*y(i+1:n,1))./a(i,i);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y = ForwardSub(a,b)
% The function solves a system of linear equations ax=b
% where a is lower triangular by using forward substitution.
% Input variables:
% a The matrix of coefficients.
% b A column vector of constants.
% Output variable:
% y A colum vector with the solution.
n = length(b);
y(1,1) = b(1)/a(1,1);
for i = 2:n
y(i,1)=(b(i)-a(i,1:i-1)*y(1:i-1,1))./a(i,i);
end
2. 高斯主元消去(Gauss Pivoting)
可以发现,当主元为(0)或非常小时,使用高斯消去法会失去精度,因此,在必要时可以交换行的顺序: $$egin{bmatrix}
0 & 2 & 5 3 & 4& 5
end{bmatrix}, $$ 高斯主元法就是通过交换行的顺序,来保证主元不为0,并且尽可能取绝对值较大的值,下面通过一个案例,说明主元很小时,高斯消去法会造成计算精度上出现较大误差,而高斯主元法的计算结果会更为精确。
例: $$ egin{cases}0.0003x_1 + 12.34x_2=12.343 .4321 x_1+x_2=5.321end{cases}$$ 使用高斯消去:
- (frac{0.4321}{0.0003}=1440) (round-off error, since (0.0003*1440 = 0.4320))
-
[egin{cases}0.0003x_1 + 12.34x_2=12.343 \ 0.0001 x_1+17770x_2=-17760end{cases} ightarrow x_2 = -frac{17760}{17770}=0.9994 ightarrow x_1 = frac{12.343-12.34*0.9994}{0.0003}=33.33 ]
MATLAB运算结果:
>> A = [0.0003 12.34;0.4321 1]
A =
0.000300000000000 12.340000000000000
0.432100000000000 1.000000000000000
>> b = [12.343 5.321]'
b =
12.343000000000000
5.321000000000000
>> A
ans =
10.000000000000000
1.000000000000000
若交换行序(高斯主元消去):
- (frac{0.0003}{0.4321}=0.0006943)(round-off error, since (0.4321*0.0006943 = 0.0003))
-
[egin{cases}0.4321 x_1+x_2=5.321 \ 0x_1 + 12.34x_2=12.34end{cases} ightarrow x_2 = -frac{12.34}{12.34}=1 ightarrow x_1 = frac{5.321-1}{0.4321}=10 ]
MATLAB实现(高斯主元消去)
function x=GaussPivot(a,b)
ab=[a,b];
[R,C]=size(ab);
for j=1:R-1
if ab(j,j)==0
for k=j+1:R
if ab(k,j)~=0
abTemp=ab(j,:);
ab(j,:)=ab(k,:);
ab(k,:)=abTemp;
break
end
end
end
for i=j+1:R
ab(i,j:C)=ab(i,j:C)-ab(i,j)/ab(j,j)*ab(j,j:C);
end
end
x=zeros(R,1);
x(R)=ab(R,C)/ab(R,R);
for i=R-1:-1:1
x(i)=(ab(i,C)-ab(i,i+1:R)*x(i+1:R))/ab(i,i);
end
end
3. 高斯约当法(Gauss-Jordan Method, 同时适用于求矩阵的逆)
高斯约当法一般采用高斯主元消去法,不同的是,高斯约当法每次会将主元规范化为(1),且最终将原矩阵转化为单位矩阵。
例:$$egin{bmatrix}4 & 1 & 2 & 21 2 & -2 & 2 & 8 1 & -2 & 4 & 16end{bmatrix}
ightarrow egin{bmatrix}1 & -2 & 4 & 16 2 & -2 & 2 & 8 4 & 1 & 2 & 21 end{bmatrix}
ightarrow egin{bmatrix}1 & -2 & 4 & 16 0 & 2 & -6 & -24 & 9 & -14 & -43 end{bmatrix}$$ $$
ightarrow egin{bmatrix}1 & -2 & 4 & 16 0 & 1 & -3 & -12 & 9 & -14 & -43 end{bmatrix}
ightarrow egin{bmatrix}1 & 0 & -2 & -8 0 & 1 & -3 & -12 & 0 & 13 & 65 end{bmatrix}
ightarrow egin{bmatrix}1 & 0 & -2 & -8 0 & 1 & -3 & -12 & 0 & 1 & 5 end{bmatrix}$$ $$
ightarrow egin{bmatrix}1 & 0 & 0 & 2 0 & 1 & 0 & 3 & 0 & 1 & 5 end{bmatrix}$$
可以看到,高斯约当法可以很方便地用来求解一个矩阵的逆:(AB=I), (Aegin{bmatrix}b_1 & b_2 & b_3 end{bmatrix} ightarrow egin{bmatrix}a_{11} & a_{12} & a_{13} & 1 & 0 & 0\ a_{21} & a_{22} & a_{23} & 0 & 1 & 0\ a_{31} & a_{32} & a_{33} & 0 & 0 & 1end{bmatrix} ightarrow egin{bmatrix}1 & 0 & 0 &c_{11} & c_{12} & c_{13} & \ 0 & 1 & 0 & c_{21} & c_{22} & c_{23} &\ 0 & 0 & 1& c_{31} & c_{32} & c_{33} end{bmatrix})
4. LU分解法(LU-decomposition Method)
LU分解是将一个矩阵分解为一个下三角矩阵和上三角矩阵的乘积:(A=LU), 其中(L) 表示下三角矩阵, (U)表示上三角矩阵。LU分解之后,求解线性方程(Ax=b) 就等价于求解如下两个方程:$$Ux=y,Ly=b$$
LU分解主要有两种方法来实现:
- 高斯消去(Gauss Elimination);
- Crout's Method
4.1 Crout's Method
例:
采用待定系数法求解:$$L_{11}=a_{11},L_{21}=a_{21},L_{31}=a_{31}$$ $$U_{12} = frac{a_{12}}{L_{11}}, U_{13}=frac{a_{13}}{L_{11}}$$ $$L_{22}=a_{22}-L_{21}U_{12},L_{32}=a_{32}-L_{31}U_{12}$$ $$U_{23}=frac{a_{23}-L_{21}U_{13}}{L_{22}}$$ $$L_{33}=a_{33}-L_{31}U_{13}-L_{32}U_{23}$$
结合计算过程,可以推导出,对于(n imes n)矩阵,使用Crout's Method进行LU分解的通用形式为:
- (L)矩阵的第一列 $$L_{i1}=a_{i1},i=1,2,...,n$$
- (U)矩阵的对角元素 $$U_{ii}=1,i=1,2,...,n$$
- (U)矩阵的第(1)行 $$U_{1j} = frac{a_{1j}}{L_{11}},j=2,...,n$$
- (L)矩阵的其他元素 $$L_{ij}=a_{ij}-sum^{j-1}{k=1}L{ik}U_{kj},i=2,3,...,n,j=2,3,...,i$$
- (U)矩阵的其他元素 $$U_{ij} = frac{a_{ij}-sum^{j-1}{k=1}L{ik}U_{kj}}{L_{ii}},i=2,3,...,n,j=i+1,i+2,...,n$$
MATLAB实现(LU分解,Crout's Method)
function [L, U] = LUdecompCrout(A)
% The function decomposes the matrix A into a lower triangular matrix L
% and an upper triangular matrix U, using Crout's method such that A=LU.
% Input variables:
% A The matrix of coefficients.
% Output variable:
% L Lower triangular matrix.
% U Upper triangular matrix.
[R, C] = size(A);
for i = 1:R
L(i,1) = A(i,1);
U(i,i) = 1;
end
for j = 2:R
U(1,j)= A(1,j)/L(1,1);
end
for i = 2:R
for j = 2:i
L(i,j)=A(i,j)-L(i,1:j-1)*U(1:j-1,j);
end
for j = i+1:R
U(i,j)=(A(i,j)-L(i,1:i-1)*U(1:i-1,j))/L(i,i);
end
end
4.2 高斯消去
general formulation: (m_{ij}=frac{a'_{ij}}{a'_{jj}})
例: (egin{bmatrix}1 & 2\3 & 4end{bmatrix}=egin{bmatrix}1 & 0\3 & 1end{bmatrix}egin{bmatrix}1 & 2\0 & -2end{bmatrix})
高斯消去法进行LU分解也可以通过待定系数法求解出矩阵参数:
- (U)矩阵第一行:$$U_{1j}=a_{1j},j=1,2,...,n$$
- (L)矩阵对角元素:$$L_{ii}=1,i=1,2,...,n$$
- (L)矩阵第一列:$$L_{i1}=frac{a_{i1}}{U_{11}}$$
- (U)矩阵其他元素:$$U_{ij}=a_{ij}-sum^{i-1}{k=1}L{ik}U_{kj},j=2,...,n,i=2,...,j$$
- (L)矩阵其他元素:$$L_{ij}=frac{a_{i,j}-sum^{j-1}{k=1}L{ik}U_{kj}}{L_{jj}},j=2,...,n,i=j+1,...,n$$
MATLAB实现(LU分解,高斯法,待定系数)
function [L,U] = LUdecopGauss(A)
[R,C] = size(A);
for i = 1:C
L(i,i) = 1;
U(1,i) = A(1,i);
end
for i = 2:R
L(i,1) = A(i,1)/U(1,1);
end
for j = 2:C
for i = 2:j
U(i,j) = A(i,j)-L(i,1:i-1)*U(1:i-1,j);
end
for i= (j+1):R
L(i,j) = (A(i,j)-L(i,1:j-1)*U(1:j-1,j))/U(j,j);
end
end
end
高斯消去的待定系数法与Crout's 方法的待定系数法相比,要难以理解,可读性差,主要在于通解形式与常规的求解顺序不一致。为便于理解,还有一种拉格朗日形式来实现高斯消去法:
回顾高斯消去的过程:
- 第一步运算 $$A = egin{bmatrix}
a_{11} & a_{12} & cdots & a_{1n}
a_{21} & a_{22} & cdots & a_{2n}
vdots
a_{n1} & a_{n2} & cdots & a_{nn}
end{bmatrix} ightarrow A'=egin{bmatrix}
a_{11} & a_{12} & cdots & a_{1n}
0 & a'{22} & cdots & a'{2n}
vdots
0 & a'{n2} & cdots & a'{nn}
end{bmatrix}=L_1A,L_1=egin{bmatrix}
1 & 0 & cdots & 0
-frac{a_{21}}{a_{11}} & 1 & cdots & 0
vdots
-frac{a_{n1}}{a_{11}} & 0 & cdots & 1
end{bmatrix}$$ 这是一步行变换 - 第二步运算 $$A'=egin{bmatrix}
a_{11} & a_{12} & cdots & a_{1n}
0 & a'{22} & cdots & a'{2n}
vdots
0 & a'{n2} & cdots & a'{nn}
end{bmatrix} ightarrow A''=egin{bmatrix}
a_{11} & a_{12} & cdots & a_{1n}
0 & a'{22} & cdots & a'{2n}
vdots
0 & 0 & cdots & a''{nn}
end{bmatrix}=L_2A‘,L_2=egin{bmatrix}
1 & 0 & cdots & 0
0 & 1 & cdots & 0
vdots
0 & -frac{a'{n2}}{a'_{22}} & cdots & 1
end{bmatrix}$$ - 第(n-1)步运算后 $$egin{bmatrix}U_{11} & U_{12} & cdots & U_{1n} & U_{22} & cdots & U_{2n} & & ddots & vdots & & & U_{nn}end{bmatrix}=L_{n-1}L_{n-2}...L_1A$$ 此时即为LU分解中的U矩阵,而(L=L^{-1}_{1}...L^{-1}_{n-2}L^{-1}_{n-1}=egin{bmatrix} 1 & 0 & cdots & 0\ frac{a_{21}}{a_{11}} & 1 & cdots & 0\ vdots\ frac{a_{n1}}{a_{11}} & 0 & cdots & 1 end{bmatrix}egin{bmatrix} 1 & 0 & cdots & 0\ 0 & 1 & cdots & 0\ vdots\ 0 & frac{a'_{n2}}{a'_{22}} & cdots & 1 end{bmatrix}...=egin{bmatrix} 1 & 0 & cdots & 0\ frac{a_{21}}{a_{11}} & 1 & cdots & 0\ vdots\ frac{a_{n1}}{a_{11}} & frac{a'_{n2}}{a'_{22}} & cdots & 1 end{bmatrix})可以看到 (L_{ij}=frac{a_{ij}}{a'_{jj}})
MATLAB实现(LU分解,高斯消去,拉格朗日形式)
function [L,U]=LUGauss2(A)
[R,C]=size(A);
L = zeros(size(A));
for i = 1:C
L(i,i) = 1;
end
for j=1:R-1
for i=j+1:R
L(i,j) = A(i,j)/A(j,j);
A(i,j:C)=A(i,j:C)-A(i,j)/A(j,j)*A(j,j:C);
end
end
U = A;
end
5. 总结
本节课主要介绍了求解线性方程组的几种直接法,包括高斯消去、高斯主元消去、高斯约当、LU分解法等,其中高斯约当法可以用来求矩阵的逆。 矩阵的LU分解又可以通过高斯消去和Crout 方法来实现,对参数矩阵进行LU分解后,很容易可以求解方程。下节课介绍求解线性方程组的迭代法。