首页 > 其他分享 >[数值分析]解线性方程组的直接法（从线性代数到数值分析~）

[数值分析]解线性方程组的直接法（从线性代数到数值分析~）

时间：2022-10-02 18:22:05浏览次数：87  
标签：begin end 线性方程组 数值 print 线性代数 bmatrix np array

解线性方程组的直接法

例1

\[\begin{cases} 3x_1 + 2x_2 + 5x_3 = 6 \\ -x_1 + 4x_2 + 3x_3 = 5 \\ x_1 - x_2 + 3x_3 = 1 \end{cases} \]

普通解法

代码如下：

import numpy as np
A = np.array([[3, 2, 5], [-1, 4, 3], [1, -1, 3]])
b = np.array([6, 5, 1])
x = np.linalg.solve(A, b)
print(x)

import numpy as np
A = np.array([[3, 2, 5], [-1, 4, 3], [1, -1, 3]])
b = np.array([6, 5, 1])
x = np.linalg.solve(A, b)
print(x)

[0.5 1.  0.5]

高斯若尔当消元法

高斯若尔当消元法是一种解线性方程组的方法，它是一种直接法，即不需要迭代求解，而是直接求解。高斯若尔当消元法的基本思想是将线性方程组转化为一个上三角矩阵，然后逐个求解。

另注：高斯若尔当消元法是一种直接法，或者说成是一种朴素想法，但如果初学者在网上搜的时候可能会看到一大堆定义，其中不乏这样的说法

这里采用最朴素的解法就好（采取说法二），毕竟我们的最终目的只是解方程而已

\[\begin{cases} 3x_1 + 2x_2 + 5x_3 = 6 \\ -x_1 + 4x_2 + 3x_3 = 5 \\ x_1 - x_2 + 3x_3 = 1 \end{cases} \]

对于上面的例子，我们可以将其转化为增广矩阵：

\[\begin{bmatrix} 3 & 2 & 5 & 6 \\ -1 & 4 & 3 & 5 \\ 1 & -1 & 3 & 1 \end{bmatrix} \]

然后进行初等行变换，得到：

\[\begin{bmatrix} 3 & 2 & 5 & 6 \\ 0 & 14 & 14 & 21 \\ 0 & 5 & -4 & 3 \end{bmatrix} \dots \text{第一步} \\ \]

\[\begin{bmatrix} 3 & 2 & 5 & 6 \\ 0 & 14 & 14 & 21 \\ 0 & 0 & 2 & 1 \end{bmatrix} \dots \text{第二步}\\ \]

最后利用回代法求解，得到：

\[\begin{cases} x_3 = \frac{1}{2} \\ x_2 = \frac{21 - 14 * x_3}{14} \\ x_1 = \frac{6 - 2 * x_2 - 5 * x_3}{3} \\ \end{cases} \Rightarrow \begin{cases} x_3 = 0.5 \\ x_2 = 1 \\ x_1 = 0.5 \\ \end{cases} \]

import numpy as np

# 高斯消去法
def gaussElimination(A, b):
    n = len(b)
    # 前進消去
    for k in range(0, n-1):  # 举例第一列 k = 0
        factor = A[k, k] / A[k+1:, k] # 求出之后每一行的系数
        A[k+1:, :] = A[k, :] - factor[:,np.newaxis] * A[k+1:, :]  #更新之后每一行 np.newaxis的作用是增加一个维度 (2,1)
        b[k+1:] = b[k, :] - factor[:,np.newaxis] * b[k+1:, :]  #更新之后每一行的b

    # 回代
    x = np.zeros(n)
    x[n-1] = b[n-1] / A[n-1, n-1]
    for k in range(n-2, -1, -1):
        x[k] = (b[k] - np.dot(A[k, k+1:], x[k+1:])) / A[k, k]
    return x

A = np.array([[3, 2, 5], [-1, 4, 3], [1, -1, 3]], dtype=float)
b = np.array([[6, 5, 1]], dtype=float)
b = b.T
gaussElimination(A, b)

array([0.5, 1. , 0.5])

随机测试一个方程组

\[\begin{cases} -23x_1 + 11x_2 + x_3 = 0 \\ 11x_1 - 3x_2 - 2x_3 = 3 \\ x_1 - 2x_2 + 2x_3 = -1 \end{cases} \]

A1 = np.array([[-23, 11 , 1], [11, -3, -2], [1, -2, 2]], dtype=float)
b1 = np.array([[0, 3, -1]], dtype=float)
b1 = b1.T
gaussElimination(A1, b1)

array([1., 2., 1.])

逐步求解代码如下

A = np.array([[3, 2, 5], [-1, 4, 3], [1, -1, 3]], dtype=float)
b = np.array([[6, 5, 1]], dtype=float)
# 装置b
b = b.T
tmp_f=A[0,0]/A[1:,0] # [-3.  3.]
A[1:,:] = A[0,:] - tmp_f[:,np.newaxis] * A[1:,:] # A = [[ 3.  2.  5.][ 0.  14.  14.][ 0.  5.  -4.]]
b[1:,:] = b[0,:] - tmp_f[:,np.newaxis] * b[1:,:] # b = [[ 6.][ 21.][ 3.]]

tmp_f1=A[1,1]/A[2:,1]
A[2:,:] = A[1,:] - tmp_f1[:,np.newaxis] * A[2:,:] # A = [[ 3.  2.  5.][ 0.  14.  14.][ 0.  0.  25.2]]
b[2:,:] = b[1,:] - tmp_f1[:,np.newaxis] * b[2:,:] # b = [[ 6.][ 21.][ 12.6]]

# 回代
x = np.zeros((3,1))
x[2,0] = b[2,0]/A[2,2] # x = [[ 0.][ 0.][ 0.5]]
x[1,0] = (b[1,0] - A[1,2]*x[2,0])/A[1,1] # x = [[ 0.][ 1.][ 0.5]]
x[0,0] = (b[0,0] - A[0,1]*x[1,0] - A[0,2]*x[2,0])/A[0,0] # x = [[ 0.5][ 1.][ 0.5]]
print(x)

[[0.5]
 [1. ]
 [0.5]]

print((-4)*(-2.8) + 14)
tmpa = (-4)*(-2.8) + 14
print(3*(-2.8) + 21)
tmpb = 3*(-2.8) + 21
print(tmpb/tmpa)

25.2
12.600000000000001
0.5000000000000001

所需前置numpy索引的知识

# demo 重新理解numpy的一些操作
import numpy as np
A = np.array([[3, 2, 5], [-1, 4, 3], [1, -1, 3]])

# 学习
print(A[0, 1]) # 第一个参数是行，第二个参数是列
print(A[0, 1].shape) # ()
print(A[0, 1].size) # 1
print(A[0, 1].ndim) # 0
print(A[0, 1].dtype) # int32
print(A[0, 1].itemsize) # 4
print(A[0, 1].nbytes) # 4

print(A[0:2, 0:2]) # 第一个参数是行，第二个参数是列，这里是取前两行，前两列 [[3 2] [-1 4]]
print(A[0:2, 0:2].shape) # shape是一个元组，表示矩阵的行数和列数 (2, 2)
print(A[0:2, 0:2].size) # size是一个整数，表示矩阵的元素个数 4
print(A[0:2, 0:2].ndim) # ndim是一个整数，表示矩阵的维数 2
print(A[0:2, 0:2].dtype) # dtype是一个对象，表示矩阵的数据类型 int32
print(A[0:2, 0:2].itemsize) # itemsize是一个整数，表示矩阵中每个元素的字节数 4
print(A[0:2, 0:2].nbytes) # nbytes是一个整数，表示矩阵中所有元素的字节数 16



2
()
0
1
int32
4
4
[[ 3  2]
 [-1  4]]
(2, 2)
4
2
int32
4
16

例2

\[\begin{cases} 3x_1 - x_2 + 4x_3 = 7 \\ -x_1 + 2x_2 - 2x_3 = -1 \\ 2x_1 - 3x_2 - 2x_3 = 0 \end{cases} \]

普通解法

代码如下：

import numpy as np
A = np.array([[3, -1, 4], [-1, 2, -2], [2, -3, -2]])
b = np.array([7, -1, 0])
x = np.linalg.solve(A, b)
print(x)

# 普通解法先确定结果
import numpy as np
A = np.array([[3, -1, 4], [-1, 2, -2], [2, -3, -2]])
b = np.array([7, -1, 0])
x = np.linalg.solve(A, b)
print(x)

[2.  1.  0.5]

# 高斯若儿当消元法验证
import numpy as np
A = np.array([[3, -1, 4], [-1, 2, -2], [2, -3, -2]], dtype=float)
b = np.array([[7, -1, 0]], dtype=float)
b = b.T
gaussElimination(A, b)

array([2. , 1. , 0.5])

高斯消元法与克劳特消元法（Gauss and Crout）

\[\begin{cases} 3x_1 - x_2 + 4x_3 = 7 \\ -x_1 + 2x_2 - 2x_3 = -1 \\ 2x_1 - 3x_2 - 2x_3 = 0 \end{cases} \]

相比于高斯消元法，克劳特消元法的基本思想其实也差不多，前者是保证L矩阵的对角线元素为1，而后者是保证U矩阵的对角线元素为1。

对于上面的例子，我们可以将其转化为增广矩阵：

\[\begin{bmatrix} 3 & -1 & 4 & 7 \\ -1 & 2 & -2 & -1 \\ 2 & -3 & -2 & 0 \end{bmatrix} \]

对于为什么要进行LU分解，一开始自己也不是很理解，但其实就是将方程组转化为两个方程组，然后分别求解，最后再合并，这样就可以避免矩阵的逆运算，从而提高计算效率。

这里权当复习整理数值分析里的知识，重新写一下LU分解的过程。

高斯消元法

我们实际上是将\(AX = B\)中的\(A\)分解为\(LU\)，其中\(U\)是一个上三角矩阵，\(L\)是一个下三角矩阵，\(L\)的对角线元素都是1，

\[A = LU \]

\[\begin{bmatrix} 3 & -1 & 4 \\ -1 & 2 & -2 \\ 2 & -3 & -2 \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 \\ -\frac{1}{3} & 1 & 0 \\ \frac{2}{3} & -1.4 & 1 \end{bmatrix} \begin{bmatrix} 3 & -1 & 4 \\ 0 & \frac{5}{3} & -\frac{2}{3} \\ 0 & 0 & -5.6 \end{bmatrix} \]

import numpy as np
A = np.array([[3, -1, 4], [-1, 2, -2], [2, -3, -2]], dtype=float)
b = np.array([[7, -1, 0]], dtype=float)
b = b.T

# 高斯消元法
def gauss(A, b,isprint = True):
    n = A.shape[0]
    L = np.zeros((n, n))
    U = np.zeros((n, n))
    for i in range(n):
        L[i, i] = 1
        for j in range(i, n):
            U[i, j] = A[i, j] - L[i, :i] @ U[:i, j] # @表示矩阵乘法
        for j in range(i, n):
            L[j, i] = (A[j, i] - L[j, :i] @ U[:i, i])/U[i, i]
    y = np.linalg.solve(L, b)
    x = np.linalg.solve(U, y)
    if isprint:
        # 输出解方程的过程
        print("LU分解：")
        print("A = LU")
        print("L = ")
        print(L)
        print("其中我们有：Ly = b，解得")
        print("y = ")
        print(y)
        print("U = ")
        print(U)
        print("其中我们有：Ux = y，解得")
        print("x = ")
        print(x)
    return x

x = gauss(A, b)

LU分解：
A = LU
L = 
[[ 1.          0.          0.        ]
 [-0.33333333  1.          0.        ]
 [ 0.66666667 -1.4         1.        ]]
其中我们有：Ly = b，解得
y = 
[[ 7.        ]
 [ 1.33333333]
 [-2.8       ]]
U = 
[[ 3.         -1.          4.        ]
 [ 0.          1.66666667 -0.66666667]
 [ 0.          0.         -5.6       ]]
其中我们有：Ux = y，解得
x = 
[[2. ]
 [1. ]
 [0.5]]

克劳特消元法

import numpy as np
A = np.array([[3, -1, 4], [-1, 2, -2], [2, -3, -2]], dtype=float)
b = np.array([[7, -1, 0]], dtype=float)
b = b.T

def crout(A, b,isprint = True):
    n = A.shape[0]
    L = np.zeros((n, n))
    U = np.zeros((n, n))
    for i in range(n):
        L[i, i] = A[i, i] - L[i, :i] @ U[:i, i]
        for j in range(i, n):
            U[i, j] = (A[i, j] - L[i, :i] @ U[:i, j])/L[i, i]
        for j in range(i, n):
            L[j, i] = (A[j, i] - L[j, :i] @ U[:i, i])
    y = np.linalg.solve(L, b)
    x = np.linalg.solve(U, y)
    if isprint:
        # 输出解方程的过程
        print("LU分解：")
        print("A = LU")
        print("L = ")
        print(L)
        print("其中我们有：Ly = b，解得")
        print("y = ")
        print(y)
        print("U = ")
        print(U)
        print("其中我们有：Ux = y，解得")
        print("x = ")
        print(x)
    return x

x = crout(A, b)

LU分解：
A = LU
L = 
[[ 3.          0.          0.        ]
 [-1.          1.66666667  0.        ]
 [ 2.         -2.33333333 -5.6       ]]
其中我们有：Ly = b，解得
y = 
[[2.33333333]
 [0.8       ]
 [0.5       ]]
U = 
[[ 1.         -0.33333333  1.33333333]
 [ 0.          1.         -0.4       ]
 [ 0.          0.          1.        ]]
其中我们有：Ux = y，解得
x = 
[[2. ]
 [1. ]
 [0.5]]

补充LU分解知识

LU分解

设有线性方程组

\[Ax = b \]

其中，\(A\)为\(n \times n\)的矩阵，\(b\)为\(n\)维列向量，\(x\)为\(n\)维列向量。举个具体例子：

\[\begin{cases} a_{11}x_1 + a_{12}x_2 + a_{13}x_3 = b_1 \\ a_{21}x_1 + a_{22}x_2 + a_{23}x_3 = b_2 \\ a_{31}x_1 + a_{32}x_2 + a_{33}x_3 = b_3 \end{cases} \]

其中，\(a_{ij}\)为矩阵\(A\)的元素，\(b_i\)为列向量\(b\)的元素，\(x_i\)为列向量\(x\)的元素。也即：

\[\begin{bmatrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} = \begin{bmatrix} b_1 \\ b_2 \\ b_3 \end{bmatrix} \]

首先，对于第一个方程，我们选取一个常数\(l_{11}\)遍除第一行，得到：

\[\frac{a_{11}}{l_{11}}x_1 + \frac{a_{12}}{l_{11}}x_2 + \frac{a_{13}}{l_{11}}x_3 = \frac{b_1}{l_{11}} \]

对于上式，我们令

\[\begin{bmatrix} \frac{a_{11}}{l_{11}} & \frac{a_{12}}{l_{11}} & \frac{a_{13}}{l_{11}} \\ \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} = \begin{bmatrix} \frac{b_1}{l_{11}} \\ \end{bmatrix} \\ \Rightarrow \begin{bmatrix} u_{11} & u_{12} & u_{13} \\ \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} = \begin{bmatrix} y_1 \\ \end{bmatrix} \]

其中，\(u_{ij}\)为矩阵\(U\)的元素，\(y_i\)为列向量\(y\)的元素(待会也会求解它)。

接下来，我们对第二个方程和第三个方程进行第一步的消元，令\(l_{21} = \frac{a_{21}}{u_{11}} \quad l_{31} = \frac{a_{31}}{u_{11}}\)

\[\begin{bmatrix} a_{21} & a_{22} & a_{23} & b_2 \\ a_{31} & a_{32} & a_{33} & b_3 \end{bmatrix} - \begin{bmatrix} l_{21} \\ l_{31} \end{bmatrix} \begin{bmatrix} u_{11} & u_{12} & u_{13} & y_1 \\ \end{bmatrix} = \begin{bmatrix} 0 & a_{22} - l_{21}u_{12} & a_{23} - l_{21}u_{13} & b_2 - l_{21}y_1 \\ 0 & a_{32} - l_{31}u_{12} & a_{33} - l_{31}u_{13} & b_3 - l_{31}y_1 \end{bmatrix} \]

注意：

  • 这里的\(y_1\)是减去上面求解的\(y_1\)即可
  • 另外上面可以看到一个很有趣的地方，就是l和u的下标，比如\(a_{23} - l_{21}u_{13}\)把\(l_{21}u_{13}\)中的1去掉我们就得到了\(a_{23}\)的下标，这刚好又是我们当前操作元素的位置；

对于等式右边的矩阵（我们重新写成方程的形式，下面对它继续进行消元工作），

\[\begin{cases} (a_{22} - l_{21}u_{12})x_2 + (a_{23} - l_{21}u_{13})x_3 = b_2 - l_{21}y_1 \\ (a_{32} - l_{31}u_{12})x_2 + (a_{33} - l_{31}u_{13})x_3 = b_3 - l_{31}y_1 \end{cases} \]

我们选取第二个常数\(l_{22}\)遍除第一行，得到：

\[\begin{bmatrix} 0 & \frac{a_{22} - l_{21}u_{12}}{l_{22}} & \frac{a_{23} - l_{21}u_{13}}{l_{22}}\\ \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} = \begin{bmatrix} \frac{b_2 - l_{21}y_1}{l_{22}} \\ \end{bmatrix} \\ \Rightarrow \begin{bmatrix} 0 & u_{22} & u_{23} \\ \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} = \begin{bmatrix} y_2 \\ \end{bmatrix} \]

继续对第三个方程进行消元，令\(l_{32} = \frac{a_{32} - l_{31}u_{12}}{u_{22}}\)

\[\begin{align} &\begin{bmatrix} 0 & a_{32} - l_{31}u_{12} & a_{33} - l_{31}u_{13} & b_3 - l_{31}y_1 \end{bmatrix} - \begin{bmatrix} l_{32} \end{bmatrix} \begin{bmatrix} 0 & u_{22} & u_{23} & y_2 \end{bmatrix}\\ = &\begin{bmatrix} 0 & 0 & a_{33} - l_{31}u_{13} - l_{32}u_{23} & b_3 - l_{31}y_1 - l_{32}y_2 \end{bmatrix} \end{align} \]

对于等式右边的矩阵，我们选取第三个常数\(l_{33}\)遍除第一行，最终得到：

\[\begin{bmatrix} 0 & 0 & \frac{a_{33} - l_{31}u_{13} - l_{32}u_{23}}{l_{33}}\\ \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} = \begin{bmatrix} \frac{b_3 - l_{31}y_1 - l_{32}y_2}{l_{33}} \\ \end{bmatrix} \\ \Rightarrow \begin{bmatrix} 0 & 0 & u_{33} \\ \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} = \begin{bmatrix} y_3 \\ \end{bmatrix} \]

至此，我们按照消元的一定顺序步骤，得到了\(u\)和\(l\)矩阵，也得到了如下的方程组（矩阵形式）：

\[\begin{bmatrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} = \begin{bmatrix} b_1 \\ b_2 \\ b_3 \end{bmatrix} \\ \Rightarrow \begin{bmatrix} l_{11} & 0 & 0 \\ l_{21} & l_{22} & 0 \\ l_{31} & l_{32} & l_{33} \end{bmatrix} \begin{bmatrix} u_{11} & u_{12} & u_{13} \\ 0 & u_{22} & u_{23} \\ 0 & 0 & u_{33} \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} = \begin{bmatrix} b_1 \\ b_2 \\ b_3 \end{bmatrix} \\ \Rightarrow \begin{bmatrix} l_{11} & 0 & 0 \\ l_{21} & l_{22} & 0 \\ l_{31} & l_{32} & l_{33} \end{bmatrix} \begin{bmatrix} y_1 \\ y_2 \\ y_3 \end{bmatrix} = \begin{bmatrix} b_1 \\ b_2 \\ b_3 \end{bmatrix} \quad \begin{bmatrix} u_{11} & u_{12} & u_{13} \\ 0 & u_{22} & u_{23} \\ 0 & 0 & u_{33} \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} = \begin{bmatrix} y_1 \\ y_2 \\ y_3 \end{bmatrix} \\ \]

简洁起见，我们统一将\(u\)矩阵记为\(U\)，将\(l\)矩阵记为\(L\)，则上述解方程的过程可以写为：

\[Ax = b \\ \Rightarrow LUx = b \\ \Rightarrow Ly = b \quad \text{（回代法解出y）} \\ \Rightarrow Ux = y \quad \text{（回代法解出x）} \]

代码实现

import numpy as np

def lu_decomposition(A):
    """
    LU分解
    :param A: 矩阵
    :return: L, U
    """
    n = A.shape[0]
    L = np.zeros((n, n))
    U = np.zeros((n, n))
    for i in range(n):
        for j in range(i, n):
            U[i, j] = A[i, j] - np.dot(L[i, :i], U[:i, j])
        for j in range(i, n):
            if i == j:
                L[i, i] = 1 # 这里的选定将决定是高斯消元法还是克劳特消元法
            else:
                L[j, i] = (A[j, i] - np.dot(L[j, :i], U[:i, i])) / U[i, i]
    return L, U

def forward_substitution(L, b):
    """
    前代法求解Ly = b
    :param L: 下三角矩阵
    :param b: 常数向量
    :return: y
    """
    n = L.shape[0]
    y = np.zeros(n)
    for i in range(n):
        y[i] = (b[i] - np.dot(L[i, :i], y[:i])) / L[i, i]
    return y

def backward_substitution(U, y):
    """
    回代法求解Ux = y
    :param U: 上三角矩阵
    :param y: 常数向量
    :return: x
    """
    n = U.shape[0]
    x = np.zeros(n)
    for i in range(n - 1, -1, -1):
        x[i] = (y[i] - np.dot(U[i, i + 1:], x[i + 1:])) / U[i, i]
    return x

def solve(A, b):
    """
    求解Ax = b
    :param A: 矩阵
    :param b: 常数向量
    :return: x
    """
    L, U = lu_decomposition(A)
    y = forward_substitution(L, b)
    x = backward_substitution(U, y)
    return x

测试代码

A = np.array([[3, -1, 4], [-1, 2, -2], [2, -3, -2]], dtype=float)
b = np.array([7, -1, 0], dtype=float) # 注意：这里的b是一维数组，如果要写成矩阵形式即二维向量，需要转置一下

lu_decomposition(A)

(array([[ 1.        ,  0.        ,  0.        ],
        [-0.33333333,  1.        ,  0.        ],
        [ 0.66666667, -1.4       ,  1.        ]]),
 array([[ 3.        , -1.        ,  4.        ],
        [ 0.        ,  1.66666667, -0.66666667],
        [ 0.        ,  0.        , -5.6       ]]))

# b = b.T
x = solve(A, b)
print(x)

[2.  1.  0.5]

A = np.array([[3, 2, 5], [-1, 4, 3], [1, -1, 3]])
b = np.array([6, 5, 1])
L,U = lu_decomposition(A)
print(L)
print(U)
y = forward_substitution(L, b)
print(y)
x = backward_substitution(U, y)
print(x)

[[ 1.          0.          0.        ]
 [-0.33333333  1.          0.        ]
 [ 0.33333333 -0.35714286  1.        ]]
[[3.         2.         5.        ]
 [0.         4.66666667 4.66666667]
 [0.         0.         3.        ]]
[6.  7.  1.5]
[0.5 1.  0.5]

-5/14

-0.35714285714285715

小结

LU分解

总而言之，LU分解的过程就是将方程组\(Ax = b\)转化为\(LUx = b\)的过程，其中\(U\)是上三角矩阵，\(L\)是下三角矩阵。

我们逐步完成这个过程的同时也逐渐实现了消元的目的，整个过程简洁如下：

  • 自己选定\(l_{ii}\)这个常数，若 \(l_{ii} = 1\)那么我们就是用高斯消元法，若\(l_{ii} = a_{ii} - \sum_{k=1}^{i-1}l_{ik}u_{ki}\)则我们就是用克劳特消元法，

  • 其他的常数都是通过消元得到的，如：

\[\begin{align} l_{ij} &= \frac{a_{ij} - \sum_{k=1}^{i-1}l_{ik}u_{kj}}{u_{ii}}\\ u_{ij} &= \frac{a_{ij} - \sum_{k=1}^{i-1}l_{ik}u_{kj}}{l_{ii}}\\ \end{align} \]

LU分解的主要部分写成代码就是这七行：

for i in range(n):
    # L[i, i] = A[i, i] - L[i, :i] @ U[:i, i] 克劳特消元法
    L[i, i] = 1 # 高斯消元法
    for j in range(i, n):
        U[i, j] = (A[i, j] - L[i, :i] @ U[:i, j])/L[i, i]
    for j in range(i, n):
        L[j, i] = (A[j, i] - L[j, :i] @ U[:i, i])/U[i, i]
回代法

回代其实就非常简单了

\[\begin{align} y_i &= \frac{b_i - \sum_{k=1}^{i-1}l_{ik}y_k}{l_{ii}}\\ x_i &= \frac{y_i - \sum_{k=i+1}^{n}u_{ik}x_k}{u_{ii}}\\ \end{align} \]

回代法的主要部分写成代码就是这四行：

for i in range(n):
    y[i] = (b[i] - L[i, :i] @ y[:i])/L[i, i]
for i in range(n-1, -1, -1):
    x[i] = (y[i] - U[i, i+1:] @ x[i+1:])/U[i, i]

注：上面出现过的@符号是矩阵乘法的符号，是python3.5之后才支持的，如果你的python版本低于3.5，可以使用np.dot()代替。

例3

\[\begin{cases} 3x_1 - x_2 + 2x_3 = 7 \\ -x_1 + 2x_2 - 2x_3 = -1 \\ 2x_1 - 2x_2 + 4x_3 = 0 \end{cases} \]

普通解法

import numpy as np
A = np.array([
    [3, -1, 2], 
    [-1, 2, -2], 
    [2, -2, 4]
    ])
b = np.array([
    7, -1, 0
    ])
x = np.linalg.solve(A, b)
print(x)

[ 3.5  -1.   -2.25]

平方根法(Cholesky decomposition)

当A是一个对称正定矩阵，我们在选取\(l_{ii}\)的值时就可以做一些对称的操作，比如直接令\(l_{ii} = u_{ii}\)，这样就可以减少一些计算量。

在平方根法中，我们选取\(l_{ii} = \sqrt{a_{ii} - \sum_{k=1}^{i-1}l_{ik}^2}\)作为\(l_{ii}\)，这样我们就可以将\(A\)分解为\(LL^T\)，其中\(U = L^T\)。

\[\begin{align} A &= LL^T\\ A &= \begin{bmatrix} l_{11} & 0 & 0 \\ l_{21} & l_{22} & 0 \\ l_{31} & l_{32} & l_{33} \end{bmatrix} \begin{bmatrix} l_{11} & l_{21} & l_{31} \\ 0 & l_{22} & l_{32} \\ 0 & 0 & l_{33} \end{bmatrix}\\ &= \begin{bmatrix} l_{11}^2 & l_{11}l_{21} & l_{11}l_{31} \\ l_{21}l_{11} & l_{21}^2 + l_{22}^2 & l_{21}l_{31} + l_{22}l_{32} \\ l_{31}l_{11} & l_{31}l_{21} + l_{32}l_{22} & l_{31}^2 + l_{32}^2 + l_{33}^2 \end{bmatrix}\\ \end{align} \]

代码实现

def cholesky(A):
    n = A.shape[0]
    L = np.zeros((n, n))
    for i in range(n):
        L[i, i] = np.sqrt(A[i, i] - L[i, :i] @ L[i, :i])
        for j in range(i+1, n):
            L[j, i] = (A[j, i] - L[j, :i] @ L[i, :i])/L[i, i]
    return L

测试代码

A = np.array([
    [3, -1, 2], 
    [-1, 2, -2], 
    [2, -2, 4]
    ])
cholesky(A)

array([[ 1.73205081,  0.        ,  0.        ],
       [-0.57735027,  1.29099445,  0.        ],
       [ 1.15470054, -1.03279556,  1.26491106]])

b = np.array([
    7, -1, 0
    ])
L = cholesky(A)
y = np.linalg.solve(L, b)
print(y)
x = np.linalg.solve(L.T, y)
print(x)

[ 4.04145188  1.03279556 -2.84604989]
[ 3.5  -1.   -2.25]

补充numpy中矩阵的一系列自带函数

import numpy as np
A = np.array([
    [3, -1, 2], 
    [-1, 2, -2], 
    [2, -2, 4]
    ])
# 求A的转置
print(A.T)
'''
[[ 3 -1  2]
 [-1  2 -2]
 [ 2 -2  4]]
 '''
# 求A的逆矩阵 
A_inv = np.linalg.inv(A)
print(A_inv)
'''
[[ 0.5    0.    -0.25 ]
 [ 0.     1.     0.5  ]
 [-0.25   0.5    0.625]]
'''
# 求A的行列式
A_det = np.linalg.det(A)
print(A_det)
'''
8.000000000000002
'''
# 求A的特征值和特征向量
A_eig = np.linalg.eig(A)
print(A_eig)
'''
(array([6.61185871, 1.65867531, 0.72946598]), 
array([[ 0.52548211, -0.83393758, -0.16857242],
       [-0.43184431, -0.43214537,  0.7916823 ],
       [ 0.73306142,  0.34321785,  0.58721586]]))
'''
# 求A的秩
A_rank = np.linalg.matrix_rank(A)
print(A_rank) # 3
# 求A的范数
A_norm = np.linalg.norm(A)
print(A_norm) # 6.855654600401044
# 求A的条件数
A_cond = np.linalg.cond(A)
print(A_cond) # 9.063971310181588
# 求A的奇异值分解
A_svd = np.linalg.svd(A)
print(A_svd)
'''
(array([
    [-0.52548211, -0.83393758, -0.16857242],
    [ 0.43184431, -0.43214537,  0.7916823 ],
    [-0.73306142,  0.34321785,  0.58721586]
    ]), 
 array([6.61185871, 1.65867531, 0.72946598]), 
 array([
    [-0.52548211,  0.43184431, -0.73306142],
    [-0.83393758, -0.43214537,  0.34321785],
    [-0.16857242,  0.7916823 ,  0.58721586]
    ]))
'''
# 求A的QR分解
A_qr = np.linalg.qr(A)
print(A_qr)
'''
(
array([[-0.80178373, -0.5179324 , -0.2981424 ],
       [ 0.26726124, -0.75697812,  0.59628479],
       [-0.53452248,  0.39840954,  0.74535599]]), 
array([[-3.74165739,  2.40535118, -4.27617987],
       [ 0.        , -1.79284291,  2.07172959],
       [ 0.        ,  0.        ,  1.19256959]])
)
'''

总结

在数值分析这门课中的解线性方程组的直接法，个人认为掌握LU分解一个点就够了；平方根法的话，如果你的矩阵是对称正定的，那么可以考虑使用平方根法，因为它的计算量比LU分解小很多。

此外，我们还学习了其它的解方程的直接法，比如列主元素法，全主元素法，个人觉得这些方法都是在高斯消元法的基础上做了一些优化，所以我就不再赘述了。（因为优化程度不太大，只是选取一个最优的参数来给\(l_{ii}\)操作而已）

还有针对特殊的矩阵，比如对称矩阵，我们还可以使用平方根法，这个方法在上面已经有所介绍了。针对三对角矩阵，我们可以使用Thomas算法（课本写作追赶法），但还是那句话，都是在高斯消元法的基础上做了一些优化而已；

下一章我们将学习迭代法来求解线性方程组。

一些测试

A = np.array([
    [2, -1, 0, 0],
    [-1, 2, -1, 0],
    [0, -1, 2, -1],
    [0, 0, -1, -2]
    ])
b = np.array([
    0, 1, 0, 2.5
    ])
# 求解线性方程组
x = np.linalg.solve(A, b)
print(x) 

[ 0.22727273  0.45454545 -0.31818182 -1.09090909]

lu_decomposition(A)

(array([[ 1.        ,  0.        ,  0.        ,  0.        ],
        [-0.5       ,  1.        ,  0.        ,  0.        ],
        [ 0.        , -0.66666667,  1.        ,  0.        ],
        [ 0.        ,  0.        , -0.75      ,  1.        ]]),
 array([[ 2.        , -1.        ,  0.        ,  0.        ],
        [ 0.        ,  1.5       , -1.        ,  0.        ],
        [ 0.        ,  0.        ,  1.33333333, -1.        ],
        [ 0.        ,  0.        ,  0.        , -2.75      ]]))

tmp = np.array([
    [2, -1, 0, 0],
    [0,1.5, -1, 0],
    [0, 0, 4/3, -1],
    [0, 0, 0, -2.75]
    ])
tmpy = np.array([
    0, 1, 2/3, 3
    ])
np.linalg.solve(tmp, tmpy)

array([ 0.22727273,  0.45454545, -0.31818182, -1.09090909])

标签：begin,end,线性方程组,数值,print,线性代数,bmatrix,np,array
From： https://www.cnblogs.com/Linkdom/p/16749185.html

相关文章

  • wsprintf 将一系列的字符和数值输入到缓冲区
    wsprintf将一系列的字符和数值输入到缓冲区函数wsprintf[1] ()将一系列的字符和数值输入到缓冲区。输出缓冲区里的的值取决于格式说明符（即"%"）。如果写入的是文字，此函......
  • 函数式编程（纯函数、避免改变参数值和全局变量）、使用.map()、.filter()使代码更简洁
    函数式编程是：独立函数——不依赖于程序（里面含有可能更改的全局变量）的状态。纯函数——相同的输入总是产生相同的输出，不产生副作用（不更改参数值和全局变量值）。副作用尽量小......
  • Sass 变量使用、数值运算以及@import 导入
    前言本文主要记录了Sass基本用法，包含变量的使用、简单的数值运算、@import导入Sass文件以及Sass@import与CSS@import的区别。变量使用在Sass中，支持使用变量，变量......
  • JS奇淫技巧：数值的七种写法
    JS奇淫技巧：数值的七种写法JS奇淫技巧：挑战前端黑科技，数值的七种写法，能全看懂的一定是高手你知道吗？在JS编程中，数值可以有很多种写法。第一种写法：一般情况而言，数值就是数值。比......
  • 算法数学笔记-二、线性代数
    目录二、线性代数矩阵模板高斯消元二、线性代数矩阵模板namespaceMatrix{ structmatrix{ inthang,lie; vector<vector<int>>num; matrix(intx=0,......
  • .NET教程 - 数值 & 随机数（Number & Random）
    更新记录转载请注明出处：2022年9月29日发布。2022年9月28日从笔记迁移到博客。System.Numerics.BigIntegerBigInteger说明BigInteger类型用于表示任意大的整数数......
  • 对数值类网游反外挂的一些思路（想法）
    今天在某群聊的时候,看见有人在谈论某网游的外挂..在有一个人问是“什么挂”的时候，我却想到另一个问题：我要是这公司的人，我该怎么反外挂？几乎是一瞬间的，我就想到2个方法1......
  • python 数值类型-浮点数，复数。数值之间的运算
    1.浮点数1.1python中的浮点数和C语言一样有误差`0.1+0.20.300000000000000040.3==0.1+0.2False`那如何精确的计算浮点数？需要借助decimal(十进制的)模块`impor......
  • python 线性代数：解多元一次方程
    因为在程序化交易策略中使用了网格算法进行交易，因为在网格中想设置动态资源大小的问题，所以就想到使用抛物线的分布方法来对网格资金配置进行分配。比如我的网格最大值设置......
  • 数值计算(1) --求解连续微分系统和混沌系统
    前言微分系统在工程项目中很常见，通过物理建模之后，基本都需要求解微分方程得到其结果，混沌系统属于特殊的一类微分系统，在某些项目上也很常见，同时可以引申出分岔图、李雅普诺夫......

赞助商

阅读排行

网站主页关于我们联系我们网站地图

本网站内容转载自其他媒体，侵权联系[admin##ips99.com]。

Copyright © 2020-2023 IPS99 版权所有 IPS99