LU分解-三角分解PYTHON实现

2021-11-15 03:36:36

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代码。该法为列主元的三角分解法，避免大数吃小数的情况使结果更加接近精确值。

‘’‘大数吃小数’‘’指的是在计算机计算中相差数量级较大的两数相加时，直接忽略小数的存在，该现象的存在会导致某些矩阵解误差较大。

下面举的例子说明了大数吃小数现象的存在，但并未说明换主元法对该现象的改善（抱歉没找到合适的例子）