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

发布时间:2024年01月22日

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

( b 1 c 1 0 0 … 0 a 2 b 2 c 2 0 … 0 0 a 3 b 3 c 3 … 0 ? ? ? ? ? ? 0 0 0 0 a n b n ) ( x 1 x 2 x 3 ? x n ) = ( d 1 d 2 d 3 ? d n ) \begin{pmatrix}b_1 & c_1 & 0 & 0 & \ldots & 0 \\a_2 & b_2 & c_2 & 0 & \ldots & 0 \\0 & a_3 & b_3 & c_3 & \ldots & 0 \\\vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\0 & 0 & 0 & 0 & a_n & b_n\end{pmatrix}\begin{pmatrix}x_1 \\x_2 \\x_3 \\\vdots \\x_n\end{pmatrix}=\begin{pmatrix}d_1 \\d_2 \\d_3 \\\vdots \\d_n\end{pmatrix} ?b1?a2?0?0?c1?b2?a3??0?0c2?b3??0?00c3??0??an??000?bn?? ? ?x1?x2?x3??xn?? ?= ?d1?d2?d3??dn?? ?

其中, a i , b i , c i a_i, b_i, c_i ai?,bi?,ci? 是已知系数, x i x_i xi? 是未知变量, d i d_i di? 是常数项.

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

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

  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);
}

求解线性方程组

( ? 4 1 1 ? 4 1 1 ? 4 1 1 ? 4 ) ( x 1 x 2 x 3 x 4 ) = ( 1 1 1 1 ) \begin{pmatrix} -4&1&&\\1&-4&1&\\&1&-4&1\\&&1&-4 \end{pmatrix}\begin{pmatrix} x_1\\x_2\\x_3\\x_4 \end{pmatrix}=\begin{pmatrix} 1\\1\\1\\1 \end{pmatrix} ??41?1?41?1?41?1?4? ? ?x1?x2?x3?x4?? ?= ?1111? ?

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

文章来源:https://blog.csdn.net/m0_61219098/article/details/135541579
本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。