线性代数笔记

第1章 线性代数中的线性方程组

1.1线性方程组

线性方程组的定义

解集，三种情况：无解、唯一解、无穷多解

线性方程组可以用矩阵表示：系数矩阵、增广矩阵

解线性方程组思路：等价变换（化简）

化简线性方程组的三种基本变换：

• 这里用method,i,j,k表示，method<0:行变换;method>0:列变换。
• |method|=1表示用k乘某一行(列);|method|=2表示两行(列)互换;
• |method|=3表示将i行(列)乘以数k后加到j行(列)上。i,j均从0开始。

线性方程组的变换对应于增广矩阵的初等变换：

• 倍加变换、对换变换、倍乘变换

2个线性方程组的增广矩阵是行等价的，则它们有相同的解集（想一下为什么）

线性方程组的2个基本问题：是否相容，有解的话解是否唯一

练习题
3.数组(3,4,-2)是否为下列方程组的解？

为解线性方程组，编写代码如下，可通过transform_matrix函数对增广矩阵进行初等变换。解题步骤在代码后面。

import numpy as np

def transform_matrix(matrix,method=0,i=0,j=0,k=1):
    #method=0:不作变换;method<0:行变换;method>0:列变换
    #|method|=1表示用k乘某一行(列);|method|=2表示两行(列)互换;|method|=3表示将i行(列)乘以数k后加到j行(列)上
    #i,j均从0开始
    #1.生成初等矩阵
    if method==0:
        return matrix
    elif method<0:
        if i<0 or i>=matrix.shape[0]:
            print("error i")
        if j<0 or j>=matrix.shape[0]:
            print("error j")
    elif method>0:
        if i<0 or i>=matrix.shape[1]:
            print("error i")
        if j<0 or j>=matrix.shape[1]:
            print("error j")
    d=matrix.shape[0] if method<0 else matrix.shape[1]
    A=np.eye(d)
    if abs(method)==1:
        A[i,i]=k
    elif abs(method)==2:
        A[i,j]=1
        A[j,i]=1
        A[i,i]=0
        A[j,j]=0
    elif method==-3:
        A[j,i]=k
    elif method==3:
        A[i,j]=k
    else:
        print("method error")
        return matrix
    if method<0:
        ret=np.matmul(A,matrix)
        return ret
    elif method>0:
        ret=np.matmul(matrix,A)
        return ret
    else:
        return matrix

a=np.array([[1,2,3,4],[5,6,7,8],[9,10,11,12]]);a

array([[ 1,  2,  3,  4],
       [ 5,  6,  7,  8],
       [ 9, 10, 11, 12]])

b=transform_matrix(a,-1,2,0,10.5);b

array([[  1. ,   2. ,   3. ,   4. ],
       [  5. ,   6. ,   7. ,   8. ],
       [ 94.5, 105. , 115.5, 126. ]])

c=transform_matrix(b,-2,1,2,0);c

array([[  1. ,   2. ,   3. ,   4. ],
       [ 94.5, 105. , 115.5, 126. ],
       [  5. ,   6. ,   7. ,   8. ]])

d=transform_matrix(c,-3,0,2,10);d

array([[  1. ,   2. ,   3. ,   4. ],
       [ 94.5, 105. , 115.5, 126. ],
       [ 15. ,  26. ,  37. ,  48. ]])

练习题3解题步骤：

a=np.array([[5,-1,2,7],[-2,6,9,0],[-7,5,-3,-7]]);a

array([[ 5, -1,  2,  7],
       [-2,  6,  9,  0],
       [-7,  5, -3, -7]])

b=transform_matrix(a,-3,0,1,2/5);b

array([[ 5. , -1. ,  2. ,  7. ],
       [ 0. ,  5.6,  9.8,  2.8],
       [-7. ,  5. , -3. , -7. ]])

c=transform_matrix(b,-3,0,2,7/5);c

array([[ 5. , -1. ,  2. ,  7. ],
       [ 0. ,  5.6,  9.8,  2.8],
       [ 0. ,  3.6, -0.2,  2.8]])

d=transform_matrix(c,-3,1,2,-3.6/5.6);d

array([[ 5. , -1. ,  2. ,  7. ],
       [ 0. ,  5.6,  9.8,  2.8],
       [ 0. ,  0. , -6.5,  1. ]])

e=transform_matrix(d,-3,2,1,9.8/6.5);e

array([[ 5.00000000e+00, -1.00000000e+00,  2.00000000e+00,
         7.00000000e+00],
       [ 0.00000000e+00,  5.60000000e+00, -1.89421128e-15,
         4.30769231e+00],
       [ 0.00000000e+00,  0.00000000e+00, -6.50000000e+00,
         1.00000000e+00]])

出问题了，是舍入误差导致的，想了半天没找到从程序上根本解决的办法，只有迂回解决了：

e=transform_matrix(d,-1,1,0,65);e

array([[  5. ,  -1. ,   2. ,   7. ],
       [  0. , 364. , 637. , 182. ],
       [  0. ,   0. ,  -6.5,   1. ]])

f=transform_matrix(e,-1,2,0,98);f

array([[   5.,   -1.,    2.,    7.],
       [   0.,  364.,  637.,  182.],
       [   0.,    0., -637.,   98.]])

g=transform_matrix(f,-3,2,1,1);g

array([[ 5.00000000e+00, -1.00000000e+00,  2.00000000e+00,
         7.00000000e+00],
       [ 0.00000000e+00,  3.64000000e+02, -1.13686838e-13,
         2.80000000e+02],
       [ 0.00000000e+00,  0.00000000e+00, -6.37000000e+02,
         9.80000000e+01]])

还是不行，通过调试发现从d=transform_matrix(c,-3,1,2,-3.6/5.6)的执行过程中已产生了误差，d[2,2]=-6.500000000000001。百度关键词numpy 舍入误差，找到《numpy矩阵求逆舍入误差(numpy matrix inversion rounding errors)》，通过在程序中加入np.set_printoptions(suppress=True)即可改进显示效果。

np.set_printoptions(suppress=True)

g

array([[   5.,   -1.,    2.,    7.],
       [   0.,  364.,   -0.,  280.],
       [   0.,    0., -637.,   98.]])

h=transform_matrix(g,-3,2,0,2/637);h

array([[   5.        ,   -1.        ,   -0.        ,    7.30769231],
       [   0.        ,  364.        ,   -0.        ,  280.        ],
       [   0.        ,    0.        , -637.        ,   98.        ]])

i=transform_matrix(h,-3,1,0,1/364);i

array([[   5.        ,    0.        ,   -0.        ,    8.07692308],
       [   0.        ,  364.        ,   -0.        ,  280.        ],
       [   0.        ,    0.        , -637.        ,   98.        ]])

j=transform_matrix(i,-1,0,0,1/5);j

array([[   1.        ,    0.        ,   -0.        ,    1.61538462],
       [   0.        ,  364.        ,   -0.        ,  280.        ],
       [   0.        ,    0.        , -637.        ,   98.        ]])

k=transform_matrix(j,-1,1,0,1/364);k

array([[   1.        ,    0.        ,   -0.        ,    1.61538462],
       [   0.        ,    1.        ,   -0.        ,    0.76923077],
       [   0.        ,    0.        , -637.        ,   98.        ]])

l=transform_matrix(k,-1,2,0,-1/637);l

array([[ 1.        ,  0.        , -0.        ,  1.61538462],
       [ 0.        ,  1.        , -0.        ,  0.76923077],
       [ 0.        ,  0.        ,  1.        , -0.15384615]])

print("x1=",l[0,3]);print("x2=",l[1,3]);print("x3=",l[2,3])

x1= 1.6153846153846154
x2= 0.7692307692307689
x3= -0.15384615384615366

验证答案是否正确：

x1=1.6153846153846154
x2=0.7692307692307689
x3=-0.15384615384615366
    
y1=np.float64(5*x1-x2+2*x3-7)
y2=np.float64(-2*x1+6*x2+9*x3)
y3=np.float64(-7*x1+5*x2-3*x3+7)
print("y1=",y1,"y2=",y2,"y3=",y3)

y1= 0.0 y2= -6.661338147750939e-16 y3= -2.6645352591003757e-15

数组(3,4,-2)不是方程组的解 □