数值分析第五章节 用Python实现解线性方程组的直接方法

参考书籍：数值分析 第五版 李庆杨 王能超 易大义编 第5章 解线性方程组的直接解法
文章声明：如有发现错误，欢迎批评指正

引言与预备知识

5.1.1：关于线性方程组的数值解法一般有两类：直接法和迭代法。5.1.2向量和矩阵：非奇异矩阵：设$A\in R^{n\times n}$，$B\in R^{n\times n}$。如果$AB=BA=I$，则称$B$是$A$的逆矩阵，记为$A^{-1}$，且$(A^{-1}{)}^T=(A^T{)}^{-1}$。如果$A^{-1}$存在，则称$A$为非奇异矩阵。如果$A，B\in R^{n\times n}$均为非奇异矩阵，则$(AB{)}^{-1}=B^{-1}{A}^{-1}$。其他跳过。5.1.3矩阵的特征值与谱半径：$A$的全体特征值称为$A$的谱，记作$\sigma (A)$，即 σ ( A ) = { λ 1 , λ 2 , … , λ n } \sigma(A)=\{\lambda_1,\lambda_2,\dots,\lambda_n\} σ(A)={λ1​,λ2​,…,λn​}。记$\rho (A)=\underset{1\le i\le n}{\max }|{\lambda }_i|$，称为矩阵$A$的谱半径。其他跳过。5.1.4特殊矩阵：直接跳过。都是高等代数知识，所以我们直接跳过。

高斯消去法

5.2.1高斯消去法：十分简单大家都会，注意一下使用条件。定理：约化的主元素$a_{ii}^{(i)}\ne 0(i=1,2,\dots ,k)$的充要条件是矩阵$A$的顺序主子式$D_i\ne 0(i=1,2,\dots ,k)$(挺重要的，挺好证的)。矩阵的初等行变换其实就是左乘一个矩阵。$L^{(1)}{A}^{(0)}=A^{(1)},\dots ,L^{(n)}{A}^{(n-1)}=A^{(n)}。L=L^{(n-1)}\times \cdots \times L^{(1)}，U=A^{(n)}$。 5.2.3列主元法：我们使用列主元法进行演示。5.2.2矩阵的三角分解：定理(矩阵的LU分解)：设$A$为$n$阶矩阵，如果$A$的顺序主子式$D_i\ne 0(i=0,1,\dots ,n-1)$，则$A$可分解为一个单位下三角矩阵$L$和一个上三角矩阵$U$的乘积，且这种分解是唯一的。5.2.1，5.2.3，5.2.2的顺序更有条理。

列主元消去法

我感觉我写得挺好，可以算作通用代码，前提必须保证有解。caozuo1是交换，caozuo2是归0。你把caozuo1删除了，就是高斯消去。这是第五章节例4，书上每步全都采用4位近似，我的更加准确一点。

def caozuo1(lt,x):
    global m
    a,b=x,abs(lt[x][x])
    for i in range(x+1,m):
        if abs(lt[i][x])>b:
            a,b=i,abs(lt[i][x])
    lt[x],lt[a]=lt[a],lt[x]
    return lt
def caozuo2(lt,x):
    global m,n
    for i in range(x+1,m):
        k=-lt[i][x]/lt[x][x];lt[i][x]=0
        for j in range(x+1,n):
            lt[i][j]+=lt[x][j]*k
    return lt
def solve(lt):
    global m,n
    x=[]
    for i in range(m-1,-1,-1):
        cnt=lt[i][n-1]
        for j in range(len(x)):
            cnt-=lt[i][n-j-2]*x[j]
        x.append(cnt/lt[i][i])
    return x[::-1]
m,n=map(eval,input().strip().split())
lt=[]
for _ in range(m):
    lt.append([eval(_) for _ in input().strip().split()])
for i in range(m-1):
    lt=caozuo1(lt,i)
    lt=caozuo2(lt,i)
print(solve(lt))

矩阵三角分解法

5.3.1直接三角分解法：1.不选主元的三角分解法。(杜利特尔分解法)2.选主元的三角分解法。我们使用杜利特尔分解法进行演示。
我感觉我写得挺好，可以算作通用代码，前提必须保证有解。这是第五章节例5，与书上的一摸一样。

杜利特尔分解法

def L(lt1,lt2,x):
    for i in range(x+1,n):
        cnt=0
        for j in range(x):
            cnt+=lt1[i][j]*lt2[j][x]
        lt1[i][x]=(lt[i][x]-cnt)/lt2[x][x]
    return lt1
def U(lt1,lt2,x):
    global lt,n
    for i in range(x,n):
        cnt=0
        for j in range(x):
            cnt+=lt1[x][j]*lt2[j][i]
        lt2[x][i]=lt[x][i]-cnt
    return lt2
n=int(input("系数方阵的行列数"))
lt1=[[0 for _ in range(n)] for _ in range(n)]
for _ in range(n):
    lt1[_][_]=1
lt2=[[0 for _ in range(n)] for _ in range(n)]
lt=[]
for _ in range(n):
    lt.append([eval(_) for _ in input().strip().split()])
for i in range(n-1):
    lt2=U(lt1,lt2,i)
    lt1=L(lt1,lt2,i)
lt2=U(lt1,lt2,n-1)
print("===============L===============")
for i in lt1:
    for j in i:
        print("{:>8.2f}".format(j),end="")
    print()
print("===============U===============")
for i in lt2:
    for j in i:
        print("{:>8.2f}".format(j),end=" ")
    print()

我们使用之前写的列主元消去法求解。y，x如下。
y：

x：

5.3.2平方根法。定理：设$A$为$n$阶对称矩阵，且$A$的所有顺序主子式均不为零，则$A$可分解为$A=LD{L}^T$，其中$A$为单位下三角矩阵，$D$为对角矩阵。定理：如果$A$为$n$阶对称正定矩阵，则存在一个实的非奇异下三角矩阵$L$，使$A=L{L}^T$，当限定$L$的对角元素为正时，这种分解是唯一的。我们对它进行演示。

平方根法

自己随便捏的对称正定矩阵。代码变量LT是L的转置意思。

def caozuo(LT,x):
    global lt,n
    cnt=0
    for i in range(x):
        cnt+=pow(LT[i][x],2)
    LT[x][x]=pow(lt[x][x]-cnt,0.5)
    for i in range(x+1,n):
        cnt=0
        for j in range(x):
            cnt+=LT[j][x]*LT[j][i]
        LT[x][i]=(lt[x][i]-cnt)/LT[x][x]
    return LT
n=int(input("系数矩阵的行列数"))
lt=[]
for _ in range(n):
    lt.append([eval(_) for _ in input().strip().split()])
LT=[[0 for _ in range(n)] for _ in range(n)]
for i in range(n):
    LT=caozuo(LT,i)
for i in LT:
    for j in i:
        print("{:>8.2f}".format(j),end=" ")
    print()

5.3.3 追赶法：追赶法更简单，我们就不演示。

向量和矩阵的范数

5.4.1向量范数：向量的$p$范数：$||x|{|}_p=(\sum _{i=1}^n|x_i{|}^p{)}^{\frac{1}{p}}，p\in [1,\infty )$。其他东西我们不管。5.4.2 矩阵范数：$||A|{|}_{\infty }=\underset{1\le i\le n}{\max }\sum _{j=1}^n|a_{ij}|$(行范数)。$||A|{|}_1=\underset{1\le j\le n}{\max }\sum _{i=1}^n|a_{ij}|$(列范数)。$||A|{|}_2=\sqrt{{\lambda }_{\max }(A^{T}A)}$。其他东西我们不管。求范数的程序简单所以我们不作演示。

误差分析

定义：设$A$为非奇异矩阵，称数$cond(A{)}_v=||A^{-1}|{|}_v||A|{|}_v(v=1,2,\infty )$为矩阵A的条件数。$cond(A{)}_2=||A^{-1}|{|}_2||A|{|}_2=\sqrt{\frac{{\lambda }_{\max }(A^{T}A)}{{\lambda }_{\min }(A^{T}A)}}$当$A$为对称矩阵时$cond(A{)}_2=\frac{|{\lambda }_1|}{|{\lambda }_2|}$，$cond(A{)}_{\infty }=||A^{-1}|{|}_{\infty }||A|{|}_{\infty }$。$cond(A)$越大方程越病态。