三对角方程组的求解(C++)

本文我们讨论一类特殊的线性方程组, 即三对角方程组

$$\left(\begin{array}{cccccc}{b}_1& c_1& 0& 0& \dots & 0\\ a_2& b_2& c_2& 0& \dots & 0\\ 0& a_3& b_3& c_3& \dots & 0\\ ?& ?& ?& ?& ?& ?\\ 0& 0& 0& 0& a_n& b_n\end{array}\right)\left(\begin{array}{c}{x}_1\\ x_2\\ x_3\\ ?\\ x_n\end{array}\right)=\left(\begin{array}{c}{d}_1\\ d_2\\ d_3\\ ?\\ d_n\end{array}\right)$$

其中, $a_i,b_i,c_i$ 是已知系数, $x_i$ 是未知变量, $d_i$ 是常数项.

三对角方程组是一类特殊的线性方程组, 其系数矩阵具有三对角的性质, 即除了主对角线上的元素外, 只有上对角线和下对角线上的非零元素. 这种形式的方程组在科学和工程领域中经常出现, 因此求解三对角方程组具有重要的应用价值.

求解三对角方程组的方法有多种, 其中一种常用的方法是追赶法(也称为托马斯算法). 追赶法通过一系列前向和后向的替代操作, 将三对角方程组转化为一个更容易求解的上三角或下三角方程组, 然后使用回代或前代方法求解. 这个方法通常具有较高的数值稳定性和效率. 以下是追赶法的主要步骤:

1. 前向消元
   1. 追赶法首先通过前向消元来将原始三对角方程组转化为上三角形式. 在这一步, 通过逐行操作, 将非零元素下移到主对角线下, 同时更新右侧常数项.
   2. 这一步骤通常从第一行开始, 通过一系列线性组合将上对角线上的元素归零. 然后, 将结果用于更新下一行, 依此类推, 直到将原方程组转化为上三角形式.
2. 后向代入(或回代) 一旦方程组被转化为上三角形式, 就可以使用回代或后向代入的方法来求解方程组. 这是一个从最后一行开始的过程, 逐步计算出每个未知变量的值, 通过从下到上的线性组合.

算法可表示为伪代码如下:

其C++实现如下:

#include <armadillo>
using namespace arma;
#define THROW_DIV_0 throw "方阵为奇异矩阵！";
#define THROW_DIV_0_DEL \
    {                   \
        free(B);        \
        THROW_DIV_0     \
    }
// #define THROW_DIV_0_ALL {free(A);free(B);free(C);THROW_DIV_0}
#define JUDGE_0(x) \
    if (!(x))      \
    THROW_DIV_0
#define JUDGE_0_DEL(x) \
    if (x)             \
    THROW_DIV_0_DEL
// #define JUDGE_0_ALL(x) if(x)THROW_DIV_0_ALL
// 追赶过程, 要求b[0]不为零, n>1
// 结果输出在d中
bool chase(unsigned n, const double *a, double *b, const double *c, double *d)
{
    unsigned k(n);
    while (--k)
    {
        if (!*b)
            return true;
        double m(*a++ / *b), &pre_d(*d++);
        *++b -= m * *c++;
        *d -= m * pre_d;
    }
    if (!*b)
        return true;
    *d /= *b;
    while (--n)
    {
        double &pre_d(*d--);
        (*d -= *--c * pre_d) /= *--b;
    }
    return false;
}
// b[0]为零的追赶, n>2
bool chase_0(unsigned n, const double *a, double *b, const double *c, double *d)
{
    if (!*c || !*a)
        return true;
    const double t(*d / *c++);
    *++d -= t * *(++b)++;
    *++d -= t * *++a;
    unsigned k(n);
    while (--k)
    {
        if (!*b)
            return true;
        const double m(*++a / *b), &pre_d(*d++);
        *++b -= m * *++c;
        *d -= m * pre_d;
    }
    if (!*b)
        return true;
    *d /= *b;
    a -= n;
    while (--n)
    {
        const double &pre_x(*d--);
        (*d -= *c-- * pre_x) /= *--b;
    }
    const double &x3(*d);
    double &x2(*--d);
    *--d = (x2 - *c * x3) / *a;
    x2 = t;
    return false;
}
/*
 * 追赶法
 * R:解向量
 * a:左下方的副对角线
 * b:主对角线
 * c:右上方的副对角线
 * d:常数项
 */
void Thomas_algorithm(vec &R, const vec &a, const vec &b, const vec &c, const vec &d)
{
    if (b.n_elem != d.n_elem)
        throw "系数矩阵阶数与常数向量维数不相等！";
    if (a.n_elem != b.n_elem - 1 || a.n_elem != c.n_elem)
        throw "主、副对角线维数不匹配！";
    if (a.empty())
    {
        double t(d.at(0) / b.at(0));
        JUDGE_0(t)
        R = {t};
        return;
    }
    unsigned n(b.n_elem);
    R.set_size(n);
    const double *Bp(&b.at(--n)), *Dp(&d.at(n));
    double *B((double *)malloc(sizeof(double) * b.n_elem)), *bp(B + n), *dp(&R.at(n));
    *bp = *Bp, *dp = *Dp;
    do
        *--bp = *--Bp, *--dp = *--Dp;
    while (--n);
    if (*bp) // 第1个元素不为0, 可以直接进入追赶过程
    {
        JUDGE_0_DEL(chase(b.n_elem, &a.at(0), B, &c.at(0), dp))
        free(B);
        return;
    }
    if (!--(n = a.n_elem)) // 系数矩阵只有2阶
    {
        const double &C(c.at(0)), &A(a.at(0));
        double &ele_2(*(dp + 1));
        JUDGE_0_DEL(!C || !A)
        *dp = (Dp[1] - B[1] * (ele_2 = *Dp / C)) / A;
        free(B);
        return;
    }
    // const unsigned A_size(sizeof(double)*n);
    // const double *AP(&a.at(0)),*CP(&c.at(0));
    // JUDGE_0_DEL(!*CP)
    // double *A((double*)malloc(A_size)),*C((double*)malloc(A_size)),last_sol(*Dp/ *CP);
    JUDGE_0_DEL(chase_0(n, &a.at(0), B, &c.at(0), dp))
    // free(A);
    free(B);
    // free(C);
}

例 求解线性方程组

$$\left(\begin{array}{cccc}?4& 1& & \\ 1& ?4& 1& \\ & 1& ?4& 1\\ & & 1& ?4\end{array}\right)\left(\begin{array}{c}{x}_1\\ x_2\\ x_3\\ x_4\end{array}\right)=\left(\begin{array}{c}1\\ 1\\ 1\\ 1\end{array}\right)$$

代入程序求得解向量为$x=(?0.3636,?0.4545,?0.4545,?0.3636{)}^T$.