【数值计算】追赶法

有心友兴

已于 2022-11-12 19:26:59 修改

阅读量8.2k

点赞数 5

分类专栏：【数学基础】文章标签：算法矩阵线性代数

于 2022-08-10 22:43:38 首次发布

本文链接：https://blog.csdn.net/yxyx13120297/article/details/126275562

版权

【数学基础】专栏收录该内容

1 篇文章 0 订阅

订阅专栏

记录一下追赶法解特殊方程组，设有 $n$ 阶方程组 $\mathbf{A}x=d$ ，其中 $\mathbf{A}$ 为三对角矩阵，即
$\mathbf{A} = \left[ \begin{array}{ccc} 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], \quad d = \left[ \begin{array}{ccc} d_{1}\\ d_{2}\\ \vdots\\ d_{n-1}\\ d_{n}\\ \end{array} \right]$
对矩阵 $\mathbf{A}$ 进行克洛托分解（ $L U$ 分解），得到
$\mathbf{L} = \left[ \begin{array}{ccc} 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[ \begin{array}{ccc} 1 & u_{1}\\ & 1 & u_{2}\\ & & \ddots & \ddots\\ & & & 1 & u_{n-1}\\ & & & & 1 & u_{n}\\ \end{array} \right]$
因此 $\mathbf{A}x=d$ 等价于求解 $L y = d$ 和 $U x = y$ ，于是可以推得追赶法算法过程如下:
1.计算 $l_{i}$ ， $u_{i}$ 递推公式。
$l_{1} = b_{1}, \quad l_{i} = b_{i} - a_{i}u_{i-1}, \quad i=2,3,4,\dots,n\\ u_{1} = c_{1} / b_{1},\quad u_{i} = c_{i} / l{i},\quad\quad i=2,3,4,\dots,n$
2.求解 $L y = d$ 。
$y_{1} = d_{1}/b_{1}, \quad y_{i} = (d_{i} - a_{i}y_{i-1})/l_{i}, \quad i=2,3,4,\dots,n$
3.求解 $U x = y$ 。
$x_{n} = y_{n}, \quad x_{i} = y_{i} - u_{i}x_{i+1}, \quad 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

追赶法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;
}