我终于学会了将矩阵化为行阶梯形和行最简形(使用python) |
您所在的位置:网站首页 › 行阶梯形矩阵第一个数必须是1吗 › 我终于学会了将矩阵化为行阶梯形和行最简形(使用python) |
如下图:看似简单的矩阵,要化成行阶梯形时,我们直接看参考答案,就发现手算不是个容易的事,于是我就思考怎么用python将一个一般的矩阵转化为行阶梯形和行最简形。 题目答案首先,做事情要有一个框架。 让我们先复习一下将一个一般的矩阵转化为行阶梯形和行最简形的基本思路。 一般的矩阵Am*n,我们常使用三种初等行变换先化成行阶梯形,再化成行最简形。 化成行阶梯形的方式是若矩阵第一行第一个元素非0,用此元素通过初等行变换消去第一列的其他元素,然后如果第二行第二个元素非0,用此元素通过初等行变换消去第二列中第三及以后的元素,以此类推;若目标位置出现为0时,向下寻找非0元素的行,然后交换这两行,转为非0的情况;若目标位置及下面的元素都为0时,将目标位置向右移动一位。以此类推直至到达最后一列或最后一行。 化成行最简形的方式是在行阶梯形的基础上,从最后一行开始,用非零首元通过初等行变换消去该列非零首元以上的元素,直至第一行。 接着,要实现这个功能,我们尝试使用自顶向下设计的思路。 我们需要定义化行阶梯形和行最简形这两个大的函数。子函数行阶梯形需要定义某元素以下清零的函数,用于将主元以下的元素化零;子函数行最简形需要定义某元素以上清零的函数,用于将主元以上的元素化零。定义某列清零的函数要使用到初等行变换,因此需要定义三个初等行变换的函数。 大致思路如上,接下来自底向上编写代码。 考虑数据类型以及之间的转化,python常用组合数据类型有列表、字典、元组、集合。最合适的当然是可随意改变元素的列表。因此使用二维列表表示矩阵。 同时考虑以下两个方面: ①为了增强程序的可读性,在每次初等变换之后要将矩阵打印出来(毕竟作业本上把过程书写好再交会更好)。由于本次使用列表套列表的形式来表示矩阵,正常打印是一串长列表,因此要把格式化打印矩阵列表的函数定义出来。 ②python的数据类型没有分数,经过触发运算之后结果都是浮点数,其中有很多我们不能直接我们看出是几分之几的数字。(意味着我们不能将这个结果写到作业本上)因此,我将使用先倍乘再倍加的方式化0,这样就可以保证都是整数。为了简化程序,我将行倍加和行倍乘变成一个函数。 行倍加乘变换但是用此方法会使计算的数字增大,因此我们增加行倍乘变换的变体函数:行约分函数,而“约分”又需要先求出这一行的所有数字的最大公因数,用到求最大公因数的函数(这个程序设计课应该都学过,直接上图) 求最大公因数行“约分”三个初等行变换,还剩下一个互换两行,代码如下: 互换两行函数然后定义用主元以下/上化0的函数。 主元以下化零函数主元以上化零函数最后定义主元化一的函数 主元化一函数准备工作完毕,定义化行阶梯形和行最简形的函数 行阶梯形函数行最简形函数一个main函数将上面包装一下。 main函数这样一个行阶梯形和行最简形化简的函数就整完了。 总代码见后文,我们先测试以下文章开头的矩阵化简,结果如下。 化行阶梯形测试测试结果和参考答案一致,甚至发现参考答案中第三行没有化到最简。下面是行最简形化简 化行最简形测试可见这个代码基本的功能还是可以的,接下来就是输入输出以及提示信息的完善了。 #线性代数初等行变换、化一般矩阵为行阶梯形、行最简形工具 import copy def P(A): '''格式化打印矩阵函数''' m=len(A)#A的行数 n=len(A[0])#A的列数 for i in range(m): for j in range(n): print(A[i][j],end='\t') print()
def 最大公因数(m,n): m,n=abs(m),abs(n) d=1 for i in range(m,0,-1): if m%i==0 and n%i==0: d=i break return d def 行约分(A,i): '''尝试对第i行进行“约分”操作''' n=len(A[0])#A的列数 g=A[i-1][0] for j in range(1,n): if A[i-1][j]!=0: g=最大公因数(A[i-1][j],g) if g!=1 and g!=0: #行全零时,未进行最大公因数求解操作,g=0; #g=1时说明没有可约分的公因数 for j in range(n): if A[i-1][j]!=0: A[i-1][j]=int(A[i-1][j]/g) print('{}:1/{}*r{}结果为'.format('矩阵',g,i)) P(A)
def 行倍加乘变换(A,i,j,k1,k2): '''第i行的k1倍加第j行的k2倍''' m=len(A)#A的行数 n=len(A[0])#A的列数 for t in range(n):#遍历第i行的所有元素 A[i-1][t]=k1*A[i-1][t]+k2*A[j-1][t] if k1==1: if k2>0: print('{}:r{}+{}*r{}结果为'.format('矩阵',i,k2,j)) else: print('{}:r{}-{}*r{}结果为'.format('矩阵',i,-k2,j)) else: if k2>0: print('{}:{}*r{}+{}*r{}结果为'.format('矩阵',k1,i,k2,j)) else: print('{}:{}*r{}-{}*r{}结果为'.format('矩阵',k1,i,-k2,j)) P(A)#格式化输出变换后的矩阵 行约分(A,i)
def 互换两行(A,i,j): n=len(A[0])#A的列数 for t in range(n): A[i-1][t],A[j-1][t]=A[j-1][t],A[i-1][t] print('{}:r{}⬅→r{}结果为'.format('矩阵',i,j)) P(A)
def 主元以下化零(A,i,j):#第i行第j列 '''将第i行第j列元素以下的元素化为零''' m=len(A)#A的行数 n=len(A[0])#A的列数 k1=A[i-1][j-1]#A中第i行第j列的元素 for s in range(i+1,m+1):#当i=m时不执行此循环 k2=A[s-1][j-1] 行倍加乘变换(A,s,i,k1,-k2) #第s行的k1(第i行第j列的元素)倍减去第i行的k2(第s行第j列的元素)倍 #s遍历第i行至最后一行
def 主元以上化零(A,i,j):#第i行第j列 '''将第i行第j列元素以上的元素化为零''' m=len(A)#A的行数 n=len(A[0])#A的列数 k1=A[i-1][j-1]#A中第i行第j列的元素 for s in range(i-1,0,-1):#当i=1时不执行此循环 k2=A[s-1][j-1] 行倍加乘变换(A,s,i,k1,-k2) #第s行的k1(第i行第j列的元素)倍减去第i行的k2(第s行第j列的元素)倍 #s遍历第i-1行至第一行 def 主元化一(A): m=len(A)#A的行数 n=len(A[0])#A的列数 A1=copy.deepcopy(A) for i in range(m): f=0 for j in range(n): if A[i][j]!=0: SouYuan=A[i][j] f=1 break if f==0: pass else: if SouYuan==1: pass elif SouYuan==-1: A1[i][j]=1 for s in range(j+1,n): A1[i][s]=-A[i][s] else: A1[i][j]=1 for s in range(j+1,n): if A1[i][s]==0: pass else: g=最大公因数(A[i][j],A[i][s]) fz,fm=int(A[i][s]/g),int(A[i][j]/g) if fmm: LQ0=1 break 互换两行(A,i,i+k) k+=1 if LQ0==1: continue 主元以下化零(A,i,t) i=i+1 if i+1>m: break
def 行最简形(A): print('-'*10+'以下为行最简形化简'+'-'*10) m=len(A)#A的行数 n=len(A[0])#A的列数 i=m while i!=1: for t in range(1,n+1): if A[i-1][t-1]==0: continue 主元以上化零(A,i,t) break i=i-1 主元化一(A) def main(A): 行阶梯形(A) 行最简形(A)
|
今日新闻 |
点击排行 |
|
推荐新闻 |
图片新闻 |
|
专题文章 |
CopyRight 2018-2019 实验室设备网 版权所有 win10的实时保护怎么永久关闭 |