【数值计算】追赶法

   记录一下追赶法解特殊方程组，设有$n$阶方程组$\mathbf{A}x=d$，其中$\mathbf{A}$为三对角矩阵，即
A = [ b 1 c 1 a 2 b 2 c 2 ⋱ ⋱ ⋱ a n − 1 b n − 1 c n − 1 a n b n ] , d = [ d 1 d 2 ⋮ d n − 1 d n ] \mathbf{A} = \left[

\right], \quad d = \left[ \right] A=⎣ ⎡​b1​a2​​c1​b2​⋱​c2​⋱an−1​​⋱bn−1​an​​cn−1​bn​​⎦ ⎤​,d=⎣ ⎡​d1​d2​⋮dn−1​dn​​⎦ ⎤​
   对矩阵$\mathbf{A}$进行克洛托分解（$LU$分解），得到
L = [ l 1 v 2 l 2 ⋱ ⋱ v n − 1 l n − 1 v n l n ] , U = [ 1 u 1 1 u 2 ⋱ ⋱ 1 u n − 1 1 u n ] \mathbf{L} = \left[

l1v2l2⋱⋱vn−1ln−1vnln" role="presentation" style="position: relative;">l1v2l2⋱⋱vn−1ln−1vnln$$\begin{array}{cc}{l}_1& \\ v_2& l_2\\ & \ddots & \ddots \\ & & v_{n-1}& l_{n-1}\\ & & & v_n& l_n\end{array}$$

\right], \quad \mathbf{U} = \left[

1u11u2⋱⋱1un−11un" role="presentation" style="position: relative;">1u11u2⋱⋱1un−11un$$\begin{array}{cc}1& u_1\\ & 1& u_2\\ & & \ddots & \ddots \\ & & & 1& u_{n-1}\\ & & & & 1& u_n\end{array}$$

\right] L=⎣ ⎡​l1​v2​​l2​⋱​⋱vn−1​​ln−1​vn​​ln​​⎦ ⎤​,U=⎣ ⎡​1​u1​1​u2​⋱​⋱1​un−1​1​un​​⎦ ⎤​
   因此$\mathbf{A}x=d$等价于求解$Ly=d$和$Ux=y$，于是可以推得追赶法算法过程如下:
   1.计算$l_i$，$u_i$递推公式。
$$\begin{gathered} l_1=b_1,l_i=b_i-a_i{u}_{i-1},i=2,3,4,\dots ,n \\ {u}_1=c_1/b_1,u_i=c_i/li,i=2,3,4,\dots ,n \end{gathered}$$
   2.求解$Ly=d$。
$$y_1=d_1/b_1,y_i=(d_i-a_i{y}_{i-1})/l_i,i=2,3,4,\dots ,n$$
   3.求解$Ux=y$。
$$x_n=y_n,x_i=y_i-u_i{x}_{i+1},i=n-2,n-2,\dots ,1$$
   故可将上述算法过程转换成代码：

追赶法MATLAB程序代码

function x = tridiagsolver(a, b, c, d)
n = length(b);
l(1) = b(1);
y(1) = d(1) / l(1);
u(1) = c(1) / l(1);
for ii = 2 : n
    l(ii) = b(ii) - a(ii - 1) * u(ii - 1);
    y(ii) = (d(ii) - y(ii - 1) * a(ii - 1)) / l(ii);
    if ii <= n-1
    	u(ii) = c(ii) / l(ii);
    end
end

x(n) = y(n);
for ii = n-1 : -1 : 1
    x(ii) = y(ii) - u(ii) * x(ii + 1);
end
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

追赶法C++程序代码

/*************************************************************************
 void tridiagsolver(double *dst, double *a, double *b, double *c, double *d, int num)
 功能： 追赶法求解三角矩阵
 参数： 
        dst：   输出结果
        a:      稀疏矩阵最下边部分的元素
        b：     稀疏矩阵中间部分的元素
        c：     稀疏矩阵最上边部分的元素
        d:      右端向量
        num:    矩阵大小
返回值： 无
 *************************************************************************/
void tridiagsolver(double *dst, double *a, double *b, double *c, double *d, int num)
{
    double *l = new double[num];
    double *y = new double[num];
    double *u = new double[num - 1];

    // 追过程
    l[0] = b[0];
    y[0] = d[0] / l[0];
    u[0] = c[0] / l[0];
    for (int i = 1; i < num; i++)
    {
        l[i] = b[i] - a[i-1] * u[i-1];
        y[i] = (d[i] - y[i-1] * a[i-1]) / l[i];
        if (i < num - 1)
                u[i] = c[i] / l[i];
    }

    // 赶过程
    dst[num-1] = y[num-1];
    for (int i = num-2; i >= 0; i--)
    {
        dst[i] = y[i] - u[i] * dst[i+1];
    }

    delete l;
    delete y;
    delete u;
}
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

   需要注意的是，程序中的代码和算法求解过程存在略微差异，是由于存储对三角a，b，c的数组大小不一致引起。如果输入的数组大小一致，则需要对代码进行相应修改。如图1所示的是对三角方程的求解过程。

图1 求解结果示意

参考文献：
   [1] 现代数值计算/同济大学计算数学教研室编著.——2版. ——北京：人民邮电出版社，2021.12冲印