高斯消去法
高斯消去法(Gaussian elimination)是求解线性方程组的一种常用方法。它的基本思想是通过一系列的行变换将增广矩阵(包含系数矩阵和常数矩阵)化为阶梯形矩阵,然后通过回代求解未知量。以下是高斯消去法的实现步骤:
设线性方程组为 A x = b Ax=b Ax=b,其中 A A A 为 n × n n\times n n×n 的系数矩阵, b b b 为 n × 1 n\times 1 n×1 的常数矩阵。
将增广矩阵 [ A ∣ b ] [A|b] [A∣b] 写成增广矩阵形式。
对矩阵 A A A 进行初等行变换,消元为主要步骤,步骤如下:
①. 将第一行首个非零元素作为主元素,用主元素消掉第一列中下面所有元素;
②. 将第二行首个非零元素作为主元素,用主元素消掉第二列中下面所有元素;
③. 重复此步骤,直至所有行变成上三角矩阵,即矩阵 U U U。
对
U
U
U 进行回带求解,步骤如下:
①. 对最后一个未知量
x
n
x_n
xn 进行求解;
②. 对第 k k k 个未知量 x k x_k xk 进行求解,其中 k = n − 1 , n − 2 , . . . , 1 k=n-1,n-2,...,1 k=n−1,n−2,...,1。
高斯消去法的公式式如下:
初等行变换:
用第 i i i 行的第一个非零元素将第 j j j 行消元: R j ← R j − a j i a i i R i R_j\leftarrow R_j-\dfrac{a_{ji}}{a_{ii}}R_i Rj←Rj−aiiajiRi。
回代求解:
x n = b n a n n x_n=\dfrac{b_n}{a_{nn}} xn=annbn, x k = 1 a k k ( b k − ∑ i = k + 1 n a k i x i ) x_k=\dfrac{1}{a_{kk}}(b_k-\sum_{i=k+1}^n a_{ki}x_i) xk=akk1(bk−∑i=k+1nakixi)。
通过这种方法可以求得线性方程组的解。
例如,求解如下线性方程组的解:
{ x 1 + x 2 + x 3 = 6 , 2 x 1 + 3 x 2 + 4 x 3 = 17 , 5 x 1 + 6 x 2 + 8 x 3 = 33. \begin{cases} x_1+x_2+x_3=6,\\ 2x_1+3x_2+4x_3=17,\\ 5x_1+6x_2+8x_3=33. \end{cases} ⎩ ⎨ ⎧x1+x2+x3=6,2x1+3x2+4x3=17,5x1+6x2+8x3=33.
应用高斯消去法,有如下步骤:
(1)写成增广矩阵形式:
[ 1 1 1 ∣ 6 2 3 4 ∣ 17 5 6 8 ∣ 33 ] \begin{bmatrix} 1 &1 &1 &| &6\\ 2 &3 &4 &| &17\\ 5 &6 &8 &| &33 \end{bmatrix} 125136148∣∣∣61733
(2)初等行变换:
[ 1 1 1 ∣ 6 0 1 2 ∣ 5 0 0 3 ∣ 3 ] \begin{bmatrix} 1 &1 &1 &| &6\\ 0 &1 &2 &| &5\\ 0 &0 &3 &| &3 \end{bmatrix} 100110123∣∣∣653
(3)回代求解:
x 3 = 1 , x 2 = 1 , x 1 = 4 x_3=1, x_2=1, x_1=4 x3=1,x2=1,x1=4
因此,该线性方程组的解为: x 1 = 4 , x 2 = 1 , x 3 = 1 x_1=4, x_2=1, x_3=1 x1=4,x2=1,x3=1。
以下是高斯消去法的Python代码实现:
import numpy as np
def gaussian_elimination(A, b):
"""
A: 系数矩阵
b: 常数矩阵
"""
n = len(b) # 线性方程组的未知数的数量
# 拼接增广矩阵
Ab = np.concatenate((A, b.reshape(n, 1)), axis=1)
# 初等行变换,将Ab变为上三角矩阵
for i in range(n-1):
# 将第i+1行到第n行的第(i+1)列消成0
for j in range(i+1, n):
Ab[j] -= Ab[i] * Ab[j][i] / Ab[i][i]
# 回代求解
x = np.zeros(n)
for i in range(n-1, -1, -1):
x[i] = (Ab[i][-1] - np.dot(Ab[i][i+1:n], x[i+1:n])) / Ab[i][i]
return x
这段代码是在做矩阵运算中的高斯消元算法,具体实现的是将矩阵中第i列的元素通过与第j列的运算,将第j列第i行位置的元素变为0。其中Ab表示增广矩阵,Ab[j][i]为第j行第i列元素,Ab[i][i]为第i行第i列元素,Ab[j][i] / Ab[i][i]表示第j行第i列元素除以第i行第i列元素, Ab[j] -= Ab[i] * Ab[j][i] / Ab[i][i] 则为第j行的各个元素分别减去第i行的各个元素乘以Ab[j][i] / Ab[i][i]的结果。
“Ab[i][i+1:n]”这一行代码的作用是截取了矩阵Ab中第i行的一部分,即从第i+1个元素到第n-1个元素的切片,用来进行高斯消元的计算(通常用于计算矩阵消元中对角线下方的元素)。在矩阵计算中,切片操作十分常见,可以快速地提取所需数据并进行针对性的运算处理。
高斯列主元消去法
高斯列主元消去法(Gaussian Elimination with Partial Pivoting)是一种用于解线性方程组的数值方法,它通过将系数矩阵变换为上三角矩阵,然后回代求解出未知数的值。
以下是该方法的公式、实现步骤和代码实现:
公式:
假设有一个n阶线性方程组 A x = b Ax=b Ax=b,其中 A A A 是一个可逆矩阵, x x x 和 b b b 是列向量。将 A A A 变换为上三角矩阵 U U U 的过程可以表示为以下一系列高斯消去步骤:
[ a 11 a 12 a 13 ⋯ a 1 n a 21 a 22 a 23 ⋯ a 2 n ⋮ ⋮ ⋮ ⋱ ⋮ a n 1 a n 2 a n 3 ⋯ a n n ] [ x 1 x 2 ⋮ x n ] [ b 1 b 2 ⋮ b n ] \begin{bmatrix} a_{11} & a_{12} & a_{13} & \ \cdots & a_{1n}\\ a_{21} & a_{22} & a_{23} & \ \cdots & a_{2n}\\ \vdots & \vdots & \vdots & \ \ddots & \vdots\\ a_{n1} & a_{n2} & a_{n3} & \ \cdots & a_{nn}\\ \end{bmatrix} \begin{bmatrix} x_{1}\\ x_{2}\\ \vdots\\ x_{n}\ \end{bmatrix} \begin{bmatrix} b_{1}\\ b_{2}\\ \vdots\\ b_{n}\\ \end{bmatrix} a11a21⋮an1a12a22⋮an2a13a23⋮an3 ⋯ ⋯ ⋱ ⋯a1na2n⋮ann x1x2⋮xn b1b2⋮bn
Step 1:
对于第1列 A A A,找到绝对值最大的元素 a i 1 a_{i1} ai1,即
i = arg max j = 1 n ∣ a j 1 ∣ i = \arg \max_{j=1}^{n} |a_{j1}| i=argj=1maxn∣aj1∣
如果 i = 1 i=1 i=1,则继续下一步;否则交换第1行和第 i i i 行,即
[ a 11 a 12 a 13 ⋯ a 1 n a 21 a 22 a 23 ⋯ a 2 n ⋮ ⋮ ⋮ ⋱ ⋮ a n 1 a n 2 a n 3 ⋯ a n n ] ∼ [ a i 1 a i 2 a i 3 ⋯ a i n a 21 a 22 a 23 ⋯ a 2 n ⋮ ⋮ ⋮ ⋱ ⋮ a n 1 a n 2 a n 3 ⋯ a n n ] \begin{bmatrix} a_{11} & a_{12} & a_{13} & \ \cdots & a_{1n}\\ a_{21} & a_{22} & a_{23} & \ \cdots & a_{2n}\\ \vdots & \vdots & \vdots & \ \ddots & \vdots\\ a_{n1} & a_{n2} & a_{n3} & \ \cdots & a_{nn}\\ \end{bmatrix} \sim \begin{bmatrix} a_{i1} & a_{i2} & a_{i3} & \ \cdots & a_{in}\\ a_{21} & a_{22} & a_{23} & \ \cdots & a_{2n}\\ \vdots & \vdots & \vdots & \ \ddots & \vdots\\ a_{n1} & a_{n2} & a_{n3} & \ \cdots & a_{nn}\\ \end{bmatrix} a11a21⋮an1a12a22⋮an2a13a23⋮an3 ⋯ ⋯ ⋱ ⋯a1na2n⋮ann ∼ ai1a21⋮an1ai2a22⋮an2ai3a23⋮an3 ⋯ ⋯ ⋱ ⋯aina2n⋮ann
Step 2:
对于第1列 A A A,将第 2 ∼ n 2 \sim n 2∼n 行乘以一个系数 m i 1 = − a i 1 a 11 m_{i1}=-\dfrac{a_{i1}}{a_{11}} mi1=−a11ai1 并加上第1行,即
[ a 11 a 12 a 13 ⋯ a 1 n a 21 a 22 a 23 ⋯ a 2 n ⋮ ⋮ ⋮ ⋱ ⋮ a n 1 a n 2 a n 3 ⋯ a n n ] ∼ [ a 11 a 12 a 13 ⋯ a 1 n 0 a 22 ( 1 ) a 23 ( 1 ) ⋯ a 2 n ( 1 ) ⋮ ⋮ ⋮ ⋱ ⋮ 0 a n 2 ( 1 ) a n 3 ( 1 ) ⋯ a n n ( 1 ) ] \begin{bmatrix} a_{11} & a_{12} & a_{13} &\ \cdots & a_{1n}\\ a_{21} & a_{22} & a_{23} &\ \cdots & a_{2n}\\ \vdots & \vdots & \vdots &\ \ddots & \vdots\\ a_{n1} & a_{n2} & a_{n3} &\ \cdots & a_{nn}\\ \end{bmatrix} \sim \begin{bmatrix} a_{11} & a_{12} & a_{13} &\ \cdots & a_{1n}\\ 0 & a_{22}^{(1)} & a_{23}^{(1)} & \cdots & a_{2n}^{(1)}\\ \vdots & \vdots & \vdots & \ \ddots & \vdots\\ 0 & a_{n2}^{(1)} & a_{n3}^{(1)} &\ \cdots & a_{nn}^{(1)}\\ \end{bmatrix} a11a21⋮an1a12a22⋮an2a13a23⋮an3 ⋯ ⋯ ⋱ ⋯a1na2n⋮ann ∼ a110⋮0a12a22(1)⋮an2(1)a13a23(1)⋮an3(1) ⋯⋯ ⋱ ⋯a1na2n(1)⋮ann(1)
其中 a i j ( 1 ) = a i j − m i 1 a 1 j a_{ij}^{(1)} = a_{ij} - m_{i1}a_{1j} aij(1)=aij−mi1a1j。
Step 3:
重复执行Step 1和Step 2,将矩阵 A A A 变换为上三角矩阵 U U U,即
[ a 11 a 12 a 13 ⋯ a 1 n 0 a 22 a 23 ⋯ a 2 n ⋮ ⋮ ⋮ ⋱ ⋮ 0 0 0 ⋯ a n n ] [ x 1 x 2 ⋮ x n ] [ b 1 b 2 ⋮ b n ] \begin{bmatrix} a_{11} & a_{12} & a_{13} & \cdots & a_{1n}\\ 0 & a_{22} & a_{23} & \cdots & a_{2n}\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & 0 & \cdots & a_{nn}\\ \end{bmatrix} \begin{bmatrix} x_{1}\\ x_{2}\\ \vdots\\ x_{n}\\ \end{bmatrix} \begin{bmatrix} b_{1}\\ b_{2}\\ \vdots\\ b_{n}\\ \end{bmatrix} a110⋮0a12a22⋮0a13a23⋮0⋯⋯⋱⋯a1na2n⋮ann x1x2⋮xn b1b2⋮bn
回代求解:
从后往前依次求解未知数 x i x_i xi,具体步骤如下:
初始化
x
n
=
b
n
a
n
n
x_n = \dfrac{b_n}{a_{nn}}
xn=annbn。
对于
i
=
n
−
1
,
n
−
2
,
⋯
,
1
i=n-1,n-2,\cdots,1
i=n−1,n−2,⋯,1,计算
x
i
=
1
a
i
i
(
b
i
−
∑
j
=
i
+
1
n
a
i
j
x
j
)
x_i = \frac{1}{a_{ii}} \left( b_i - \sum_{j=i+1}^{n} a_{ij}x_j \right)
xi=aii1(bi−j=i+1∑naijxj)
代码实现:
下面是使用Python实现高斯列主元消去法的示例代码。其中,函数gauss_elimination_with_partial_pivoting实现了高斯列主元消去法的核心算法,函数solve_linear_system则是调用该函数并输出解向量的代码。
import numpy as np
def Gaussian_row_elimination(A, b):
n = len(A)
for k in range(n):
# 找主元
max_row = k
for i in range(k+1, n):
if abs(A[i][k]) > abs(A[max_row][k]):
max_row = i
# 行列交换
A[k], A[max_row] = A[max_row], A[k]
b[k], b[max_row] = b[max_row], b[k]
# 消去变量
for i in range(k+1, n):
factor = A[i][k] / A[k][k]
for j in range(k+1, n):
A[i][j] -= factor * A[k][j]
b[i] -= factor * b[k]
# 回代
x = [0] * n
for i in range(n-1, -1, -1):
x[i] = (b[i] - sum(A[i][j] * x[j] for j in range(i+1, n))) / A[i][i]
return x
其中,A为线性方程组的系数矩阵,b为常数矩阵,n为方程组的阶数,x是一个用于存储方程组解的一维数组。在该代码中,我们首先使用列主元素法将系数矩阵进行了局部缩放和单位主元化,然后通过反向代回求解该方程组。在使用该代码时,请注意处理可能的奇异矩阵情况。
奇异矩阵:
在线性代数中,如果一个矩阵的行列式为零,则该矩阵被称为奇异矩阵,否则被称为非奇异矩阵。奇异矩阵的一个重要特点是它们的逆矩阵不存在。这是因为,如果一个矩阵的行列式为零,它的伴随矩阵也将具有非零的行列式,从而无法表示成系数矩阵的逆矩阵的形式。
在解线性方程组时,奇异矩阵会导致使用高斯消元法、LU分解等方法失败,因为这些方法需要求解矩阵的逆矩阵。如果一个线性方程组的系数矩阵是奇异的,这意味着该方程组的解可能不存在,或者存在无穷多个解。因此,我们需要在解线性方程组之前检查矩阵是否是奇异的,以便及时识别出无法求解的情况。
高斯全主元消去法
高斯全主元消元法是一种线性代数中求解线性方程组的一种算法,相较于其他消元方法,通过适当的选取主元素可以有效防止算法过程中数值误差的积累,提高求解结果的准确性。
算法步骤如下:
设线性方程组为 A x = b Ax=b Ax=b,其中 A A A 为 n × n n\times n n×n 的系数矩阵, b b b 为 n × 1 n\times 1 n×1 的常数向量。
对矩阵
A
A
A 进行行列变换,选取
A
A
A 的主元素,消元过程中将主元素所在行列的其他元素全部归零。
对向量
b
b
b 进行相同的行变换。
重复步骤 1.2 直到
A
A
A 被归一为对角矩阵。
依次回代求解方程组。
算法公式如下:
1.一般的矩阵消元公式:
m
i
j
=
a
i
j
(
k
)
a
k
k
(
k
)
,
i
=
k
+
1
,
.
.
.
,
n
;
j
=
k
+
1
,
.
.
.
,
n
m_{ij}=\frac{a_{ij}^{(k)}}{a_{kk}^{(k)}}, \quad i=k+1,...,n; j=k+1,...,n
mij=akk(k)aij(k),i=k+1,...,n;j=k+1,...,n
a
i
j
(
k
+
1
)
=
a
i
j
(
k
)
−
m
i
j
a
k
j
(
k
)
,
i
=
k
+
1
,
.
.
.
,
n
;
j
=
k
+
1
,
.
.
.
,
n
a_{ij}^{(k+1)}=a_{ij}^{(k)}-m_{ij}a_{kj}^{(k)}, \quad i=k+1,...,n; j=k+1,...,n
aij(k+1)=aij(k)−mijakj(k),i=k+1,...,n;j=k+1,...,n
b
i
(
k
+
1
)
=
b
i
(
k
)
−
m
i
k
b
k
(
k
)
,
i
=
k
+
1
,
.
.
.
,
n
b_i^{(k+1)}=b_i^{(k)}-m_{ik}b_k^{(k)}, \quad i=k+1,...,n
bi(k+1)=bi(k)−mikbk(k),i=k+1,...,n
2.全主元消元公式:
k
=
arg
max
k
≤
i
≤
n
∣
a
i
k
∣
k=\underset{k \leq i \leq n}{\operatorname{arg}\max} |a_{ik}|
k=k≤i≤nargmax∣aik∣
l
=
arg
max
k
≤
j
≤
n
∣
a
l
j
∣
l=\underset{k \leq j \leq n}{\operatorname{arg}\max} |a_{lj}|
l=k≤j≤nargmax∣alj∣
m
i
j
=
a
i
j
a
k
l
,
i
=
k
+
1
,
.
.
.
,
n
;
j
=
l
+
1
,
.
.
.
,
n
m_{ij}=\frac{a_{ij}}{a_{kl}},\quad i=k+1,...,n;j=l+1,...,n
mij=aklaij,i=k+1,...,n;j=l+1,...,n
a
i
j
=
a
i
j
−
m
i
j
a
k
j
,
i
=
k
+
1
,
.
.
.
,
n
;
j
=
l
+
1
,
.
.
.
,
n
a_{ij}=a_{ij}-m_{ij}a_{kj},\quad i=k+1,...,n;j=l+1,...,n
aij=aij−mijakj,i=k+1,...,n;j=l+1,...,n
b
i
=
b
i
−
m
i
k
b
k
,
i
=
k
+
1
,
.
.
.
,
n
b_{i}=b_{i}-m_{ik}b_{k},\quad i=k+1,...,n
bi=bi−mikbk,i=k+1,...,n
其中
m
i
j
m_{ij}
mij 为消元矩阵,k 和 l 分别为当前消元所在的行和列的主元素。
实现中需要注意,若主元素为 0 0 0,则需要进行行/列变换选取非零主元素,同时若方程组无解、有无穷多解的情况也需要进行处理。
以下是高斯全主元消元法的Python代码实现:
import numpy as np
def gauss_full_pivoting(A,B):
n = len(A)
# 转换为增广矩阵
AB = np.hstack([A, B.reshape(-1,1)])
# 记录列交换的情况
col_swaps = np.arange(n)
for i in range(n-1):
# 在子矩阵中寻找主元(最大元素)的索引
maxindex = np.argmax(np.abs(AB[i:, i:-1]))
# 将索引转换为行列表达式
k, j = divmod(maxindex, n-i-1)
# 将子矩阵的行列还原为原矩阵的行列
k += i
j += i
# 行列交换
AB[[i, k], :] = AB[[k, i], :]
AB[:, [i, j]] = AB[:, [j, i]]
col_swaps[[i, j]] = col_swaps[[j, i]]
# 消去过程
for m in range(i+1, n):
r = AB[m,i]/AB[i,i]
AB[m,i:-1] -= r*AB[i,i:-1]
AB[m,-1] -= r*AB[i,-1]
# 回代过程
x = np.zeros(n)
for i in range(n-1, -1, -1):
x[i] = (AB[i,-1] - AB[i,i+1:-1].dot(x[i+1:])) / AB[i,i]
# 还原解
for i in range(n-1, -1, -1):
if col_swaps[i] != i:
x[i], x[col_swaps[i]] = x[col_swaps[i]], x[i]
return x
为什么需要解的复原?
1.高斯全主元消元法是一种通过行变换和列变换将线性方程组转化为行阶梯形矩阵的方法。在得到行阶梯形矩阵之后,可以通过回代的方式求出方程组的解。
然而,由于进行了行变换和列变换,得到的行阶梯形矩阵并不是原始矩阵的简化形式,其中可能包含了一些额外的信息。因此,需要进行还原操作,将行阶梯形矩阵转化为原始矩阵的简化形式,以方便对原始方程组的分析。
具体来说,高斯全主元消元法中的还原操作包括两个方面:
行还原:将行阶梯形矩阵中的主元位置(即每一列第一个非零元素所在的行)上方的元素全部清零,从而得到一个上三角矩阵。
列还原:将列变换的信息还原,即通过列变换将上三角矩阵中每一列的主元移动到对角线上。
通过行还原和列还原操作,可以得到一个简化的矩阵,它的对角线上是原始方程组的系数矩阵的行阶梯形式,从而方便进行求解。
2.高斯列主元消元法是一种通过列变换和行变换将线性方程组转化为列阶梯形矩阵的方法。在得到列阶梯形矩阵之后,可以通过回代的方式求出方程组的解。
与高斯全主元消元法不同的是,高斯列主元消元法在进行列变换的时候,总是选择主元所在列的绝对值最大的元素作为主元。这样,列变换过程中主元所在的列不会发生变化,从而不需要进行列还原操作。
因此,在高斯列主元消元法中,只需要进行行还原操作,将列阶梯形矩阵中的主元位置(即每一行第一个非零元素所在的列)左侧的元素全部清零,从而得到一个左上角到右下角的对角线为主元的上三角矩阵。这个上三角矩阵即为原始方程组的简化形式,可以直接进行回代求解。
因此,高斯列主元消元法不需要解的还原,是由于在进行列变换时,选择主元的方式保证了主元所在列不会发生变化,从而省去了列还原操作。
高斯约当消去法
高斯约当消去法(Gauss-Jordan elimination)是用于求解线性方程组的一种算法,它通过一系列列变换将系数矩阵变为单位矩阵,从而得到方程组的解。以下是该算法的公式和步骤。
假设有一个
n
n
n 元线性方程组:
a
1
,
1
x
1
+
a
1
,
2
x
2
+
⋯
+
a
1
,
n
x
n
=
b
1
a
2
,
1
x
1
+
a
2
,
2
x
2
+
⋯
+
a
2
,
n
x
n
=
b
2
⋮
a
n
,
1
x
1
+
a
n
,
2
x
2
+
⋯
+
a
n
,
n
x
n
=
b
n
\begin{align} a_{1,1}x_1 + a_{1,2}x_2 + \dots + a_{1,n}x_n &= b_1 \\ a_{2,1}x_1 + a_{2,2}x_2 + \dots + a_{2,n}x_n &= b_2 \\ \vdots \\ a_{n,1}x_1 + a_{n,2}x_2 + \dots + a_{n,n}x_n &= b_n \\ \end{align}
a1,1x1+a1,2x2+⋯+a1,nxna2,1x1+a2,2x2+⋯+a2,nxn⋮an,1x1+an,2x2+⋯+an,nxn=b1=b2=bn
其中 a i , j a_{i,j} ai,j 和 b i b_i bi 都是已知的常数。
步骤:
- 构建增广矩阵,将系数矩阵和常数向量组成一个 ( n × n + 1 ) (n\times n+1) (n×n+1) 的矩阵 A:
A = [ a 1 , 1 a 1 , 2 … a 1 , n ∣ b 1 a 2 , 1 a 2 , 2 … a 2 , n ∣ b 2 ⋮ ⋮ ⋱ ⋮ ∣ ⋮ a n , 1 a n , 2 … a n , n ∣ b n ] A = \begin{bmatrix} a_{1,1} & a_{1,2} & \dots & a_{1,n} & | & b_1 \\ a_{2,1} & a_{2,2} & \dots & a_{2,n} &| & b_2 \\ \vdots & \vdots & \ddots & \vdots & | & \vdots \\ a_{n,1} & a_{n,2} & \dots & a_{n,n} &| & b_n \end{bmatrix} A= a1,1a2,1⋮an,1a1,2a2,2⋮an,2……⋱…a1,na2,n⋮an,n∣∣∣∣b1b2⋮bn
- 对 A A A 进行行变换,将矩阵 A A A 变为阶梯形矩阵(echelon form)。
- 选取最左边非零的列,将此列上方所有元素变为 0 0 0。
- 重复以上步骤,直到 A A A 变为阶梯形矩阵。
- 对 A A A 进行列变换,将阶梯形矩阵 A A A 变为行简化阶梯形矩阵(reduced row echelon form)。
- 每一行第一个非零元素为 1 1 1。
- 每个 1 1 1 下面都是 0 0 0。
- 每个 1 1 1 所在列的其他元素也都是 0 0 0。
- 对于行简化阶梯形矩阵 A A A 中的每个未知数,可以通过其它已知解求出。
- 为每个未知数选择适当的数值。
- 对于已知的符号不同的元素,我们选择不同的数值,例如对于
x
i
x_i
xi 的值
∣
x
i
∣
>
∣
x
i
+
1
∣
|x_i|>|x_{i+1}|
∣xi∣>∣xi+1∣ 时,我们选择
x
i
x_i
xi 的值为正值。
以下是消元的数学公式:
假设我们需要将矩阵
A
A
A 变为阶梯形矩阵。即将第
r
r
r 行以下所有行的第
k
k
k 个元素都变成
0
0
0,其中
k
≤
r
k\leq r
k≤r,则可以使用以下公式:
A
i
,
j
′
=
A
i
,
j
−
A
r
,
j
A
r
,
r
A
i
,
r
A_{i,j}^{'} = A_{i,j} - \frac{A_{r,j}}{A_{r,r}}A_{i,r}
Ai,j′=Ai,j−Ar,rAr,jAi,r
而对于行简化阶梯形矩阵
A
A
A,则可以使用以下公式将每个未知数表示为其他已知解的线性组合:
x
k
=
∑
j
=
k
+
1
n
(
−
A
k
,
j
x
j
)
+
b
k
A
k
,
k
x_k = \frac{\sum_{j=k+1}^{n}(-A_{k,j}x_{j})+b_k}{A_{k,k}}
xk=Ak,k∑j=k+1n(−Ak,jxj)+bk
高斯约旦消去法是线性方程组求解中的一种常用方法,其代码实现如下:
import numpy as np
def GaussJordan(A, b):
n = len(b)
for k in range(n):
# 找到该列中绝对值最大的元素及其行数
max_val = abs(A[k][k])
max_row = k
for i in range(k+1, n):
if abs(A[i][k]) > max_val:
max_val = abs(A[i][k])
max_row = i
# 如果对角线元素为0,则无法继续运算
if A[max_row][k] == 0:
raise ValueError("此矩阵为奇异矩阵,无法进行消去")
# 交换当前行和绝对值最大的行,确保对角线元素最大
A[k], A[max_row] = A[max_row], A[k]
b[k], b[max_row] = b[max_row], b[k]
# 将当前行消成对角线元素为1
for i in range(k+1, n):
factor = A[i][k]/A[k][k]
for j in range(k+1, n):
A[i][j] -= factor * A[k][j]
b[i] -= factor * b[k]
# 将剩余行消成对角线为0
for i in range(n):
if i != k:
factor = A[i][k]/A[k][k]
for j in range(k+1, n):
A[i][j] -= factor * A[k][j]
b[i] -= factor * b[k]
# 将对角线元素变为1
scale = A[k][k]
for i in range(k, n):
A[k][i] /= scale
b[k] /= scale
return b
其中,我们使用了numpy包来进行矩阵运算,函数的输入为矩阵 A A A和列向量 b b b。函数的输出是高斯约当列主元消去法运算后的解向量。
Doolittle分解法
Doolittle分解法是一种矩阵分解方法,将一个 n n n阶方阵 A A A分解成一个下三角矩阵 L L L和一个上三角矩阵 U U U的乘积: A = L U A=LU A=LU。
给定线性方程组 Ax = b (其中 A 是一个 n×n 的矩阵,x 和 b 都是 n×1 的向量),Doolittle分解法将矩阵 A 分解为下三角矩阵 L 和上三角矩阵 U 的乘积,即 A = LU。
则可以使用以下公式来解线性方程组:
-
首先,将 A 分解成 L 和 U,其中 L 是一个下三角矩阵,U 是一个上三角矩阵,且有 A = LU。
实现步骤:a. 将 L L L和 U U U的对角线元素初始化为1。
b. 针对 A A A的第一列,计算 U U U的第一行元素,以及 L L L的第2到 n n n行第一列元素
c. 针对 A A A的第二列到第 n n n列,每列都计算出 U U U的对应行和 L L L的对应列。
d. 在每一列的计算中, U U U中每个元素都需要借助 L L L中的相应元素计算。例如, U i j U_{ij} Uij是由 L i k L_{ik} Lik和 U k j U_{kj} Ukj计算得出的,其中 k k k可以是 1 1 1到 j − 1 j-1 j−1的任何值。Doolittle分解的公式如下:
L 和 U 的计算:U(i,k) = A(i,k) - ∑(j=1)^(k-1) L(i,j)*U(j,k)
L(i,k) = (A(i,k) - ∑(j=1)^(k-1) L(i,j)*U(j,k)) / U(k,k)
其中,1<=i,j<=n,1<=k<=n。
-
然后令 Ly = b,使用前向替换求解 y;接着令 Ux = y,使用后向替换求解 x。
具体公式为:
y 的计算:
y 1 = b 1 / l 11 y i = 1 l i i ( b i − ∑ j = 1 i − 1 l i j y j ) ( i = 2 , 3 , ⋯ , n ) \begin{aligned} y_1 &= b_1/l_{11} \\ y_i &= \frac{1}{l_{ii}}(b_i - \sum_{j=1}^{i-1}l_{ij}y_j) & & (i=2,3,\cdots,n) \\ \end{aligned} y1yi=b1/l11=lii1(bi−j=1∑i−1lijyj)(i=2,3,⋯,n)
x 的计算:
x n = y n / u n n x i = 1 u i i ( y i − ∑ j = i + 1 n u i j x j ) ( i = n − 1 , n − 2 , ⋯ , 1 ) \begin{aligned} x_n &= y_n/u_{nn} \\ x_i &= \frac{1}{u_{ii}}(y_i - \sum_{j=i+1}^{n}u_{ij}x_j) & & (i=n-1,n-2,\cdots,1) \end{aligned} xnxi=yn/unn=uii1(yi−j=i+1∑nuijxj)(i=n−1,n−2,⋯,1)
其中,lij 表示 L 矩阵的第 i 行第 j 列元素,uij 表示 U 矩阵的第 i 行第 j 列元素。
具体的Doolittle分解的实现代码示例如下:
import numpy as np
def Doolittle(A,b):
# Doolittle分解
n = A.shape[0]
L = np.zeros((n, n))
U = np.zeros((n, n))
for k in range(n):
U[k, k] = A[k, k] - np.sum(L[k, :k] * U[:k, k])
L[k, k] = 1
for i in range(k + 1, n):
U[k, i] = A[k, i] - np.sum(L[k, :k] * U[:k, i])
L[i, k] = (A[i, k] - np.sum(L[i, :k] * U[:k, k])) / U[k, k]
# 求解Ly=b和Ux=y
y = np.zeros(n)
x = np.zeros(n)
for i in range(n):
y[i] = b[i] - np.sum(L[i, :i] * y[:i])
for i in range(n - 1, -1, -1):
x[i] = (y[i] - np.sum(U[i, i + 1:] * x[i + 1:])) / U[i, i]
return x
向量、矩阵的范数
向量和矩阵的范数是用于测量它们的大小的数学概念。向量范数通常称为向量的长度或大小,矩阵范数类似于测量线性变换的大小。
下面是一些向量和矩阵的范数定义以及对应的公式:
向量范数:
- ∣ x ∣ p = ( ∑ i = 1 n ∣ x i ∣ p ) 1 p |x|_{ p} = \left(\sum_{i=1 }^{n} |x_i|^p\right)^{\frac{1}{p}} ∣x∣p=(∑i=1n∣xi∣p)p1
当 p = ∞ p=\infty p=∞时,这个范数称为最大范数或 L ∞ L_{\infty} L∞范数,表示向量中各个分量绝对值的最大值
矩阵范数:
-
行和范数:
∥ A ∥ ∞ = max 1 ≤ i ≤ n ∑ j = 1 n ∣ a i j ∣ \|\boldsymbol{A}\|_{\infty}=\max _{1 \leq i \leq n} \sum_{j=1}^n\left|a_{i j}\right| ∥A∥∞=1≤i≤nmaxj=1∑n∣aij∣ -
列和范数:
∥ A ∥ 1 = max 1 ≤ j ≤ n ∑ i = 1 n ∣ a i j ∣ \|\boldsymbol{A}\|_1=\max _{1 \leq j \leq n} \sum_{i=1}^n\left|a_{i j}\right| ∥A∥1=1≤j≤nmaxi=1∑n∣aij∣ -
谱范数:
∥ A ∥ 2 = ρ ( A T A ) \|\boldsymbol{A}\|_2=\sqrt{\rho\left(A^{\mathrm{T}} \boldsymbol{A}\right)} ∥A∥2=ρ(ATA)
满足的性质
向量范数和矩阵范数满足下列性质:
-
非负性质: ∣ x ∣ ≥ 0 \left|\boldsymbol{x}\right| ≥ 0 ∣x∣≥0,当且仅当 x = 0 \boldsymbol{x} = \boldsymbol{0} x=0 时等号成立。
-
齐次性质: ∣ k x ∣ = ∣ k ∣ ∣ x ∣ \left|k\boldsymbol{x}\right| = |k|\left|\boldsymbol{x}\right| ∣kx∣=∣k∣∣x∣ 对于所有 x ∈ R n , k ∈ R \boldsymbol{x} \in \mathbb{R}^n, k \in \mathbb{R} x∈Rn,k∈R 成立。
-
三角不等式: ∣ x + y ∣ ≤ ∣ x ∣ + ∣ y ∣ \left|\boldsymbol{x} + \boldsymbol{y}\right| \leq \left|\boldsymbol{x}\right| + \left|\boldsymbol{y}\right| ∣x+y∣≤∣x∣+∣y∣ 对于所有 x , y ∈ R n \boldsymbol{x,y} \in \mathbb{R}^n x,y∈Rn 成立。
对于第三个性质:
向量范数具有三角不等式的性质,即对于任意的向量x和y,有以下不等式成立:
||x + y|| ≤ ||x|| + ||y||
这意味着向量的长度满足三角形中任意两条边之和大于第三边的性质。具体地,如果将x和y表示为空间中的两个向量,则x + y表示从x到y的直线段的长度,而||x||和||y||则分别表示x和y向量的长度,因此上述不等式表明从起点到终点的路径长度不能短于直接从起点到终点的长度。
矩阵范数满足的性质(包含上述向量范数的三个性质):
- ∣ A B ∣ ≤ ∣ A ∣ ⋅ ∣ B ∣ \left| \boldsymbol{AB} \right| \leq \left| \boldsymbol{A} \right| \cdot \left| \boldsymbol{B} \right| ∣AB∣≤∣A∣⋅∣B∣
谱范数又称矩阵2-范数,是一种衡量矩阵大小的方法,定义为矩阵的最大奇异值。谱范数在机器学习、数据分析、信号处理等领域中广泛应用。
数学公式:
∣ A ∣ 2 = σ max ( A ) |A|2 = \sigma{\max}(A) ∣A∣2=σmax(A)
其中, ∣ A ∣ 2 |A|2 ∣A∣2表示矩阵A的谱范数, σ max ( A ) \sigma{\max}(A) σmax(A)表示矩阵A的最大奇异值。
实现步骤:
1.求出矩阵A的奇异值分解,得到矩阵U、V和对角矩阵 Σ \Sigma Σ;
2.计算矩阵A的谱范数 ∣ A ∣ 2 = σ max ( A ) = max j ( σ j ) |A|2 = \sigma{\max}(A) = \max_{j}(\sigma_j) ∣A∣2=σmax(A)=maxj(σj),其中 σ j \sigma_j σj表示矩阵A的第j个奇异值。
线性方程组的误差分析
设 A 非奇异,已知线性方程组:
A
x
⃗
=
b
⃗
A \vec{x}=\vec{b}
Ax=b
-
设 A A A 精确, b ⃗ \vec{b} b 有误差 δ b ⃗ \delta \vec{b} δb ,得到的解为 x ⃗ + δ x ⃗ \vec{x}+\delta \vec{x} x+δx ,即:
A ( x ⃗ + δ x ⃗ ) = b ⃗ + δ b ⃗ A(\vec{x}+\delta \vec{x})=\vec{b}+\delta \vec{b} A(x+δx)=b+δb
有右端项误差估计表达式:∥ δ x ⃗ ∥ ∥ x ⃗ ∥ ≤ ∥ A ∥ ⋅ ∥ A − 1 ∥ ⋅ ∥ δ b ⃗ ∥ ∥ b ⃗ ∥ \frac{\|\delta \vec{x}\|}{\|\vec{x}\|} \leq \| A\|\cdot\| A^{-1} \| \cdot \frac{\|\delta \vec{b}\|}{\|\vec{b}\|} ∥x∥∥δx∥≤∥A∥⋅∥A−1∥⋅∥b∥∥δb∥
-
设 b ⃗ \vec{b} b 精确, A A A 有误差 δ A \delta A δA, 得到的解为 x ⃗ + δ x ⃗ \vec{x}+\delta \vec{x} x+δx, 即:
( A + δ A ) ( x ⃗ + δ x ˉ ) = b ⃗ (A+\delta A)(\vec{x}+\delta \bar{x})=\vec{b} (A+δA)(x+δxˉ)=b
有误差估计表达式:
∥ δ x ˙ ∥ ∥ x ⃗ ∥ ≤ ∥ A − 1 ∥ ⋅ ∥ δ A ∥ 1 − ∥ A − 1 ∥ ⋅ ∥ δ A ∥ = ∥ A ∥ ⋅ ∥ A − 1 ∥ ⋅ ∥ δ A ∥ ∥ A ∥ 1 − ∥ A ∥ ⋅ ∥ A − 1 ∥ ⋅ ∥ δ A ∥ ∥ A ∥ \frac{\|\delta \dot{x}\|}{\|\vec{x}\|} \leq \frac{\left\|A^{-1}\right\| \cdot\|\delta A\|}{1-\left\|A^{-1}\right\| \cdot\|\delta A\|} =\frac{\|A\| \cdot\left\|A^{-1}\right\| \cdot \frac{\|\delta A\|}{\|A\|}}{1-\|A\| \cdot\left\|A^{-1}\right\| \cdot \frac{\|\delta A\|}{\|A\|}} ∥x∥∥δx˙∥≤1−∥A−1∥⋅∥δA∥ A−1 ⋅∥δA∥=1−∥A∥⋅∥A−1∥⋅∥A∥∥δA∥∥A∥⋅ A−1 ⋅∥A∥∥δA∥ -
设 A A A 有误差 δ A \delta A δA, b ⃗ \vec{b} b 有误差 δ b ⃗ \delta \vec{b} δb , 得到的解为 x ⃗ + δ x ⃗ \vec{x}+\delta \vec{x} x+δx, 即:
( A + δ A ) ( x ⃗ + δ x ˉ ) = b ⃗ + δ b ⃗ (A+\delta A)(\vec{x}+\delta \bar{x})=\vec{b}+\delta \vec{b} (A+δA)(x+δxˉ)=b+δb有误差估计表达式:
∥ δ x ⃗ ∥ ∥ x ⃗ ∥ ≤ cond ( A ) 1 − cond ( A ) ∥ δ A ∥ / ∥ A ∥ ( ∥ δ A ∥ ∥ A ∥ + ∥ δ b ⃗ ∥ ∥ b ⃗ ∥ ) \frac{\|\delta \vec{x}\|}{\|\vec{x}\|} \leq \frac{\operatorname{cond}(A)}{1-\operatorname{cond}(A)\|\delta A\| /\|A\|}\left(\frac{\|\delta A\|}{\|A\|}+\frac{\|\delta \vec{b}\|}{\|\vec{b}\|}\right) ∥x∥∥δx∥≤1−cond(A)∥δA∥/∥A∥cond(A)(∥A∥∥δA∥+∥b∥∥δb∥)
其中, ∥ A ∥ ⋅ ∥ A − 1 ∥ \|A\| \cdot\left\|A^{-1}\right\| ∥A∥⋅ A−1 是关键的误差放大因子, 称为 A A A 的条件数, 记为 cond ( A ) (A) (A)。越大则 A \boldsymbol{A} A 越病态,难得准确解。
常用条件数:
cond
(
A
)
1
\operatorname{cond}(A)_1
cond(A)1、cond
(
A
)
∞
(A)_{\infty}
(A)∞、cond
(
A
)
2
(A)_2
(A)2
特别地,若A对称,则:
cond
(
A
)
2
=
max
∣
λ
∣
min
∣
λ
∣
\operatorname{cond}(A)_2=\frac{\max |\lambda|}{\min |\lambda|}
cond(A)2=min∣λ∣max∣λ∣