计算方法——三角分解法与高斯迭代法求解线性方程组

计算方法——三角分解法与高斯迭代法求解线性方程组

• 三角分解法
• 高斯赛德尔迭代法

一、三角分解法

【问题描述】为求解一个线性方程组，首先采用偏序选主元策略的三角分解法构造矩阵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
import math
n=input()
n=int(n)
A=np.zeros((n,n))
L=np.zeros((n,n))
U=np.zeros((n,n))
P=np.zeros((n,n))
A0=np.zeros((1,n))
B=np.zeros((n,1))
Y=np.zeros((n,1))
X=np.zeros((n,1))
for i in range(0,n):
    Enter=input()
    Enter=Enter.split(" ")
    A[i]=Enter
for i in range(0,n):
    Enter=input()
    B[i][0]=Enter
for i in range(0,n):
    P[i][i]=1
for i in range(0,n):
    max = abs(A[i][i])
    maxc=i
    for j in range(i+1,n):
        if(max<=abs(A[j][i])):
            max=abs(A[j][i])
            maxc=j
    A0[0]=A[maxc]
    A[maxc]=A[i]
    A[i] = A0[0]
    A0[0]=P[maxc]
    P[maxc]=P[i]
    P[i]=A0[0]
    for j in range(i + 1, n):
        A[j][i+1:n]=A[j][i+1:n]-A[j][i]/A[i][i]*A[i][i+1:n]
        A[j][i:i+1]=A[j][i]/A[i][i]
for i in range(0,n):
    L[i][i]=1
    for j in range(0,i):
        L[i][j]=A[i][j]
for i in range(0,n):
    for j in range(0,i+1):
        U[n-i-1][n-j-1]=A[n-i-1][n-j-1]
B=np.matmul(P,B)
Y[0]=B[0]
i=1
while(i!=n):
    Y[i]=(B[i]-np.matmul(L[i,0:i],Y[0:i]))
    i=i+1
i=n-2
X[n-1]=Y[n-1]/U[n-1][n-1]
while(i!=-1):
    X[i]=(Y[i]-np.matmul(U[i,i+1:n],X[i+1:n]))/U[i,i]
    i=i-1
print(A)
print(X)

二、高斯赛德尔迭代法

【问题描述】为求解一个线性方程组，使用高斯赛德尔迭代法，采用欧几里得距离判定是否收敛。精度delta为1E-9，最大迭代次数为20。

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

【输出形式】输出实际迭代次数，然后每一行输出一个根。

【样例1输入】

3

4 -1 1

4 -8 1

-2 1 5

7

-21

15

1

2

2

【样例1输出】

10

[[2.]

[4.]

[3.]]

【样例1说明】输入：第1行为方阵阶数3，第2行至4行为系数矩阵A，第5行至7行为常数矩阵B，第8行至10行为起始点。输出：实际迭代次数为10，然后每行依次输出方程解：x1, x2, x3。

【评分标准】根据输入得到的输出准确

import numpy as np
n=input()
n=int(n)
A=np.zeros((n,n))
B=np.zeros((n,1))
X=np.zeros((n,1))
P=np.zeros((n,1))
err=0
reller=0
count=0
for i in range(0,n):
    Enter=input()
    Enter=Enter.split(" ")
    A[i]=Enter
for i in range(0,n):
    Enter=input()
    B[i][0]=Enter
for i in range(0,n):
    Enter=input()
    P[i][0]=Enter
    X[i][0]=Enter
for i in range(0,20):
    for j in range(0,n):
        P[j][0]=X[j][0]
        X[j][0]=(B[j]-np.matmul(A[j,0:j],X[0:j])-np.matmul(A[j,j+1:n],X[j+1:n]))/A[j][j]
    temp=np.subtract(P,X)
    temp=np.fabs(temp)
    err=temp.sum(axis=0)
    reller=err/(X.sum(axis=0)+1e-12)
    if(err<1e-9 or reller<1e-9):
        break
    count=count+1
print(count)
print(X)