C++实现托马斯算法求解三对角矩阵输出x1=0错误排查

问题现象 #

使用托马斯算法求解三对角线性方程组时，手动计算得到的正确解为{-1, 2, -3, 4}，但程序输出的x1始终为0，不符合预期。

原始代码 #

#include<iostream>
#include<cmath>
using namespace std;

int main()
{
    float mat[100][100];
    int n;
    cout << "Enter the dimention of the matrix: ";
    cin >> n;
    cout << "Enter the matrix: ";
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < n; j++)
        {
            cin >> mat[i][j];
        }
    }
    float d[100];
    cout << "Enter the values of di's: ";
    for (int i = 0; i < n; i++)
    {
        cin >> d[i];
    }
    for (int i = 0; i < n; i++)
    {
        cout << "The value of d" << i + 1 << " is" << d[i];
    }
    cout << " \n";

    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < n; j++)

            cout << mat[i][j] << " ";
        {
            cout << endl;
        }

    }

    float a[100], b[100], c[100];
    for (int i = 0; i < n + 1; i++)
    {
        b[i] = mat[i][i];
    }
    for (int l = 0; l < n; l++)
    {
        cout << "The value of b" << l + 1 << " is " << b[l] << " " << endl;
    }
    cout << "\n";

    for (int j = 0; j < n; j++)
    {
        a[j] = mat[j][j - 1];
    }

    for (int i = 1; i < n - 1; i++)
    {
        cout << "The value of a" << i + 1 << " is: " << a[i] << " " << endl;
    }
    cout << "\n";
    for (int k = 0; k < n - 1; k++)
    {
        c[k] = mat[k][k + 1];
    }
    for (int i = 0; i < n - 1; i++)
    {
        cout << "The value of c" << i + 1 << " is: " << c[i] << " " << endl;
    }
//to find alpha
    cout << "\n";
    float alpha[100];
    alpha[0] = b[0];
    for (int i = 1; i < n; i++)
    {
        alpha[i] = b[i] - ((a[i] * c[i - 1]) / alpha[i - 1]);
    }
    for (int i = 0; i < n; i++)
    {
        cout << "The value of alpha" << i + 1 << " is: " << alpha[i] << endl;
    }
    cout << "\n";
//to find beta
    float beta[100];
    beta[0] = (d[0] / b[0]);
    for (int i = 1; i < n; i++)
    {
        beta[i] = ((d[i] - a[i] * beta[i - 1]) / alpha[i]);
    }
    for (int i = 0; i < n; i++)
    {
        cout << " The value of beta" << i + 1 << " is: " << beta[i] << "\n";
    }

    //finding the solutions
    float x[100];
    x[n - 1] = beta[n - 1];
    for (int i = n - 2; i > 0; i--)
    {
        x[i] = beta[i] - ((c[i] * x[i + 1]) / alpha[i]);
    }
    cout << "the solutions are ";
    for (int i = 0; i < n - 1; i++)
    {
        cout << x[i] << "\t";
    }
    cout << beta[n - 1] << "\t";
}

运行输出 #

PS C:\Users\jitub\desktop> ./a.exe
Enter the dimention of the matrix: 4                                                                           
Enter the matrix: 3 -1 0 0 -1 3 -1 0 0 -1 3 -1 0 0 -1 3
Enter the values of di's: -5 10 -15 15
The value of d1 is-5The value of d2 is10The value of d3 is-15The value of d4 is15 
3 -1 0 0 
-1 3 -1 0 
0 -1 3 -1 
0 0 -1 3 
The value of b1 is 3 
The value of b2 is 3 
The value of b3 is 3 
The value of b4 is 3 

The value of a2 is: -1
The value of a3 is: -1 

The value of c1 is: -1
The value of c2 is: -1
The value of c3 is: -1

The value of alpha1 is: 3
The value of alpha2 is: 2.66667
The value of alpha3 is: 2.625
The value of alpha4 is: 2.61905

 The value of beta1 is: -1.66667
 The value of beta2 is: 3.125
 The value of beta3 is: -4.52381
 The value of beta4 is: 4
 the solutions are 0     2       -3      4

问题原因与修复方案 #

1. 反向求解循环逻辑错误 #

原始代码反向计算x的循环为：

for (int i = n - 2; i > 0; i--)

循环终止条件为i>0，当i递减到0时直接跳出循环，从未计算x[0]的值，x[0]是未初始化的栈内存值，本次运行刚好为0。
修复：将循环条件改为i >= 0

2. a数组赋值越界 #

原始代码给a数组赋值的循环为：

for (int j = 0; j < n; j++)
{
    a[j] = mat[j][j - 1];
}

j=0时访问mat[0][-1]属于数组越界，会触发未定义行为。三对角矩阵的下对角线元素a[i]从第二行（索引1）开始存在，a[0]无实际意义。
修复：将a数组的赋值循环改为：

for (int j = 1; j < n; j++)
{
    a[j] = mat[j][j - 1];
}

3. 打印逻辑冗余 #

原始代码打印解时单独打印beta[n-1]，实际x[n-1]已经提前赋值为beta[n-1]，直接遍历整个x数组即可，避免遗漏元素。
修复：将打印逻辑改为：

cout << "the solutions are ";
for (int i = 0; i < n; i++)
{
    cout << x[i] << "\t";
}

修复后运行结果 #

输出为-1 2 -3 4，与手动计算结果一致。

内容的提问来源于stack exchange，提问作者Jitu Biswas