【国科大——矩阵分析与应用】使用高斯消元法，测试二元一次方程系数产生的误差

使用高斯消元法，测试二元一次方程系数产生的误差

原始方程

$$\begin{array}{c}\{\begin{array}{l}47x+28y=19\\ 89x+53y=36\end{array}\end{array}$$

误差来源

对于上述方程组来说，肉眼可解：x=1，y=-1

然而，高斯消元会出现分数 $\frac{89}{47}$ 。 由于计算机的精度不同，会出现误差，下面测试不同的精度导致使用 高斯消去法 得到的结果不同。

保留有效数字

print(float(format(52.91, '.3g')))  # 保留3位有效数字
1

代码

import numpy as np
from numpy import *


def Gauss(Mat, Bais, num):
    if (Mat.shape[0] != Mat.shape[1]) or (Mat.shape[0] != Bais.shape[0]):
        raise ValueError("size of matrix doesn't match!")
    A = Mat.copy()
    B = Bais.copy()
    n = A.shape[0]  # n=2
    arg = '.' + str(num) + 'g'
    l_ik = float(format(A[1, 0] / A[0, 0], arg))
    # print(l_ik)

    a1 = float(format(A[0, 0], arg))
    b1 = float(format(A[0, 1], arg))
    a2 = float(format(A[1, 0], arg))
    b2 = float(format(A[1, 1], arg))
    c1 = float(format(B[0, 0], arg))
    c2 = float(format(B[1, 0], arg))
    b2 = float(format(b2 - float(format(l_ik * b1, arg)), arg))
    # print('b2', b2)
    c2 = float(format(c2 - float(format(l_ik * c1, arg)), arg))
    y = float(format(c2 / b2, arg))
    # print(y)
    c1 = float(format(c1 * l_ik, arg))
    b1 = float(format(b1 * l_ik * y, arg))
    # print('b1',b1)
    # print('c1',c1)
    temp = float(format(c1 - b1, arg))
    # print('temp:',temp)
    x = float(format(temp / a2, arg))
    # print('x:',x)
    return [x, y]


if __name__ == '__main__':
    x = np.mat([[47, 28],
                [89, 53]])
    y = np.mat([[19],
                [36]])
    print(Gauss(x, y, 4))

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43