s=list(eval(input('lu'))) i=len(s)#增广矩阵行数 j=len(s[0])#增广矩阵列数 def primary(s1,a): i=len(s1)#该步运算行数 j=len(s1[0])#该步运算列数:j=1+1 global s #该步运算第一列元素按照大小编号,0对应绝对值最大元素,以用来做本步运算主元 ls=[] ls1=[] cnt={} for item in range(i): ls.append(abs(s1[item][0])) ls1.append(s1[item][0]) ls.sort(reverse=True) for item in ls1: cnt[item]=ls.index(abs(item)) #换行:绝对值最大元素放到主元位置 s11=s.copy() for item in range(i): first=ls1[item] if cnt[first]==0: s[a]=s11[item+a] s[item+a]=s11[a] break def tobe1(s1,a):#主元下方元素除以主元 for item in range(a+1,i): s1[item][a]/=s1[a][a] def mathi(s1,a,b):#计算主元所在行主元右侧元素 if a!=1: for item in range(a): s1[a][b]-=s1[item][b]*s1[a][item] return s1[a][b] else: s1[1][b]-=s1[0][b]*s1[1][0] return s1[1][b] def mathj(s1,a,b):#计算主元所在列主元下方元素 if a!=1: for item in range(b): s1[a][b]-=s1[item][b]*s1[a][item] return s1[a][b] else: s1[1][b]-=s1[0][b]*s1[1][0] return s1[1][b] def news(s1):#创造一个剔除上一步运算的矩阵 s1=s1[1:] for item in range(len(s1)): s1[item]=s1[item][1:] return s1 def makeatry(s1,a):#找到下步运算的主元 for item in range(a,i): s1[item][a]=mathj(s1,item,a) for item in range(a): s1=news(s1) primary(s1,a) #计算LU分解后的矩阵 primary(s,0) tobe1(s,0) for item in range(1,i-1): makeatry(s,item) tobe1(s,item) for jtem in range(item+1,j): mathi(s,item,jtem) for item in range(j-2,j): mathi(s,j-2,item) sc=s.copy() print(s) #输出方程组的解,由最后向最前储存在ls里 ls=[] def backout(sc,i): sc[i-1][i-1]=sc[i-1][i]/sc[i-1][i-1] ls.append(sc[i-1][i-1]) for item in range(i-1): sc[item][i-1]=sc[item][i]-sc[item][i-1]*sc[i-1][i-1] sc=sc[:i-1] for item in range(i-1): sc[item]=sc[item][:i] i=len(sc) if i==0: return ls backout(sc,i) backout(sc,i) for item in ls: print(item,end=' ')
在我们学习科学计算(数值分析)过程中,求解某个方程组时,常常用到LU分解(三角分解法)。我针对课本中常出现的方阵的增列矩阵求解问题,编写如上python代码。该法为列主元的三角分解法,避免大数吃小数的情况使结果更加接近精确值。
‘’‘大数吃小数’‘’指的是在计算机计算中相差数量级较大的两数相加时,直接忽略小数的存在,该现象的存在会导致某些矩阵解误差较大。
下面举的例子说明了大数吃小数现象的存在,但并未说明换主元法对该现象的改善(抱歉没找到合适的例子)
下面回顾一下选列主元的三角分解法计算步骤
上面的代码print的矩阵可以分解成LU两式,
输出的S下三角元素外加主对角线元素为1为L
输出的S上三角元素及输出S的主对角线元素素为U
下面为一个例子
精确解为
x4=4
x3=3
x2=2
x1=1
执行该代码我也遇见了问题:如果执行下方代码,S将发生该变,但该段代码并不直接作用于S而是作用于S.COPY()
执行软件为IDLE,读者可以评论区留言。