LU分解-三角分解PYTHON实现
作者:互联网
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,读者可以评论区留言。
标签:PYTHON,s1,LU,主元,item,range,分解,ls,sc 来源: https://blog.csdn.net/m0_56264007/article/details/121051072