涉及Crout分解、追趕法的線性方程組求解方法的Python實現。
Codes
def CroutLU(A:np.ndarray)->Tuple[np.ndarray,np.ndarray]:
"""Crout LU分解演算法,A=LU
input:
A: (n,n) np.ndarray,方陣
output:
L: (n,n) np.ndarray,下三角矩陣
U: (n,n) np.ndarray,上三角矩陣,對角線元素為1.0
usage:
```python
A=np.array([[1,2,3,4],
[1,4,9,16],
[1,8,27,64],
[1,16,81,256]])
L,U=CroutLU(A)
print("L矩陣:\n", L)
print("U矩陣:\n", U)
# 驗證分解是否正確
print("驗證LU是否等於A:\n", np.dot(L, U))
輸出:
L矩陣:
[[ 1. 0. 0. 0.]
[ 1. 2. 0. 0.]
[ 1. 6. 6. 0.]
[ 1. 14. 36. 24.]]
U矩陣:
[[1. 2. 3. 4.]
[0. 1. 3. 6.]
[0. 0. 1. 4.]
[0. 0. 0. 1.]]
驗證LU是否等於A:
[[ 1. 2. 3. 4.]
[ 1. 4. 9. 16.]
[ 1. 8. 27. 64.]
[ 1. 16. 81. 256.]]
```
"""
row,col=A.shape
n=row
L=np.zeros((n,n))
U=np.zeros((n,n))
for k in range(n):
# 迴圈列,從k+1列到n列,i=k,...n-1
for i in range(k,n):
L[i,k]=A[i,k]-sum([L[i,s]*U[s,k] for s in range(k)])
for j in range(k-1,n):
U[k,j]=(A[k,j]-sum([L[k,s]*U[s,j] for s in range(k)]))/L[k,k]
return L,U
def LUChaseMethod(A:np.ndarray,d:np.ndarray)->np.ndarray:
"""LU分解法,追趕法求解線性方程組Ax=d
求解三對角矩陣A,d的線性方程組Ax=d,其中A為三對角矩陣,d為右端常數
"""
n=A.shape[0]
# x: x1,x2...xn
x=np.zeros(n)
a=np.zeros(n)
# a:a1,a2...an
# b:b1...bn
# c:c1...cn-1
a[1:],b,c=np.diag(A,k=-1).copy(),np.diag(A,k=0).copy(),np.diag(A,k=1).copy()
L,U=CroutLU(A)
# size: (n-1,) , u0,u1...u(n-1),u0=0
us=np.zeros(n)
us[1:]=np.diag(U,k=1).copy()
# size: (n,) ,ls : l1...ln
ls=np.diag(L,k=0).copy()
# size: (n-1,) , v2...vn-1
vs=np.diag(L,k=-1).copy()
# y: y0,y1...yn , y0=0
y=np.zeros(n+1)
# i=1....n-1
for i in range(1,n):
# print(f"第{i}次迭代")
y_i_1=y[i-1]
a_i=a[i-1]
b_i=b[i-1]
c_i=c[i-1]
d_i=d[i-1]
u_i_1=us[i-1]
l_i=b_i-a_i*u_i_1
u_i=c_i/l_i
y_i=(d_i-a_i*y_i_1)/l_i
y[i]=y_i
ls[i-1]=l_i
us[i]=u_i
ls[-1]=b[n-1]-a[n-1]*us[n-1]
y[n]=(d[n-1]-y[n-1]*a[n-1])/ls[n-1]
x[n-1]=y[n]
for i in range(n-2,-1,-1):
# xi=yi-ui*x(i+1),i=n-2...1
x[i]=y[i+1]-us[i+1]*x[i+1]
return x
Vaildation
from formu_lib import *
import numpy as np
import matplotlib.pyplot as plt
A=np.array([[1,2,3,4],
[1,4,9,16],
[1,8,27,64],
[1,16,81,256]])
L,U=CroutLU(A)
print("L矩陣:\n", L)
print("U矩陣:\n", U)
# 驗證分解是否正確
print("驗證LU是否等於A:\n", np.dot(L, U))
- 輸出:
L矩陣:
[[ 1. 0. 0. 0.]
[ 1. 2. 0. 0.]
[ 1. 6. 6. 0.]
[ 1. 14. 36. 24.]]
U矩陣:
[[1. 2. 3. 4.]
[0. 1. 3. 6.]
[0. 0. 1. 4.]
[0. 0. 0. 1.]]
驗證LU是否等於A:
[[ 1. 2. 3. 4.]
[ 1. 4. 9. 16.]
[ 1. 8. 27. 64.]
[ 1. 16. 81. 256.]]
from formu_lib import *
import numpy as np
import matplotlib.pyplot as plt
a=np.array([[2,-1,0,0],[-1,3,-2,0],[0,-2,4,-3],[0,0,-3,5]])
d=np.array([6,1,-2,1])
x=LUChaseMethod(a,d)
print(f"x={x}")
# x=[5. 4. 3. 2.]
答案:[5. 4. 3. 2.],ok