追赶法实现

1.追赶法原理

解系数矩阵为对角占优的三对角线方程组
$$\left[\begin{array}{ccccc}{b}_1& c_1& & & \\ a_2& b_2& c_2& & \\ & \ddots & \ddots & \ddots & \\ & & a_{n-1}& b_{n-1}& c_{n-1}\\ & & & a_n& b_n\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_{n-1}\\ x_n\end{array}\right]=\left[\begin{array}{c}{f}_1\\ f_2\\ ⋮\\ f_{n-1}\\ f_n\end{array}\right],$$
简记为$Ax=f.$其中，当$|i-j|>1$时，$a_{ij}=1,$且:
$|b_1|>|c_1|>0;$
$|b_i|\ge |a_i|+|c_i|,a_i,c_i\ne 0,i=2,3,\dots ,n-1$
$|b_n|>|a_n|>0.$
设$$A=\left[\begin{array}{ccccc}{b}_1& c_1& & & \\ a_2& b_2& c_2& & \\ & \ddots & \ddots & \ddots & \\ & & a_{n-1}& b_{n-1}& c_{n-1}\\ & & & a_n& b_n\end{array}\right]\left[\begin{array}{cccc}{\alpha }_1& & & \\ {\gamma }_2& {\alpha }_2& & \\ & \ddots & \ddots & \\ & & {\gamma }_n& {\alpha }_n\end{array}\right]=\left[\begin{array}{cccc}1& {\beta }_1& & \\ & 1& \ddots & \\ & & \ddots & {\beta }_{n-1}\\ & & & 1\end{array}\right],$$其中${\alpha }_i,{\beta }_i,{\gamma }_i为待定系数.$
经过分析，求解$Ax=f$等价于求解两个三角形方程组:
$(1)Ly=f,求y;$ $(2)Ux=y,求x.$
从而得到解三对角线方程组的追赶法公式:
(1)计算{${\beta }_i$}的递推公式
$$\begin{array}{rl}{\beta }_1& =c_1/b_1,\\ {\beta }_i& =c_i/(b_i-{\alpha }_i{\beta }_{i-1}),i=2,3,\dots ,n-1;\end{array}$$
(2)解$Ly=f$
$$\begin{array}{rl}{y}_1& =f_1/b_1,\\ y_i& =(f_i-{\alpha }_i{y}_{i-1})/(b_i-{\alpha }_i{\beta }_{i-1}),i=2,3,\dots ,n;\end{array}$$
(3)解$Ux=y$
$$\begin{array}{rl}{x}_n& =y_n,\\ x_i& =y_i-{\beta }_i{x}_{i+1},i=n-1,n-2,\dots ,1.\end{array}$$

2.Python实现追赶法解三对角线方程组

# 自己原创，注意这里没有做矩阵类型判断，该方法仅适用于解系数矩阵为对角占优的三对角线方程组,P159
def chasing_thomas(triple_diagonal_matrix: np.ndarray, f_vector: np.ndarray):
    b_diagonal = np.diag(triple_diagonal_matrix)
    c_diagonal = np.diag(triple_diagonal_matrix, k=1)
    a_diagonal = np.diag(triple_diagonal_matrix, k=-1)
    # 计算Bata_i,dtype=np.float64,提高运算精度，否则默认最小精度
    beta_diagonal = np.ones_like(c_diagonal, dtype=np.float64)
    beta_diagonal[0] = c_diagonal[0] / b_diagonal[0]
    rows = triple_diagonal_matrix.shape[0]
    for i in range(1, rows - 1):
        beta_diagonal[i] = c_diagonal[i] / (b_diagonal[i] - a_diagonal[i] * beta_diagonal[i - 1])
    # 解Ly=f
    y = np.ones_like(f_vector, dtype=np.float64)
    y[0, 0] = f_vector[0, 0] / b_diagonal[0]
    for i in range(1, rows):
        y[i, 0] = (f_vector[i, 0] - a_diagonal[i - 1] * y[i - 1, 0]) / (
                b_diagonal[i] - a_diagonal[i - 1] * beta_diagonal[i - 1])
    # 解Ux=y,自己写的发现第一个解x0=1.,而不是0.8333333333333334，有精度损失，
    x = np.ones_like(f_vector, dtype=np.float64)
    x[rows - 1, 0] = y[rows - 1, 0]
    for i in range(rows - 2, 0, -1):
        x[i, 0] = y[i, 0] - beta_diagonal[i] * x[i + 1, 0]
    return x

3.解三对角线方程

$$\left[\begin{array}{ccccc}2& -1& 0& 0& 0\\ -1& 2& -1& 0& 0\\ 0& -1& 2& -1& 0\\ 0& 0& -1& 2& -1\\ 0& 0& 0& -1& 2\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\\ x_4\\ x_5\end{array}\right]=\left[\begin{array}{c}1\\ 0\\ 0\\ 0\\ 0\end{array}\right]$$

4.测试

if __name__ == '__main__':
	# 追赶法测试成功，来源详见李庆扬数值分析第5版P177，e.g.9
    # dtype指定为双精度浮点提高运算精度
    tri_diagonal = np.array([[2, -1, 0, 0, 0],
                             [-1, 2, -1, 0, 0],
                             [0, -1, 2, -1, 0],
                             [0, 0, -1, 2, -1],
                             [0, 0, 0, -1, 2]], dtype=np.float64)
    f = np.array([1, 0, 0, 0, 0], dtype=np.float64).reshape((5, 1))
    print(chasing_thomas(tri_diagonal, f))

5.运行截图

说明：解向量{$x_i$}为倒序$x_5,\dots ,x_1.$
其中$x_1=\frac{5}{6}$,由于自己实现方法精度不足，故只输出保留到整数位