计算方法——线性方程组的数值解法

线性方程组的数值解法（python）

• 上三角线性方程组
• 高斯消去法

上三角线性方程组

1.简介
NⅹN矩阵A=[aij]中的元素满足对所有i>j，有aij =0，则称矩阵A为上三角矩阵。如果A中的元素满足，对所有i<j，有aij =0，是称矩阵A为下三角矩阵。设AX=B是上三角线性方程组，如果akk≠0, k=1,2, …, N，则称该方程组存在唯一解，如：
a11x1+a12x2+a13x3+…+a1nxn=b1
a22x2+a23x3+…+a2nxn=b2
a33x3+…+a3nxn=b3
…
an-1n-1xn-1+annxn=bn-1
annxn=bn

2.方法
①由最后一个方程可求得：xn=bn/ann
②将xn代入倒数第二个方程可求得：xn-1=(bn-1-annxn)/an-1n-1
③所以回带公式为：xk=(bk-((akk+1xk+1)+(akk+2xk+2)+…(aknxn)))/akk

3.题目（用numpy）
【问题描述】在一个上三角线性方程组基础上，进行线性方程组求解。

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

【输出形式】每一行输出一个根

【样例1输入】

4

4 -1 2 3

0 -2 7 -4

0 0 6 5

0 0 0 3

20

-7

4

6

【样例1输出】

[[ 3.]

[-4.]

[-1.]

[ 2.]]

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

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

import numpy as np
n=input()
n=int(n)
A=np.zeros((n,n))
B=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
i=n-2
X[n-1]=B[n-1]/A[n-1][n-1]
while(i!=-1):
    X[i]=(B[i]-np.matmul(A[i,i+1:n],X[i+1:n]))/A[i,i]
    i=i-1
print(X)

高斯消去法

1.简介
求解有N个方程和N个未知数的一般方程组AX=B，如果可以构造一个等价的上三角方程组UX=Y，这样可以利用回代法进行求解。

2.主元
用系数矩阵A中的元素arr用来消去ark，其中k=r+1, r+2, …, N，这里称arr为第r个主元，第r行称为主元行。

3.消元
第p列消元：若app!=0则计算mrp=arp/app，r=p+1, …, N，将第r行减去mrp与第p行的乘积，消去从第p+1行到第N行的xp.
arc=0 if c=p
arc=acr-mrp*apc if c=p+1,…,n+1

4.平凡选主元
如果 app=0 则不能使用第p行消去主对角线之下列p的元素，有必要寻找第k行，满足 akp!=0 而且k>p，然后交换第k行和第p行，使得主元非零，这个过程称为选主元，选择行的判定条件称为选主元策略。

5.偏序选主元（可减少误差）
为了减少误差的传播，偏序选主元策略首先检查位于主对角线或主对角线下第p列的所有元素，确定行k，它的元素的绝对值最大，即|akp|=max{|app|, |ap+1p|, …, |aNp|}，如果k>p，则交换第k行和第p行。

6.题目
【问题描述】为求解一个线性方程组，首先构造增广矩阵[A|B]，采用偏序选主元策略的高斯消去法变换成上三角矩阵，再执行回代过程得到解。

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

【输出形式】首先输出上三角矩阵（变换后的增广矩阵），然后每一行输出一个根

【样例1输入】

4

1 2 1 4

2 0 4 3

4 2 2 1

-3 1 3 2

13

28

20

6

【样例1输出】

[[ 4. 2. 2. 1. 20. ]

[ 0. 2.5 4.5 2.75 21. ]

[ 0. 0. 4.8 3.6 26.4 ]

[ 0. 0. 0. 3.75 7.5 ]]

[[ 3.]

[-1.]

[ 4.]

[ 2.]]

【样例1说明】输入：第1行为方阵阶数4，第2行至5行为系数矩阵A，第6行至9行为常数矩阵B。

输出：首先输出上三角矩阵（变换后的增广矩阵），然后每行依次输出方程解：x1, x2, x3, x4。

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

import numpy as np
import math
n=input()
n=int(n)
A=np.zeros((n,n))
A0=np.zeros((1,n+1))
B=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
A=np.concatenate((A,B),axis=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]
    for j in range(i + 1, n):
        A[j][i+1:n+1]=A[j][i+1:n+1]-A[j][i]/A[i][i]*A[i][i+1:n+1]
        A[j][0:i+1]=0
X[n-1]=A[n-1][n]/A[n-1][n-1]
while(i!=-1):               #上三角矩阵的迭代方法
    X[i]=(A[i][n]-np.matmul(A[i,i+1:n],X[i+1:n]))/A[i,i]
    i=i-1
print(A)
print(X)                    #根的值