python 实现高斯消元法
2018-03-01 11:19
309 查看
步骤1 检查A,b是否行数相同步骤2 构造增广矩阵Ab步骤3 逐列转换Ab为化简行阶梯形矩阵 中文维基链接
from decimal import Decimal
from fractions import Fraction
# TODO 实现 Gaussain Jordan 方法求解 Ax = b
""" Gaussian Jordan 方法求解 Ax = b.
参数
A: 方阵
b: 列向量
decPts: 四舍五入位数,默认为4
epsilon: 判读是否为0的阈值,默认 1.0e-16
返回列向量 x 使得 Ax = b
返回None,如果 A,b 高度不同
返回None,如果 A 为奇异矩阵
"""
def gj_Solve(A, b, decPts=4, epsilon = 1.0e-16):
if len(A)!= len(b) :
return None
# 构造扩增矩阵
result = augmentMatrix(A, b)
# 转换为阶梯型矩阵
for j in range(len(result[0])-1):
row,maxNum = 0,0
for i in range(len(result)):
if i >= j:
if abs(result[i][j]) > maxNum:
maxNum = abs(result[i][j])
row = i
# 返回奇异矩阵
if(abs(maxNum) < epsilon):
return None
else:
# 找出最大放到最上面
if(row != j):
swapRows(result,j,row)
row,maxNum = 0,0
scaleRow(result,j,Fraction(1,result[j][j]))
# 消除下面
for dis in range(len(result)):
if dis != j:
if abs(-result[dis][j]) > epsilon:
addScaledRow(result,dis,j,-result[dis][j])
newResult =[]
for i in range(len(result)):
row =[]
row.append(result[i][len(result)])
newResult.append(row)
return newResult测试 def test_gj_Solve(self):
for _ in range(9999):
r = np.random.randint(low=3,high=4)
A = np.random.randint(low=-10,high=10,size=(r,r))
b = np.arange(r).reshape((r,1))
x = gj_Solve(A.tolist(),b.tolist(),epsilon=1.0e-8)
if np.linalg.matrix_rank(A) < r:
self.assertEqual(x,None,"Matrix A is singular")
else:
self.assertNotEqual(x,None,"Matrix A is not singular")
self.assertEqual(np.array(x).shape,(r,1),"Expected shape({},1), but got shape{}".format(r,np.array(x).shape))
Ax = np.dot(A,np.array(x))
loss = np.mean((Ax - b)**2)
self.assertTrue(loss<0.1,"Bad result.")
对于Ab的每一列(最后一列除外) 当前列为列c 寻找列c中 对角线以及对角线以下所有元素(行 c~N)的绝对值的最大值 如果绝对值最大值为0 那么A为奇异矩阵,返回None (你可以在选做问题2.4中证明为什么这里A一定是奇异矩阵) 否则 使用第一个行变换,将绝对值最大值所在行交换到对角线元素所在行(行c) 使用第二个行变换,将列c的对角线元素缩放为1 多次使用第三个行变换,将列c的其他元素消为0步骤4 返回Ab的最后一列注: 我们并没有按照常规方法先把矩阵转化为行阶梯形矩阵,再转换为化简行阶梯形矩阵,而是一步到位。如果你熟悉常规方法的话,可以思考一下两者的等价性。
from decimal import Decimal
from fractions import Fraction
# TODO 实现 Gaussain Jordan 方法求解 Ax = b
""" Gaussian Jordan 方法求解 Ax = b.
参数
A: 方阵
b: 列向量
decPts: 四舍五入位数,默认为4
epsilon: 判读是否为0的阈值,默认 1.0e-16
返回列向量 x 使得 Ax = b
返回None,如果 A,b 高度不同
返回None,如果 A 为奇异矩阵
"""
def gj_Solve(A, b, decPts=4, epsilon = 1.0e-16):
if len(A)!= len(b) :
return None
# 构造扩增矩阵
result = augmentMatrix(A, b)
# 转换为阶梯型矩阵
for j in range(len(result[0])-1):
row,maxNum = 0,0
for i in range(len(result)):
if i >= j:
if abs(result[i][j]) > maxNum:
maxNum = abs(result[i][j])
row = i
# 返回奇异矩阵
if(abs(maxNum) < epsilon):
return None
else:
# 找出最大放到最上面
if(row != j):
swapRows(result,j,row)
row,maxNum = 0,0
scaleRow(result,j,Fraction(1,result[j][j]))
# 消除下面
for dis in range(len(result)):
if dis != j:
if abs(-result[dis][j]) > epsilon:
addScaledRow(result,dis,j,-result[dis][j])
newResult =[]
for i in range(len(result)):
row =[]
row.append(result[i][len(result)])
newResult.append(row)
return newResult测试 def test_gj_Solve(self):
for _ in range(9999):
r = np.random.randint(low=3,high=4)
A = np.random.randint(low=-10,high=10,size=(r,r))
b = np.arange(r).reshape((r,1))
x = gj_Solve(A.tolist(),b.tolist(),epsilon=1.0e-8)
if np.linalg.matrix_rank(A) < r:
self.assertEqual(x,None,"Matrix A is singular")
else:
self.assertNotEqual(x,None,"Matrix A is not singular")
self.assertEqual(np.array(x).shape,(r,1),"Expected shape({},1), but got shape{}".format(r,np.array(x).shape))
Ax = np.dot(A,np.array(x))
loss = np.mean((Ax - b)**2)
self.assertTrue(loss<0.1,"Bad result.")
相关文章推荐
- 高斯消元法的MATLAB和PYTHON实现
- 高斯消元法(三):用Python简单实现顺序消元法
- 高斯消元法的python实现
- Python函数的周期性执行实现方法
- 在C#中实现Python的分片技术
- ICTClAS2013(NLPIR) 的python接口实现
- python实现rsa加密实例详解
- Fibonacci(斐波那契)数列的递归与非递归实现 python
- PYTHON实现刷流量工具
- python实现SimpleHTTPServer的POST方法
- 二进制中1的个数[剑指offer]之python实现
- 教你使用python实现微信每天给女朋友说晚安
- 用Python实现RSA签名和验签
- IP代理池的Python实现
- [Python数据结构] 使用List实现Stack
- 借助apktool.jar工具,使用python代码简化批量反编译apk安装包的简单实现
- Kmeans python 实现
- python学习——用lambda实现斐波那契函数
- Python3之实现单例模式de几种方式
- python实现删除文件与目录的方法