三角分解法（线性方程求解）python

 

 

【问题描述】为求解一个线性方程组，首先采用偏序选主元策略的三角分解法构造矩阵L，U和P，再用前向替换法对方程组LY=PB求解Y，最后用回代法对方程组UX=Y求解X。

【输入形式】在屏幕上依次输入方阵阶数n，系数矩阵A和常数矩阵B。

【输出形式】先输出LU分解结果，再输出方程解。

【样例1输入】

4

1 2 4 1

2 8 6 4

3 10 8 8

4 12 10 6

21

52

79

82

【样例1输出】

[[ 4.   12.   10.    6.  ]

 [ 0.5   2.    1.    1.  ]

 [ 0.25 -0.5   2.    0.  ]

 [ 0.75  0.5   0.    3.  ]]

[[1.]

[2.]

[3.]

[4.]]

【样例1说明】输入：第1行为方阵阶数4，第2行至5行为系数矩阵A，第6行至9行为常数矩阵B。输出：第1至第4行输出LU分解结果，第5行至第8行依次输出方程解：x1, x2, x3, x4。

 

代码：

import numpy as np
#print("请输入矩阵的行列数：")
input001= list(map(int, input().split()))
n = int(input001[0])
#print("请输入矩阵NxN:")
a = np.zeros((n,n),dtype = np.double)
for r in range(n):
    a[r,:] = np.array(input().split(),dtype = np.double)
#print("输入常数矩阵b:")
b = np.zeros((n,1),dtype= np.double)
for r in range(n):
    b[r] = np.array(input(),dtype=np.double)

#for循环求出LU分解结果
for i in range(n):
    max = i#换行标记符
    for i1 in range(i,n):
        #找到该列元素中的最大值
        if(a[max][i]<a[i1][i]):
            max = i1
    #若没有交换行元素
    if(max == i):
        for k in range(i,n-1):
            x = a[k+1][i]/a[i][i]
            #给下三角的一个元素赋值
            a[k+1][i] = x
            #计算上三角i+1行的右边值
            for m in range (i+1,n):
                a[k+1][m] = a[k+1][m] - x*a[i][m]
    else:
        #a矩阵交换两行，是后n-i行
        for i2 in range (i,n):
            temp = a[i][i2]
            a[i][i2] = a[max][i2]
            a[max][i2] = temp
        #b矩阵交换两行
        temp = b[max][0]
        b[max][0] = b[i][0]
        b[i][0] = temp
        #a矩阵交换两行，是前i行
        for i3 in range(0,i):
            temp = a[i][i3]
            a[i][i3] = a[max][i3]
            a[max][i3] = temp
        for k in range(i,n-1):
            x = a[k+1][i]/a[i][i]
            #a1[k+1][i] = x
            a[k+1][i] = x
            for m in range (i+1,n):
                a[k+1][m] = a[k+1][m] - x*a[i][m]
print(a)
#使用向前替代法求解方程
#第一个b[0]不用计算
for i in range (1,n):
    x = 0
    for i1 in range (0,i):
        x = a[i][i1]*b[i1][0] + x
    b[i] = (b[i][0] - x)
#注意到b的元素提取用b[i][0],b[i]代表一行
#回代法
b[n-1][0] = b[n-1][0]/a[n-1][n-1]
for i in range(n-2,-1,-1):
    x = 0
    for i1 in range (i,n-1):
        x = x + a[i][i1+1]*b[i1+1][0]
        #print(x)
    b[i][0] = (b[i][0]-x)/a[i][i]
print(b)