# Start writing code here...
实验一
实验二 线性方程组的直接解法
高斯消去法
import numpy as np
#消元过程
def gauss_elimination(A,b):#输入系数矩阵和b
n=A.shape[0]
for i in range(n-1):
for j in range(i+1,n):
b[j]=b[j]-b[i]/A[i,i]*A[j,i] #!!!b的运算应在A之前，不然会导致减的全是0
A[j,...]= A[j,...]-A[i,...]/A[i,i]*A[j,i]
print('上三角矩阵\n',A)
return A,b
#回代过程
def gauss_back(A_up,b):#输入上三角矩阵和b
n=A_up.shape[0]
for i in range(n-1,0,-1):#自下而上
for j in range(i-1,-1,-1):
b[j]=b[j]-b[i]/A_up[i,i]*A[j,i]
A_up[j,...]= A_up[j,...]-A_up[i,...]/A_up[i,i]*A[j,i]
#print(b[j],'=',b[j],'-',b[i],'/',A_up[i,i],'*',A[j,i])
print('对角阵\n',A_up)
return A_up,b
A=np.array([(10,-7,0,1),
(-3 ,2.099999, 6, 2 ),\
(5 ,-1 ,5 ,-1), \
(2 ,1 ,0 ,2)])
b=np.array([8,5.900001,5,1]).T
#消元过程
A_up,b=gauss_elimination(A,b)
#回代过程
A_diag,b=gauss_back(A_up,b)
#求根
root=b/np.diagonal(A_diag)
root
上三角矩阵
[[ 1.00000000e+01 -7.00000000e+00 0.00000000e+00 1.00000000e+00]
[ 0.00000000e+00 -1.00000000e-06 6.00000000e+00 2.30000000e+00]
[ 0.00000000e+00 0.00000000e+00 1.50000050e+07 5.74999850e+06]
[ 0.00000000e+00 0.00000000e+00 0.00000000e+00 5.07999891e+00]]
对角阵
[[ 1.00000000e+01 0.00000000e+00 0.00000000e+00 0.00000000e+00]
[ 0.00000000e+00 -1.00000000e-06 0.00000000e+00 0.00000000e+00]
[ 0.00000000e+00 0.00000000e+00 1.50000050e+07 0.00000000e+00]
[ 0.00000000e+00 0.00000000e+00 0.00000000e+00 5.07999891e+00]]
列主元高斯消去法
import numpy as np
#消元过程
def gauss_elimination(A,b):#输入系数矩阵和b
n=A.shape[0]
for i in range(n-1):
# 在此增加换主元操作
fabs_col=np.fabs(A[i:n,i])#取该列
#print(fabs_col)
colMAX = max(fabs_col)
indexMAX = np.where(fabs_col == colMAX)[0][0] #找最大元素位置
#print(indexMAX)
#交换
A[[i,i+indexMAX],:]=A[[i+indexMAX,i],:]
b[i],b[i+indexMAX]=b[i+indexMAX],b[i]
for j in range(i+1,n):
b[j]=b[j]-b[i]/A[i,i]*A[j,i] #!!!b的运算应在A之前，不然会导致减的全是0
A[j,...]= A[j,...]-A[i,...]/A[i,i]*A[j,i]
print('上三角矩阵\n',A)
return A,b
#回代过程
def gauss_back(A_up,b):#输入上三角矩阵和b
n=A_up.shape[0]
for i in range(n-1,0,-1):#自下而上
for j in range(i-1,-1,-1):
b[j]=b[j]-b[i]/A_up[i,i]*A[j,i]
A_up[j,...]= A_up[j,...]-A_up[i,...]/A_up[i,i]*A[j,i]
#print(b[j],'=',b[j],'-',b[i],'/',A_up[i,i],'*',A[j,i])
print('对角阵\n',A_up)
return A_up,b
A=np.array([(10,-7,0,1),
(-3 ,2.099999, 6, 2 ),\
(5 ,-1 ,5 ,-1), \
(2 ,1 ,0 ,2)])
b=np.array([8,5.900001,5,1]).T
#消元过程
A_up,b=gauss_elimination(A,b)
#回代过程
A_diag,b=gauss_back(A_up,b)
#求根
root=b/np.diagonal(A_diag)
root
上三角矩阵
[[10. -7. 0. 1. ]
[ 0. 2.5 5. -1.5 ]
[ 0. 0. 6.000002 2.2999994 ]
[ 0. 0. 0. 5.07999891]]
对角阵
[[10. 0. 0. 0. ]
[ 0. 2.5 0. 0. ]
[ 0. 0. 6.000002 0. ]
[ 0. 0. 0. 5.07999891]]
直接三角分解法
import numpy as np
def get_L_U(A):
L,U = np.zeros([n,n]),np.zeros([n, n])
for i in range(n):
L[i][i]=1
if i==0:
U[0][0] = A[0][0]
for j in range(1,n):
U[0][j]=A[0][j]
L[j][0]=A[j][0]/U[0][0]
else:
for j in range(i, n):#U
temp=0
for k in range(0, i):
temp = temp+L[i][k] * U[k][j]
U[i][j]=A[i][j]-temp
for j in range(i+1, n):#L
temp = 0
for k in range(0, i ):
temp = temp + L[j][k] * U[k][i]
L[j][i] = (A[j][i] - temp)/U[i][i]
print('L:\n',L,'\nU:\n',U)
return L,U
def cal_y(L,b):
#y[0]=b[0]
y=b #直接全赋值得了
for k in range(0,n):
for j in range(0,k):
y[k]=y[k]-L[k,j]*y[j]
print('y:\n',y)
return y
def cal_x(U,y):
x=y/np.diagonal(U) #这里同样是全部赋值
for k in range(n-2,-1,-1):
for j in range(k+1,n): #将公式拆分
x[k]=x[k]-U[k,j]*x[j]/U[k,k]
print('x:\n',x)
return x
A=np.array([(10,-7,0,1),
(-3 ,2.099999, 6, 2 ),\
(5 ,-1 ,5 ,-1), \
(2 ,1 ,0 ,2)])
b=np.array([8,5.900001,5,1]).T
n=A.shape[0]
L,U=get_L_U(A)
y=cal_y(L,b)
x=cal_x(U,y)
L:
[[ 1.0000000e+00 0.0000000e+00 0.0000000e+00 0.0000000e+00]
[-3.0000000e-01 1.0000000e+00 0.0000000e+00 0.0000000e+00]
[ 5.0000000e-01 -2.5000000e+06 1.0000000e+00 0.0000000e+00]
[ 2.0000000e-01 -2.4000000e+06 9.5999968e-01 1.0000000e+00]]
U:
[[ 1.00000000e+01 -7.00000000e+00 0.00000000e+00 1.00000000e+00]
[ 0.00000000e+00 -1.00000000e-06 6.00000000e+00 2.30000000e+00]
[ 0.00000000e+00 0.00000000e+00 1.50000050e+07 5.74999850e+06]
[ 0.00000000e+00 0.00000000e+00 0.00000000e+00 5.07999891e+00]]
y:
[8.00000000e+00 8.30000100e+00 2.07500035e+07 5.07999891e+00]
x:
[-9.5955556e-10 -1.0000000e+00 1.0000000e+00 1.0000000e+00]