- 前言
- LU直接分解法 — Doolittle_Decomposition
- 代码实现
- test example
- 三对角矩阵的特殊解法
- 代码实现
- test example
- 迭代法解线性方程组
- Jacobi迭代法
- 代码实现
- test example
- Gauss-Seidel迭代法
- 代码实现
- test example
- 超松弛迭代法SOR
- 代码实现
- test example
- 整合Jacobi、Gauss_Seidel和SOR超松弛迭代法
- 写在最后
首先很多内容都跟矩阵论的部分重合了，相关的关于LU分解、LDV分解、LU解线性方程组、求行列式、求逆等都可以在矩阵论专栏中对应查看，这里仅补充一些其他方法。
LU直接分解法 — Doolittle_Decomposition
- 道立特分解属于一种先设定好U矩阵第一行和L矩阵第一列的值，便可以开始迭代求解L和U的其他值的算法，分解结果如下（以三阶矩阵为例）：
[ a 11 a 12 a 13 a 21 a 22 a 23 a 31 a 32 a 33 ] = [ 1 0 0 l 21 1 0 l 31 l 32 1 ] [ u 11 u 12 u 13 0 u 22 u 23 0 0 u 33 ] \begin{bmatrix} a_{11} & a_{12} & a_{13} \ a_{21} & a_{22} & a_{23} \ a_{31} & a_{32} & a_{33} \end{bmatrix}= \begin{bmatrix} 1 & 0 & 0 \ l_{21} & 1& 0 \ l_{31} & l_{32} & 1 \end{bmatrix} \begin{bmatrix} u_{11} & u_{12} & u_{13} \ 0 & u_{22} & u_{23} \ 0 & 0 & u_{33} \end{bmatrix} ⎣⎡a11a21a31a12a22a32a13a23a33⎦⎤=⎣⎡1l21l3101l32001⎦⎤⎣⎡u1100u12u220u13u23u33⎦⎤ - 算法流程（先求解U的对应行再求解L的对应列）：
- 求解U的第i行第j个元素： u i j = a i j − ∑ k = 1 i − 1 l i k ∗ u k j u_{ij}=a_{ij}-\sum_{k=1}^{i-1}l_{ik}*u_{kj} uij=aij−∑k=1i−1lik∗ukj
- 求解L的第j列第i个元素： l j i = ( a j i − ∑ k = 1 i − 1 l j k ∗ u k i ) / u i i l_{ji}=(a_{ji}-\sum_{k=1}^{i-1}l_{jk}*u_{ki})/u_{ii} lji=(aji−∑k=1i−1ljk∗uki)/uii
@staticmethod
def domain_row_transform(arr, eps=1e-6, Test=False):
"""domain_row_transform 选取主元并进行行变换
Args:
arr ([np.darray]): [input matrix]
eps ([float]): numerical threshold.
test (bool, optional): [show checking information]. Defaults to False.
Returns:
[A]: [resort one of arr]
Author: Junno
Date: 2022-04-27
"""
assert len(arr.shape) == 2
A = np.copy(arr)
M, N = A.shape
for i in range(M):
for j in range(i, N):
if Test:
print('During process in row:{}, col:{}'.format(i, j))
if sum(abs(A[i:, j])) > eps:
zero_index = []
non_zero_index = []
for k in range(i, M):
if abs(A[k, j]) > eps:
non_zero_index.append(k)
else:
zero_index.append(k)
non_zero_index = sorted(
non_zero_index, key=lambda x: abs(A[x, j]))
sort_idnex = non_zero_index+zero_index
if Test:
print('Sorting index for {} cols:\n'.format(i), sort_idnex)
if sort_idnex[0] != i:
A[[sort_idnex[0], i], :] = A[[i, sort_idnex[0]], :]
if Test:
print('After resorting\n', A)
return A
@classmethod
def Doolittle_Decomposition(self, target, domain_selection=True, test=False):
"""Doolittle_Decomposition equal to LU_Decomposition
Args:
target ([np.darray]): [input matrix]
domain_selection (bool, optional): [operate domain selection on target]. Defaults to True.
test (bool, optional): [show testing information]. Defaults to False.
Returns:
[resort input, L, U] if domain_selection=True else [L,U]
Author: Junno
Date: 2022-04-27
"""
assert len(target.shape) == 2
arr = np.copy(target)
M, N = arr.shape
L = np.zeros_like(arr, dtype=arr.dtype)
U = np.zeros_like(arr, dtype=arr.dtype)
if domain_selection:
arr = self.domain_row_transform(arr)
for i in range(M):
if test:
print(i)
print(L, U, sep='\n')
for j in range(i, N):
if test:
print('U', i, j, arr[i, j], L[i, :i], U[:i, j])
U[i, j] = arr[i, j]-L[i, :i]@U[:i, j]
if i == j:
L[i, j] = 1
else:
if test:
print('L', j, i, arr[j, i],
L[j, :i], U[:i, i], U[i, i])
L[j, i] = (arr[j, i]-(L[j, :i]@U[:i, i]))/U[i, i]
if domain_selection:
print('return resort input matrix, L, U')
return arr, L, U
else:
return L, U
- Tips：为避免大数吃小数的问题和舍入误差问题，加入了选主元的选项，可以对原输入矩阵进行行变换，不改变结果的正确性，可以自行选择
A=np.random.randn(5,5)
A
>>>
array([[-1.50948856, -1.61033161, -0.20671825, 0.4286645 , -0.2572608 ],
[-1.28615783, 0.87991375, 1.71080037, -1.61832726, 1.52242031],
[-0.98092031, -0.32273876, 2.20229521, 0.853808 , -2.15168542],
[-2.02003405, 0.05052931, -1.01116616, 0.26437957, -1.65870163],
[-0.83327297, 0.32495716, -1.02878403, -0.96607898, -0.59284746]])
L,U=Doolittle_Decomposition(A,domain_selection=False,test=False)
L
>>>
[[ 1. 0. 0. 0. 0. ]
[ 4.688 1. 0. 0. 0. ]
[ 5.484 -3863.106 1. 0. 0. ]
[ 8.449 -8279.836 2.143 1. 0. ]
[ -8.763 10132.862 -2.623 -11.842 1. ]]
U
>>>
[[ -0.165 0.618 -0.885 -0.515 -1.011]
[ 0. 0.001 3.434 2.532 6.176]
[ 0. 0. 13269.089 9782.795 23865.142]
[ 0. 0. 0. 0.216 -4.26 ]
[ 0. 0. 0. 0. -43.798]]
print(L@U)
>>>
array([[-1.50948856, -1.61033161, -0.20671825, 0.4286645 , -0.2572608 ],
[-1.28615783, 0.87991375, 1.71080037, -1.61832726, 1.52242031],
[-0.98092031, -0.32273876, 2.20229521, 0.853808 , -2.15168542],
[-2.02003405, 0.05052931, -1.01116616, 0.26437957, -1.65870163],
[-0.83327297, 0.32495716, -1.02878403, -0.96607898, -0.59284746]])
-
对于一般的线性方程组，做高斯消去法的乘法复杂度大概为 O ( N 3 ) O(N^3) O(N3)， 当N超过100时就比较吃力了。在很多物理仿真的求解计算中，最常见的问题便是求解线性方程组，但是它们的维度一般都成千上万。但同时也存在一个比较有用的特性，即稀疏特性，现在有很多先进的算法来求解这些问题。这里我们介绍一种特殊情况下，即对角占优下三对角线性方程组的快速解法——Thomas分解或者追赶法，它的复杂度可以与N同阶，为 O ( 3 N − 2 ) O(3N-2) O(3N−2)。
-
对角占优的三对角方程组形如下：
[ b 1 c 1 0 ⋯ 0 a 2 b 2 c 2 ⋯ 0 0 ⋱ ⋱ ⋱ 0 0 0 a n − 1 b n − 1 c n − 1 0 0 0 a n b n ] [ x 1 x 2 ⋮ x n ] = [ f 1 f 2 ⋮ f n ] \begin{bmatrix} b_{1} & c_{1} & 0 & \cdots & 0 \ a_{2} & b_{2} & c_{2} & \cdots & 0\ 0 & \ddots & \ddots &\ddots&0\ 0 & 0 & a_{n-1} & b_{n-1} & c_{n-1}\ 0& 0& 0& a_{n} & b_{n} \end{bmatrix} \begin{bmatrix} x_1\ x_2\ \vdots\ x_{n} \end{bmatrix} = \begin{bmatrix} f_1\ f_2\ \vdots\ f_{n} \end{bmatrix} ⎣⎢⎢⎢⎢⎡b1a2000c1b2⋱000c2⋱an−10⋯⋯⋱bn−1an000cn−1bn⎦⎥⎥⎥⎥⎤⎣⎢⎢⎢⎡x1x2⋮xn⎦⎥⎥⎥⎤=⎣⎢⎢⎢⎡f1f2⋮fn⎦⎥⎥⎥⎤
-
对角占优条件：
- ∣ b 1 ∣ > ∣ c 1 ∣ > 0 |b_1| > |c_1| >0 ∣b1∣>∣c1∣>0
- ∣ b i ∣ ≥ ∣ a i ∣ + ∣ c i ∣ , a i , c i ≠ 0 , ( i = 2 , 3 , . . . , n − 1 ) |b_i| \geq |a_i| +|c_i|, a_i,c_i \neq0,(i=2,3,...,n-1) ∣bi∣≥∣ai∣+∣ci∣,ai,ci=0,(i=2,3,...,n−1)
- ∣ b n ∣ > ∣ a n ∣ > 0 |b_n| > |a_n|>0 ∣bn∣>∣an∣>0
-
类似于LU分解解线性方程组的做法，尝试将三对角矩阵分解为上下三角矩阵，形如：
A = [ b 1 c 1 0 ⋯ 0 a 2 b 2 c 2 ⋯ 0 0 ⋱ ⋱ ⋱ 0 0 0 a n − 1 b n − 1 c n − 1 0 0 0 a n b n ] = L U = [ α 1 0 0 ⋯ 0 γ 2 α 2 0 ⋯ 0 0 ⋱ ⋱ ⋱ 0 0 0 γ n − 1 α n − 1 0 0 0 0 γ n α n ] [ 1 β 1 0 ⋯ 0 0 1 β 2 ⋯ 0 0 ⋱ ⋱ ⋱ 0 0 0 0 1 β n − 1 0 0 0 0 1 ] A=\begin{bmatrix} b_{1} & c_{1} & 0 & \cdots & 0 \ a_{2} & b_{2} & c_{2} & \cdots & 0\ 0 & \ddots & \ddots &\ddots&0\ 0 & 0 & a_{n-1} & b_{n-1} & c_{n-1}\ 0& 0& 0& a_{n} & b_{n} \end{bmatrix}\ =LU=\begin{bmatrix} \alpha_{1} & 0 & 0 & \cdots & 0 \ \gamma_{2} & \alpha_{2} & 0 & \cdots & 0\ 0 & \ddots & \ddots &\ddots&0\ 0 & 0 & \gamma_{n-1} & \alpha_{n-1} & 0\ 0& 0& 0& \gamma_{n} & \alpha_{n} \end{bmatrix} \begin{bmatrix} 1 & \beta_{1} & 0 & \cdots & 0 \ 0 & 1 & \beta_{2} & \cdots & 0\ 0 & \ddots & \ddots &\ddots&0\ 0 & 0 & 0 &1 & \beta_{n-1}\ 0& 0& 0& 0 & 1 \end{bmatrix} A=⎣⎢⎢⎢⎢⎡b1a2000c1b2⋱000c2⋱an−10⋯⋯⋱bn−1an000cn−1bn⎦⎥⎥⎥⎥⎤=LU=⎣⎢⎢⎢⎢⎡α1γ20000α2⋱0000⋱γn−10⋯⋯⋱αn−1γn0000αn⎦⎥⎥⎥⎥⎤⎣⎢⎢⎢⎢⎡10000β11⋱000β2⋱00⋯⋯⋱10000βn−11⎦⎥⎥⎥⎥⎤ -
比较待定系数，可以逐步求解出各系数，（推导过程详见教材第七章P184）此时 A x = f Ax=f Ax=f等价于两个三角方程组 L y = f Ly=f Ly=f和 U x = y Ux=y Ux=y，得到追赶法求解过程：
- 计算 β i \beta_i βi： β 1 = c 1 / b 1 , β i = c i / ( b i − a i β i − 1 ) , i = ( 2 , 3 , . . . , n − 1 ) \beta_1=c_1/b_1,\beta_i=c_i/(b_i-a_i \beta_{i-1}), i=(2,3,...,n-1) β1=c1/b1,βi=ci/(bi−aiβi−1),i=(2,3,...,n−1)
- 解 L y = f Ly=f Ly=f： y 1 = f 1 / b 1 , y i = ( f i − a i y i − 1 ) / ( b i − a i β i − 1 ) , i = ( 2 , 3 , . . . , n ) y_1=f_1/b_1,y_i=(f_i-a_i y_{i-1})/(b_i-a_i \beta_{i-1}), i=(2,3,...,n) y1=f1/b1,yi=(fi−aiyi−1)/(bi−aiβi−1),i=(2,3,...,n)
- 解 U x = y Ux=y Ux=y： x n = y n , x i = y i − β i x i + 1 , ( i = n − 1 , n − 2 , . . . , 2 , 1 ) x_n=y_n,x_i=y_i-\beta_i x_{i+1}, (i=n-1,n-2,...,2,1) xn=yn,xi=yi−βixi+1,(i=n−1,n−2,...,2,1)
- 实际计算时只需要传入两个次对角线和主对角线的三个向量值便可以求解，内存消耗更小！
def check_SDD_vec(a, b, c, eps=1e-4, test=False):
"""check_SDD_vec 向量形式的严格对角占优条件判断
Args:
a ([np.darray]): [lower-secondary diagonal]
b ([np.darray]): [diagonal]
c ([np.darray]): [upper-secondary diagonal]
eps ([type], optional): [threshold]. Defaults to 1e-4.
test (bool, optional): [show testing information]. Defaults to False.
Returns:
[bool]: [whether or not]
Author: Junno
Date: 2022-04-28
"""
assert len(a) == len(c) and len(a)+1 == len(b)
N = len(b)
a_ = np.concatenate(([eps], a), axis=0)
c_ = np.concatenate((c, [eps]), axis=0)
if test:
print(a_, c_)
sum_a_c = np.abs(a_)+np.abs(c_)
b_ = np.abs(b)
for i in range(N):
if b_[i] < sum_a_c[i]:
return False
return True
def check_SDD_M(arr,eps=1e-4):
"""check_SDDM 判断输入矩阵是否满足严格对角占优条件
Args:
arr ([np.darray]): [input matrix]
eps ([type], optional): [threshold]. Defaults to 1e-4.
Returns:
[bool]: [whether or not]
Author: Junno
Date: 2022-04-28
"""
assert arr.shape[0]==arr.shape[1]
return check_SDD_vec(np.diagonal(arr,offset=-1),np.diagonal(arr,offset=0),np.diagonal(arr,offset=1),eps=eps)
def Thomas_Linear_Equation_Solve(a, b, c, d, eps=1e-4, test=False):
"""Thomas_Linear_Equation_Solve
Args:
a ([list or np.darray]): [lower-secondary diagonal]
b ([list or np.darray]): [diagonal]
c ([list or np.darray]): [upper-secondary diagonal]
d ([list or np.darray]): [b of Ax=b]
eps ([type], optional): [threshold]. Defaults to 1e-4.
test (bool, optional): [show testing information]. Defaults to False.
Returns:
x [list]: [solution of Ax=b]
Author: Junno
Date: 2022-04-28
"""
assert len(a) == len(c) and len(a)+1 == len(b)
if not check_SDD_vec(a, b, c, eps=eps):
print('The answer may be wrong because the input are not meet SDD condition completely')
N = len(b)
if test:
print(a, b, c, d)
beta = np.zeros_like(b)
x = np.zeros_like(b)
y = np.zeros_like(b)
beta[0] = c[0]/b[0]
y[0] = d[0]/b[0]
for i in range(1, N):
if i < N-1:
beta[i] = c[i]/(b[i]-a[i-1]*beta[i-1])
y[i] = (d[i]-a[i-1]*y[i-1])/(b[i]-a[i-1]*beta[i-1])
for i in range(N-1, -1, -1):
if i == N-1:
x[i] = y[i]
else:
x[i] = y[i]-beta[i]*x[i+1]
return x
A=np.array([[3,1,0,0],[1,2,1,0],[0,1,3,2],[0,0,3,4]],dtype=np.float32).reshape((4,4))
A
>>>
array([[3., 1., 0., 0.],
[1., 2., 1., 0.],
[0., 1., 3., 2.],
[0., 0., 3., 4.]], dtype=float32)
# create testing example
x=np.random.randn(4,1)
print(x)
>>>
[[-0.35812746]
[ 0.44621921]
[-0.45669147]
[ 1.85002579]]
f=A@x
print(f)
>>>
[[-0.62816316]
[ 0.07761949]
[ 2.77619638]
[ 6.03002876]]
Thomas_Linear_Equation_Solve(np.diagonal(A,offset=-1),np.diagonal(A,offset=0),np.diagonal(A,offset=1),f,test=False)
>>>
array([-0.35812744, 0.44621915, -0.45669138, 1.8500258 ], dtype=float32)
迭代法解线性方程组
-
对于大型的稀疏矩阵，有时候迭代求解方法在效率和成本上较直接求解更优。对于方程组 A x = b Ax=b Ax=b，可以通过变形得到 x = B x + f x=B x+f x=Bx+f，再通过迭代方法，逐步逼近真实解， x k + 1 = B x k + f x^{k+1}=B x^k +f xk+1=Bxk+f。其中，B称为迭代传递矩阵。
-
如果 lim k → ∞ x ( k ) \lim_{k \rightarrow \infty}x^{(k)} limk→∞x(k)存在，称迭代过程收敛。引进误差向量， ϵ k + 1 = x k + 1 − x ∗ \epsilon^{k+1}=x^{k+1}-x^* ϵk+1=xk+1−x∗，递推可以得到， ϵ k = B x k − 1 = B k ϵ 0 \epsilon^{k}=B x^{k-1}=B^k \epsilon^0 ϵk=Bxk−1=Bkϵ0。即若想要迭代收敛，需要知道B满足什么情况时有 B k → 0 ( k → ∞ ) B^k \rightarrow 0 \,\, (k \rightarrow \infty) Bk→0(k→∞)。
-
迭代法基本定理：当 ρ ( B ) < 1 \rho(B) <1 ρ(B)<1时，上述迭代法在任意初始向量 x 0 x_0 x0及任意f开始迭代后均收敛。其中， ρ ( B ) \rho(B) ρ(B)表示B的谱半径。可以由Jordan分解简单证明得到，可以参考教材P206-207。
-
由Jordan的证明过程还可以得到，当 ρ ( B ) \rho(B) ρ(B)越小时，迭代收敛的速度越快，若给定求解精度为 1 0 − s 10^{-s} 10−s，则至少需要的迭代次数可以由以下式子确定：
{ [ ρ ( B ) ] k ≤ 1 0 − s k ≥ s ln 10 − ln ρ ( B ) \begin{cases} [\rho(B)]^k \leq 10^{-s}\ k \geq \frac{s \ln{10}}{-\ln{\rho(B)}}\ \end{cases} {[ρ(B)]k≤10−sk≥−lnρ(B)sln10 -
由于任意阶矩阵范数都是谱半径的上届，所以可以有一些容易计算的矩阵范数得到收敛性的大概估计。设某一矩阵范数表示为 ∣ ∣ B ∣ ∣ v = q < 1 ||B||_v=q <1 ∣∣B∣∣v=q<1，第k步的误差估计可以表示为：
∣ ∣ x ∗ − x k ∣ ∣ v ≤ q 1 − q ∣ ∣ x k − x k − 1 ∣ ∣ v ||x^*-x^k||_v \leq \frac{q}{1-q} ||x^k-x^{k-1}||_v ∣∣x∗−xk∣∣v≤1−qq∣∣xk−xk−1∣∣v
- 设有方程组
∑
j
=
1
n
a
i
j
x
j
=
b
i
\sum_{j=1}^n a_{ij}x_j=b_i
∑j=1naijxj=bi，记为
A
x
=
b
Ax=b
Ax=b，将A分解为
A
=
D
+
L
+
U
A=D+L+U
A=D+L+U，格式如下：
D = [ a 11 a 22 ⋱ a n n ] D=\begin{bmatrix} a_{11} & & & \ & a_{22} & & \ & &\ddots& \ & & &a_{nn} \ \end{bmatrix} D=⎣⎢⎢⎡a11a22⋱ann⎦⎥⎥⎤
L = [ 0 a 21 0 ⋮ ⋱ ⋱ a n 1 ⋯ a n , n − 1 0 ] L=\begin{bmatrix} 0 & & & \ a_{21} & 0 & & \ \vdots & \ddots &\ddots& \ a_{n1} & \cdots & a_{n,n-1} &0 \ \end{bmatrix} L=⎣⎢⎢⎢⎡0a21⋮an10⋱⋯⋱an,n−10⎦⎥⎥⎥⎤
U = [ 0 a 12 ⋯ a 1 n 0 ⋱ ⋮ ⋱ a n − 1 , n 0 ] U=\begin{bmatrix} 0 & a_{12} & \cdots & a_{1n} \ & 0& \ddots & \vdots \ & &\ddots&a_{n-1,n} \ & & &0 \ \end{bmatrix} U=⎣⎢⎢⎢⎡0a120⋯⋱⋱a1n⋮an−1,n0⎦⎥⎥⎥⎤
- 对于方程组
∑
j
=
1
n
a
i
j
x
j
=
b
i
\sum_{j=1}^n a_{ij}x_j=b_i
∑j=1naijxj=bi，可以改写为
x
i
=
1
a
i
i
(
b
i
−
∑
j
=
1
,
j
≠
i
n
a
i
j
x
j
)
,
(
i
=
1
,
2
,
.
.
.
,
n
)
x_i=\frac{1}{a_{ii}}(b_i-\sum_{j=1,j \neq i}^n a_{ij}x_j),(i=1,2,...,n)
xi=aii1(bi−j=1,j=i∑naijxj),(i=1,2,...,n)若是以
x
=
B
x
+
f
x=Bx+f
x=Bx+f，则可以对照写出矩阵形式：
{ B = − D − 1 ( L + U ) f = D − 1 b \begin{cases} B=-D^{-1}(L+U)\ f=D^{-1}b\ \end{cases} {B=−D−1(L+U)f=D−1b
def Jacobi_iteration(A, b, N=0, eps=1e-4, x_0=None, norm='infty', own_method=False, test=False):
"""Jacobi_iteration
Args:
A ([np.darray]): [A]
b ([np.darray]): [b]
N (int, optional): [iteration time, it will give a suitable one when N is 0]. Defaults to 0.
eps ([float], optional): [error threshold]. Defaults to 1e-4.
x_0 ([np.darray], optional): [initial x for iteration]. Defaults to None.
norm (str, optional): [norm method to calculate fitting error,'infty','L1']. Defaults to 'infty'.
own_method (bool, optional): [use own method to calculate inv or norm, etc]. Defaults to False.
test (bool, optional): [show testing information]. Defaults to False.
Raises:
ValueError: [The spetrum of Iteration Matrix must be smaller than 1]
Returns:
[x]: [solution of Ax=b]
Author: Junno
Date: 2022-04-29
"""
assert A.shape[0] == b.shape[0]
# use np method or your own
Inv = np.linalg.inv if not own_method else Matrix_Solutions.Primary_Row_Transformation_Inv
SVD = np.linalg.svd if not own_method else Matrix_Solutions.svd
# substract D,L,U from A
D, L, U = np.diag(np.diag(A)), np.tril(A, k=-1), np.triu(A, k=1)
# calculate B and f
B = -Inv(D)@(L+U)
f = Inv(D)@b
# initial x for iteration
if x_0 == None:
x_0 = np.zeros_like(b)
# get eig of B
u, s, v = SVD(B)
# spectrum of B
eig_B = np.diag(s)
max_eig_B = np.max(eig_B)
# check convergence condition
if max_eig_B > 1:
print('The spetrum of Iteration Matrix must be smaller than 1')
raise ValueError
# matrix norm
if norm == 'infty':
F = Matrix_Solutions.m_infinite_norm
elif norm == 'L1':
F = Matrix_Solutions.m_L1_norm
q = F(B)
if test:
print('q =', q)
# iteration times
N = int(np.ceil(np.log10(1/eps)*np.log(10) /
(-np.log(max_eig_B)))) if N == 0 else N
# iteration step
while N:
x = B@x_0+f
if q < 1 and q/(1-q)*F(x-x_0) < eps:
break
x_0 = x
N -= 1
if test:
print("The speed of iteration is:{:.2f}".format(-np.log(max_eig_B)))
return x
A=np.array([[-1.35 , 0.209, -0.218, -0.123],
[ 0.179, 1.616, -0.104, -0.497],
[ 0.333, 0.278, 2.417, -0.226],
[ 1.128, 0.346, 0.104, 1.237]],dtype=np.float32).reshape((4,4))
x=np.random.randn(A.shape[0],1)
print(x)
>>>
[[-1.872]
[ 0.693]
[ 0.809]
[-1.577]]
b=A@x
print(b)
>>>
[[ 2.691]
[ 1.485]
[ 1.881]
[-3.738]]
ans=Jacobi_iteration(A, b, N=20, norm='infty', own_method=False, test=False)
print(ans)
>>>
[[-1.872]
[ 0.693]
[ 0.809]
[-1.577]]
Matrix_Solutions.m_L1_norm(ans-x)
>>>
4.621872773391544e-09
Gauss-Seidel迭代法
- Jacobi迭代的第k步可以表示为：
x
i
k
+
1
=
1
a
i
i
(
b
i
−
∑
j
=
1
,
j
≠
i
n
a
i
j
x
j
k
)
x_i^{k+1}=\frac{1}{a_{ii}}(b_i-\sum_{j=1,j \neq i}^n a_{ij}x_j^{k})
xik+1=aii1(bi−j=1,j=i∑naijxjk) 在计算k+1步中，要用到所有第k步的结果，但是在计算第i个x的第k+1步之前，新计算出的
x
1
k
+
1
,
x
1
k
+
1
,
.
.
.
,
x
i
−
1
k
+
1
x_1^{k+1},x_1^{k+1},...,x_{i-1}^{k+1}
x1k+1,x1k+1,...,xi−1k+1没有被用到。 Gauss-Seidel迭代法便是用上这些新计算出的结果来参与计算。表示如下：
x i k + 1 = 1 a i i ( b i − ∑ j = 1 i − 1 a i j x j k + 1 − ∑ j = i + 1 n a i j x j k ) x_i^{k+1}=\frac{1}{a_{ii}}(b_i-\sum_{j=1}^{i-1} a_{ij}x_j^{k+1}-\sum_{j=i+1}^n a_{ij}x_j^{k}) xik+1=aii1(bi−j=1∑i−1aijxjk+1−j=i+1∑naijxjk)
写成矩阵形式，即：
D x k + 1 = b − L x k + 1 − U x k ( D + L ) x k + 1 = b − U x k x k + 1 = − ( D + L ) − 1 U x k + ( D + L ) − 1 b Dx^{k+1}=b-Lx^{k+1}-Ux^{k}\ (D+L)x^{k+1}=b-Ux^{k}\ x^{k+1}=-(D+L)^{-1}Ux^k+(D+L)^{-1}b Dxk+1=b−Lxk+1−Uxk(D+L)xk+1=b−Uxkxk+1=−(D+L)−1Uxk+(D+L)−1b
按照 x = B x + f x=Bx+f x=Bx+f的格式，即有：
{ B = − ( D + L ) − 1 U f = ( D + L ) − 1 b \begin{cases} B=-(D+L)^{-1}U\ f=(D+L)^{-1}b\ \end{cases} {B=−(D+L)−1Uf=(D+L)−1b
# 整合Jacobi和Gauss_Seidel到一起
def Serial_iteration(A, b, N=0, eps=1e-4, x_0=None, iter_way='Gauss_Seidel', norm='infty', own_method=False, test=False):
"""Serial_iteration:Jacobi and Gauss_Seidel
Args:
A ([np.darray]): [A]
b ([np.darray]): [b]
N (int, optional): [iteration time, it will give a suitable one when N is 0]. Defaults to 0.
eps ([float], optional): [error threshold]. Defaults to 1e-4.
x_0 ([np.darray], optional): [initial x for iteration]. Defaults to None.
iter_way (str, optional): [iteration method,'Gauss_Seidel','Jacobi']. Defaults to 'Gauss_Seidel'.
norm (str, optional): [norm method to calculate fitting error,'infty','L1']. Defaults to 'infty'.
own_method (bool, optional): [use own method to calculate inv or norm, etc]. Defaults to False.
test (bool, optional): [show testing information]. Defaults to False.
Raises:
ValueError: [The spetrum of Iteration Matrix must be smaller than 1]
Returns:
[x]: [solution of Ax=b]
Author: Junno
Date: 2022-04-29
"""
assert A.shape[0] == b.shape[0]
# use np method or your own
Inv = np.linalg.inv if not own_method else Matrix_Solutions.Primary_Row_Transformation_Inv
SVD = np.linalg.svd if not own_method else Matrix_Solutions.svd
# substract D,L,U from A
D, L, U = np.diag(np.diag(A)), np.tril(A, k=-1), np.triu(A, k=1)
# calculate B and f
if iter_way == 'Jacobi':
B = -Inv(D)@(L+U)
f = Inv(D)@b
elif iter_way == 'Gauss_Seidel':
B = -Inv(D+L)@U
f = Inv(D+L)@b
# initial x for iteration
if x_0 == None:
x_0 = np.zeros_like(b)
# get eig of B
u, s, v = SVD(B)
# spectrum of B
eig_B = np.diag(s)
max_eig_B = np.max(eig_B)
# check convergence condition
if max_eig_B > 1:
print('The spetrum of Iteration Matrix must be smaller than 1')
raise ValueError
# matrix norm
if norm == 'infty':
F = Matrix_Solutions.m_infinite_norm
elif norm == 'L1':
F = Matrix_Solutions.m_L1_norm
q = F(B)
if test:
print('q =', q)
# iteration times
N = int(np.ceil(np.log10(1/eps)*np.log(10) /
(-np.log(max_eig_B)))) if N == 0 else N
# iteration step
while N:
x = B@x_0+f
if q < 1 and q/(1-q)*F(x-x_0) < eps:
print('last iteration error: ',q/(1-q)*F(x-x_0))
break
x_0 = x
N -= 1
if test:
print("The speed of iteration is: {:.2f}".format(-np.log(max_eig_B)))
return x
# use example above
ans=Serial_iteration(A, b, N=0, eps=1e-4, iter_way='Gauss_Seidel',norm='infty', own_method=False, test=True)
print(ans)
>>>
q = 0.41703337
last iteration error: 4.200151608806934e-05
The speed of iteration is: 0.96
[[-1.873]
[ 0.694]
[ 0.809]
[-1.576]]
-
对于Gauss-Seidel迭代法， x i k + 1 = 1 a i i ( b i − ∑ j = 1 i − 1 a i j x j k + 1 − ∑ j = i + 1 n a i j x j k ) x_i^{k+1}=\frac{1}{a_{ii}}(b_i-\sum_{j=1}^{i-1} a_{ij}x_j^{k+1}-\sum_{j=i+1}^n a_{ij}x_j^{k}) xik+1=aii1(bi−j=1∑i−1aijxjk+1−j=i+1∑naijxjk)
引入 x i k x_i^k xik，有
{ x i k + 1 = x i k + Δ x i Δ x i = 1 a i i ( b i − ∑ j = 1 i − 1 a i j x j k − ∑ j = i n a i j x j k ) \begin{cases} x_i^{k+1}=x_i^k+\Delta x_i\ \Delta x_i=\frac{1}{a_{ii}}(b_i-\sum_{j=1}^{i-1} a_{ij}x_j^{k}-\sum_{j=i}^n a_{ij}x_j^{k})\ \end{cases} {xik+1=xik+ΔxiΔxi=aii1(bi−∑j=1i−1aijxjk−∑j=inaijxjk)
引入松弛因子 ω \omega ω， x i k + 1 = x i k + ω a i i ( b i − ∑ j = 1 i − 1 a i j x j k − ∑ j = i n a i j x j k ) \large x_i^{k+1}=x_i^k+\frac{\omega}{a_{ii}}(b_i-\sum_{j=1}^{i-1} a_{ij}x_j^{k}-\sum_{j=i}^n a_{ij}x_j^{k}) xik+1=xik+aiiω(bi−∑j=1i−1aijxjk−∑j=inaijxjk)，选取合适的松弛因子就可以加速收敛过程，上式即为SOR的表达公式。 -
推到SOR的矩阵形式：
- a i i x i k + 1 = ( 1 − ω ) a i i x i k + ω ( b i − ∑ j = 1 i − 1 a i j x j k − ∑ j = i + 1 n a i j x j k ) a_{ii}x_i^{k+1}=(1-\omega)a_{ii}x_i^{k}+\omega(b_i-\sum_{j=1}^{i-1} a_{ij}x_j^{k}-\sum_{j=i+1}^n a_{ij}x_j^{k}) aiixik+1=(1−ω)aiixik+ω(bi−j=1∑i−1aijxjk−j=i+1∑naijxjk)
- D x k + 1 = ( 1 − ω ) D x k + ω ( b − L x k + 1 − U x k ) Dx^{k+1}=(1-\omega)Dx^{k}+\omega(b-Lx^{k+1}-Ux^{k}) Dxk+1=(1−ω)Dxk+ω(b−Lxk+1−Uxk)
- x k + 1 = ( D + ω L ) − 1 [ ( 1 − ω ) D − ω U ] x k + w ( D + ω L ) − 1 b x^{k+1}=(D+\omega L)^{-1}[(1-\omega)D-\omega U]x^k+w(D+\omega L)^{-1}b xk+1=(D+ωL)−1[(1−ω)D−ωU]xk+w(D+ωL)−1b
-
收敛的充要条件：迭代矩阵 B = ( D + ω L ) − 1 [ ( 1 − ω ) D − ω U ] B=(D+\omega L)^{-1}[(1-\omega)D-\omega U] B=(D+ωL)−1[(1−ω)D−ωU]的谱半径<1。
-
充分条件：若SOR收敛，则 0 < ω < 2 0<\omega <2 0<ω<2，证明如下：
- SOR方法收敛，则有 ρ ( B ) < 1 \rho(B)<1 ρ(B)<1，又 ∣ d e t ( B ) ∣ = ∣ λ 1 λ 2 ⋯ λ n ∣ |det(B)|=|\lambda_1 \lambda_2 \cdots \lambda_n| ∣det(B)∣=∣λ1λ2⋯λn∣，即 ∣ d e t ( B ) ∣ < 1 |det(B)|<1 ∣det(B)∣<1
- d e t ( ( D + ω L ) − 1 ) = 1 / d e t ( ( D + ω L ) ) = ∏ i = 0 n 1 / a i i det((D+\omega L)^{-1})=1/det((D+\omega L))=\prod_{i=0}^n 1/a_{ii} det((D+ωL)−1)=1/det((D+ωL))=∏i=0n1/aii
- d e t ( ( 1 − ω ) D − ω U ) = ( 1 − ω ) n ∏ i = 0 n 1 / a i i det((1-\omega)D-\omega U)=(1-\omega)^n \prod_{i=0}^n 1/a_{ii} det((1−ω)D−ωU)=(1−ω)n∏i=0n1/aii
- d e t ( B ) = ( 1 − ω ) n det(B)=(1-\omega)^n det(B)=(1−ω)n， ( 1 − ω ) n < 1 (1-\omega)^n<1 (1−ω)n<1，
- 即 0 < ω < 2 0< \omega <2 0<ω<2
- 依照前面的代码，很容易写出SOR：
def SOR_iteration(A, b, eps=1e-4, x_0=None, omega=1, norm='infty', own_method=False, test=False):
"""Serial_iteration:Jacobi and Gauss_Seidel
Args:
A ([np.darray]): [A]
b ([np.darray]): [b]
eps ([float], optional): [error threshold]. Defaults to 1e-4.
x_0 ([np.darray], optional): [initial x for iteration]. Defaults to None.
omega (float, optional): [iteration coefficient]. Defaults to 1 means to Gauss_Seidel method.
norm (str, optional): [norm method to calculate fitting error,'infty','L1']. Defaults to 'infty'.
own_method (bool, optional): [use own method to calculate inv or norm, etc]. Defaults to False.
test (bool, optional): [show testing information]. Defaults to False.
Raises:
ValueError: [The spetrum of Iteration Matrix must be smaller than 1]
Returns:
[x]: [solution of Ax=b]
Author: Junno
Date: 2022-04-29
"""
assert A.shape[0] == b.shape[0]
assert 0<omega<2,'iteration coefficient must be during (0,2) to meet the convergence condition'
# use np method or your own
Inv = np.linalg.inv if not own_method else Matrix_Solutions.Primary_Row_Transformation_Inv
SVD = np.linalg.svd if not own_method else Matrix_Solutions.svd
# substract D,L,U from A
D, L, U = np.diag(np.diag(A)), np.tril(A, k=-1), np.triu(A, k=1)
# calculate B and f
D_ = Inv((D+omega*L))
B = D_@((1-omega)*D-omega*U)
f = omega*D_@b
# initial x for iteration
if x_0 == None:
x_0 = np.zeros_like(b)
# get eig of B
u, s, v = SVD(B)
# spectrum of B
max_eig_B = np.max(np.diag(s))
# check convergence condition
if max_eig_B > 1:
print('The spetrum of Iteration Matrix must be smaller than 1')
raise ValueError
# matrix norm
if norm == 'infty':
F = Matrix_Solutions.m_infinite_norm
elif norm == 'L1':
F = Matrix_Solutions.m_L1_norm
q = F(B)
# iteration times
least_steps = int(np.ceil(np.log10(1/eps)*np.log(10) /
(-np.log(max_eig_B))))
N=0
# iteration step
while N<least_steps:
x = B@x_0+f
if q < 1 and q/(1-q)*F(x-x_0) < eps:
print('last iteration error: ',q/(1-q)*F(x-x_0))
break
x_0 = x
N += 1
if test:
print('q =', q)
print("Iteration times:{}".format(N))
print("The speed of iteration is:{:.2f}".format(-np.log(max_eig_B)))
return x
# use example above
ans=SOR_iteration(A, b, eps=1e-4, omega=0.95,norm='infty', own_method=False, test=True)
print(ans)
>>>
q = 0.437037
last iteration error: 5.658720150643954e-05
Iteration times:6
The speed of iteration is:0.99
[[-0.709]
[-0.886]
[ 1.286]
[ 0.215]]
ans-x
>>>
array([[-0.],
[ 0.],
[ 0.],
[ 0.]])
A=array([[ 0.758, 0.28 , 0.383, 0.259],
[-0.165, 2.509, -0.542, 0.271],
[ 0.762, 1.151, 1.999, 0.842],
[-0.724, 0.515, -0.081, -1.539]],dtype=np.float32).reshape((4,4))
x=np.random.randn(A.shape[0],1)
print(x)
>>>
[[ 0.178]
[-0.471]
[-0.099]
[-0.02 ]]
b=A@x
print(b)
>>>
[[-0.04 ]
[-1.162]
[-0.622]
[-0.332]]
- 对于上面这个方程组， ω \omega ω的选择所导致的满足相邻误差条件的最少迭代数如下：
def Serial_iteration(A, b, omega=0, eps=1e-4, x_0=None, iter_way='Gauss_Seidel', norm='infty', own_method=False, test=False):
"""Serial_iteration:Jacobi and Gauss_Seidel
Args:
A ([np.darray]): [A]
b ([np.darray]): [b]
omega ([float], optional): [iteration coefficient for SOR method]. Defaults to 0.
eps ([float], optional): [error threshold]. Defaults to 1e-4.
x_0 ([np.darray], optional): [initial x for iteration]. Defaults to None.
iter_way (str, optional): [iteration method,'Gauss_Seidel','Jacobi','SOR']. Defaults to 'Gauss_Seidel'.
norm (str, optional): [norm method to calculate fitting error,'infty','L1']. Defaults to 'infty'.
own_method (bool, optional): [use own method to calculate inv or norm, etc]. Defaults to False.
test (bool, optional): [show testing information]. Defaults to False.
Raises:
ValueError: [The spetrum of Iteration Matrix must be smaller than 1]
Returns:
[x]: [solution of Ax=b]
Author: Junno
Date: 2022-04-29
"""
from Matrix_Solutions import Matrix_Solutions
assert A.shape[0] == b.shape[0]
# use np method or your own
Inv = np.linalg.inv if not own_method else Matrix_Solutions.Primary_Row_Transformation_Inv
SVD = np.linalg.svd if not own_method else Matrix_Solutions.svd
# substract D,L,U from A
D, L, U = np.diag(np.diag(A)), np.tril(A, k=-1), np.triu(A, k=1)
# calculate B and f
if iter_way == 'Jacobi':
B = -Inv(D)@(L+U)
f = Inv(D)@b
elif iter_way == 'Gauss_Seidel':
B = -Inv(D+L)@U
f = Inv(D+L)@b
elif iter_way == 'SOR':
assert 0 < omega < 2, 'iteration coefficient must be during (0,2) to meet the convergence condition'
D_ = Inv((D+omega*L))
B = D_@((1-omega)*D-omega*U)
f = omega*D_@b
# initial x for iteration
if x_0 == None:
x_0 = np.zeros_like(b)
# get eig of B
u, s, v = SVD(B)
# spectrum of B
eig_B = np.diag(s)
max_eig_B = np.max(eig_B)
# check convergence condition
if max_eig_B > 1:
print('The spetrum of Iteration Matrix must be smaller than 1')
raise ValueError
# matrix norm
if norm == 'infty':
F = Matrix_Solutions.m_infinite_norm
elif norm == 'L1':
F = Matrix_Solutions.m_L1_norm
q = F(B)
# iteration times
least_steps = int(np.ceil(np.log10(1/eps)*np.log(10) /
(-np.log(max_eig_B))))
N = 0
# iteration step
while N < least_steps:
x = B@x_0+f
if q < 1 and q/(1-q)*F(x-x_0) < eps:
print('last iteration error: ', q/(1-q)*F(x-x_0))
break
x_0 = x
N += 1
if test:
print('q =', q)
print("Iteration times:{}".format(N))
print(
"The speed of iteration is:{:.2f}".format(-np.log(max_eig_B)))
return x
- 有问题的uu请在留言板提出指正，关于代码中Matrix_Solutions的部分是在矩阵论部分写的常用函数集合，若需请挪步~
欢迎分享，转载请注明来源：内存溢出
微信扫一扫
支付宝扫一扫
评论列表（0条）