通过python判断一个矩阵是否是不相容的,并且再写一个backsubstitution函数,来解出矩阵.判断是否是不相容(inconsistent)需要用numpy.nonzero方法
import numpy as np
import warnings
def swapRows(A, i, j):
"""
interchange two rows of A
operates on A in place
"""
tmp = A[i].copy()
A[i] = A[j]
A[j] = tmp
def relError(a, b):
"""
compute the relative error of a and b
"""
with warnings.catch_warnings():
warnings.simplefilter("error")
try:
return np.abs(a-b)/np.max(np.abs(np.array([a, b])))
except:
return 0.0
def rowReduce(A, i, j, pivot):
"""
reduce row j using row i with pivot pivot, in matrix A
operates on A in place
"""
factor = A[j][pivot] / A[i][pivot]
for k in range(len(A[j])):
if np.isclose(A[j][k], factor * A[i][k]):
A[j][k] = 0.0
else:
A[j][k] = A[j][k] - factor * A[i][k]
# stage 1 (forward elimination)
def forwardElimination(B):
"""
Return the row echelon form of B
"""
A = B.copy().astype(float)
m, n = np.shape(A)
for i in range(m-1):
# Let lefmostNonZeroCol be the position of the leftmost nonzero value
# in row i or any row below it
leftmostNonZeroRow = m
leftmostNonZeroCol = n
## for each row below row i (including row i)
for h in range(i,m):
## search, starting from the left, for the first nonzero
for k in range(i,n):
if (A[h][k] != 0.0) and (k < leftmostNonZeroCol):
leftmostNonZeroRow = h
leftmostNonZeroCol = k
break
# if there is no such position, stop
if leftmostNonZeroRow == m:
break
# If the leftmostNonZeroCol in row i is zero, swap this row
# with a row below it
# to make that position nonzero. This creates a pivot in that position.
if (leftmostNonZeroRow > i):
swapRows(A, leftmostNonZeroRow, i)
# Use row reduction operations to create zeros in all positions
# below the pivot.
for h in range(i+1,m):
rowReduce(A, i, h, leftmostNonZeroCol)
return A
####################
# If any operation creates a row that is all zeros except the last element,
# the system is inconsistent; stop.
def inconsistentSystem(A):
# 判断一个矩阵是否是不相容的,只能用numpy的nonzero方法
def backsubstitution(B):
"""
return the reduced row echelon form matrix of B
"""
#####################
# 测试的值
a = np.array([[1, 2, 3,24], [5, 4, 6,24], [10, 9, 8,9]])
backsubstitution(a)
# 最后的结果应该是[1,0,0,-8],[0,1,0,1],[0,0,1,10]
只用写最后两个函数就好了,之前写完了
import numpy as np
import warnings
def swapRows(A, i, j):
"""
interchange two rows of A
operates on A in place
"""
tmp = A[i].copy()
A[i] = A[j]
A[j] = tmp
def relError(a, b):
"""
compute the relative error of a and b
"""
with warnings.catch_warnings():
warnings.simplefilter("error")
try:
return np.abs(a - b) / np.max(np.abs(np.array([a, b])))
except:
return 0.0
def rowReduce(A, i, j, pivot):
"""
reduce row j using row i with pivot pivot, in matrix A
operates on A in place
"""
factor = A[j][pivot] / A[i][pivot]
for k in range(len(A[j])):
if np.isclose(A[j][k], factor * A[i][k]):
A[j][k] = 0.0
else:
A[j][k] = A[j][k] - factor * A[i][k]
# stage 1 (forward elimination)
def forwardElimination(B):
"""
Return the row echelon form of B
"""
A = B.copy().astype(float)
m, n = np.shape(A)
for i in range(m - 1):
# Let lefmostNonZeroCol be the position of the leftmost nonzero value
# in row i or any row below it
leftmostNonZeroRow = m
leftmostNonZeroCol = n
## for each row below row i (including row i)
for h in range(i, m):
## search, starting from the left, for the first nonzero
for k in range(i, n):
if (A[h][k] != 0.0) and (k < leftmostNonZeroCol):
leftmostNonZeroRow = h
leftmostNonZeroCol = k
break
# if there is no such position, stop
if leftmostNonZeroRow == m:
break
# If the leftmostNonZeroCol in row i is zero, swap this row
# with a row below it
# to make that position nonzero. This creates a pivot in that position.
if (leftmostNonZeroRow > i):
swapRows(A, leftmostNonZeroRow, i)
# Use row reduction operations to create zeros in all positions
# below the pivot.
for h in range(i + 1, m):
rowReduce(A, i, h, leftmostNonZeroCol)
return A
####################
# If any operation creates a row that is all zeros except the last element,
# the system is inconsistent; stop.
def inconsistentSystem(A):
temp = forwardElimination(A)
return max(np.nonzero(temp[:, :-1])[0]) != max(np.nonzero(temp)[0])
# 判断一个矩阵是否是不相容的,只能用numpy的nonzero方法
##原文链接：https://blog.csdn.net/mr_jjpolarbear/article/details/88649587
def rsmat(arbmat):
""" Convert an arbitrary matrix to a simplest matrix """
arbmat = arbmat.astype(float)
row_number, column_number = arbmat.shape
if row_number == 1:
if arbmat[0, 0] != 0:
return (arbmat / arbmat[0, 0])
else:
return arbmat
else:
rc_number = min(row_number, column_number)
anarbmat = arbmat.copy()
r = 0
for n in range(rc_number):
s_row = -1
for i in arbmat[r:row_number, n]:
s_row += 1
if abs(i) > 1e-10:
anarbmat[r, :] = arbmat[s_row + r, :]
for j in range(r, row_number):
if j < s_row + r:
anarbmat[j + 1, :] = arbmat[j, :]
arbmat = anarbmat.copy()
if abs(anarbmat[r, n]) > 1e-10:
anarbmat[r, :] = anarbmat[r, :] / anarbmat[r, n]
for i in range(row_number):
if i != r:
anarbmat[i, :] -= \
anarbmat[i, n] * anarbmat[r, :]
arbmat = anarbmat.copy()
if abs(arbmat[r, n]) < 1e-10:
r = r
else:
r = r + 1
for m in range(column_number):
if abs(arbmat[-1, m]) > 1e-10:
arbmat[-1, :] = arbmat[-1, :] / arbmat[-1, m]
for i in range(row_number - 1):
arbmat[i, :] -= \
arbmat[i, m] * arbmat[-1, :]
break
return arbmat
def backsubstitution(B):
"""
return the reduced row echelon form matrix of B
"""
if not inconsistentSystem(B):
return rsmat(B)
#####################
# 测试的值
a = np.array([[1, 2, 3, 24], [5, 4, 6, 24], [10, 9, 8, 9]])
print(backsubstitution(a))
# 最后的结果应该是[1,0,0,-8],[0,1,0,1],[0,0,1,10]
def backsubstitution(B):
"""
return the reduced row echelon form matrix of B
"""
"""
return the reduced row echelon form matrix of B
"""
# Get a copy of the matrix
A = B.copy()
# Get the rows and columns
m, n = np.shape(B)
# Go through each row, starting at the bottom
for i in range(m-1, -1, -1):
# Set the pivot to 0.0 and location of pivot to -1.0
pivot = 0.0
pivotLocation = -1.0
# Go through each column
for h in range(n):
# Get if the current item is zero or not
isNonZero =np.nonzero(A[i][h])
# If it's not zero and the pivot is (not set yet)
if isNonZero and pivot == 0.0:
# Save the current item as pivot, and it's location as pivotLocation
pivot = A[i][h]
pivotLocation = h
# If the pivot is not 0.0 (it has been set)
if not(pivot == 0.0):
# Divide the current item by the pivot
A[i][h] = A[i][h] / pivot
# Go through each row above the current row (going up)
for j in range(i - 1, -1, -1):
# Reduce the the row by the current row and the pivot location
rowReduce(A, i, j, pivotLocation)
# Return the matrix
return A
在最后一个backsubstitution函数中,出现了什么错误,导致了结果不同呀.
def backsubstitution(B):
import sympy as sp
ma = sp.Matrix(B)
return ma.rref()[0]
if __name__ == '__main__':
a = np.array([[1, 2, 3,24], [5, 4, 6,24], [10, 9, 8,9]])
res = backsubstitution(a)
print(res)
def backsubstitution(B):
"""
return the reduced row echelon form matrix of B
"""
import sympy # 导入科学计算包sympy
ma = sympy.Matrix(B) # 生成B的矩阵
return ma.rref()[0] # 将矩阵简化为Row梯形形式
import sys
import numpy as np
import warnings
def swapRows(A, i, j):
"""
interchange two rows of A
operates on A in place
"""
tmp = A[i].copy()
A[i] = A[j]
A[j] = tmp
def relError(a, b):
"""
compute the relative error of a and b
"""
with warnings.catch_warnings():
warnings.simplefilter("error")
try:
return np.abs(a-b)/np.max(np.abs(np.array([a, b])))
except:
return 0.0
def rowReduce(A, i, j, pivot):
"""
reduce row j using row i with pivot pivot, in matrix A
operates on A in place
"""
factor = A[j][pivot] / A[i][pivot]
for k in range(len(A[j])):
if np.isclose(A[j][k], factor * A[i][k]):
A[j][k] = 0.0
else:
A[j][k] = A[j][k] - factor * A[i][k]
# stage 1 (forward elimination)
def forwardElimination(B):
"""
Return the row echelon form of B
"""
A = B.copy().astype(float)
m, n = np.shape(A)
for i in range(m-1):
# Let lefmostNonZeroCol be the position of the leftmost nonzero value
# in row i or any row below it
leftmostNonZeroRow = m
leftmostNonZeroCol = n
## for each row below row i (including row i)
for h in range(i,m):
## search, starting from the left, for the first nonzero
for k in range(i,n):
if (A[h][k] != 0.0) and (k < leftmostNonZeroCol):
leftmostNonZeroRow = h
leftmostNonZeroCol = k
break
# if there is no such position, stop
if leftmostNonZeroRow == m:
break
# If the leftmostNonZeroCol in row i is zero, swap this row
# with a row below it
# to make that position nonzero. This creates a pivot in that position.
if (leftmostNonZeroRow > i):
swapRows(A, leftmostNonZeroRow, i)
# Use row reduction operations to create zeros in all positions
# below the pivot.
for h in range(i+1,m):
rowReduce(A, i, h, leftmostNonZeroCol)
return A
####################
# If any operation creates a row that is all zeros except the last element,
# the system is inconsistent; stop.
def inconsistentSystem(A):
# 判断一个矩阵是否是不相容的,只能用numpy的nonzero方法
index = np.nonzero(A)
count = np.count_nonzero(index, axis=1)
if 0 in count[:-1]:
print("the matrix is inconsistent, stop")
sys.exit(-1)
def backsubstitution(B):
"""
return the reduced row echelon form matrix of B
"""
A = forwardElimination(B)
inconsistentSystem(A)
N = A.shape[0]
for i in reversed(range(0, N)):
for j in range(i+1, N):
A[i] -= A[i][j] * A[j]
A[i] /= A[i][i]
print("Solution for the system:")
print(A)
#####################
# 测试的值
a = np.array([[1, 2, 3,24], [5, 4, 6,24], [10, 9, 8,9]])
backsubstitution(a)
# 最后的结果应该是[1,0,0,-8],[0,1,0,1],[0,0,1,10]