在一些实际问题中,例如解常微分方程边值问题,解热传导方程以及船体数学放样中建立三次样条函数等,都会要求解系数矩阵为对角占优的三对角线方程组;——《数值分析(第五版)》

追赶法解三对角方程组例题Python 追赶法求三对角方程_矩阵

 即为Ax=f;

其中,当

追赶法解三对角方程组例题Python 追赶法求三对角方程_线性代数_02

>1时,

追赶法解三对角方程组例题Python 追赶法求三对角方程_矩阵_03

=0时,且满足下面三个条件时,即可对A矩阵进行LU分解,求x矩阵;1.

追赶法解三对角方程组例题Python 追赶法求三对角方程_线性代数_04\left |c _{1} \right |>0" title="\left |b _{1}\right |>\left |c _{1} \right |>0" style="width: 114px; visibility: visible;" data-type="block">

;2.

追赶法解三对角方程组例题Python 追赶法求三对角方程_c++_05\left |a _{n} \right |>0" title="\left |b _{n}\right |>\left |a _{n} \right |>0" style="width: 118px; visibility: visible;" data-type="block">

;3.

追赶法解三对角方程组例题Python 追赶法求三对角方程_矩阵_06

,

追赶法解三对角方程组例题Python 追赶法求三对角方程_i++_07


追赶法解三对角方程组例题Python 追赶法求三对角方程_i++_08

LU分解:

追赶法解三对角方程组例题Python 追赶法求三对角方程_c++_09

其中

追赶法解三对角方程组例题Python 追赶法求三对角方程_c++_10


三对角线方程组的追赶法公式:

1.计算{

追赶法解三对角方程组例题Python 追赶法求三对角方程_矩阵_11

}的递推公式;

追赶法解三对角方程组例题Python 追赶法求三对角方程_矩阵_12

2.解Ly=f;

追赶法解三对角方程组例题Python 追赶法求三对角方程_追赶法解三对角方程组例题Python_13

 3.解Ux=y; 

追赶法解三对角方程组例题Python 追赶法求三对角方程_线性代数_14

 程序主要思路:

赋值模块;

在程序中使用的赋值方式是按照对角线的方向进行赋值,由于只有三对角线,其他数值为0.主对角线为数组b,下对角线为数组a,上对角线为数组c。这种赋值方式是由追赶公式决定的。

求解函数;

程序中一共采用了4个调用函数进行求解,分别为:求解U矩阵函数、求解L矩阵函数、求解Y矩阵函数以及求解方程组的解x函数。由于追赶法不是一个迭代的过程,追赶法可以求解完β和α之后,再进行求解y,最后求解x。所以在主函数中并不需要循环语句求解。

 输出模块;

主要解释L矩阵和U矩阵的输出方式,在输出矩阵之前,我先定义了一个元素全为0的同阶矩阵(这主要是总结了LU矩阵的特点),然后只要将0矩阵相应的位置上的元素重新赋值为求解出来的数组L和数组U。不过别忘了,还需要将L矩阵的下对角线的元素赋值为数组a,U的主对角线赋值为1。

附上实现程序(Visual Studio 2019)

#include <iostream> 
#include <iomanip>//用于cout保留小数

using namespace std;
double b[100]; //定义系数矩阵的主对角线b
double a[100]; //定义系数矩阵下对角线a
double c[100]; //定义系数矩阵上对角线c
double L[100]; //定义拆解后的矩阵L
double U[100]; //定义拆解后的矩阵U
double Y[100]; //定义矩阵Y
double f[100]; //定义矩阵f
double x[100]; //定义求解结果x
double X_x[100][100];//定义输出矩阵的载体
void Func_U(double c[], double a[], double b[], int n);//定义求解U矩阵函数
void Func_L(double b[], double a[], double U[], int n);//定义求解L矩阵函数
void Func_Y(double f[], double a[], double L[], int n);//定义求解Y矩阵函数
void Func_x(double U[], double Y[], int n);//定义求解最终方程组的解x函数 

int main()
{
    int N; //定义矩阵阶数
    cout << "请输入系数矩阵的阶数:" << endl;
    cin >> N;
    cout << "请输入原始矩阵的主对角线b上的元素" << endl;
    for (int i = 1; i <= N; i++)
    {
        cin >> b[i];
    }
    cout << "请输入原始矩阵的下对角线a上的元素" << endl;
    for (int i = 2; i <= N; i++)
    {
        cin >> a[i];
    }
    cout << "请输入原始矩阵的上对角线c上的元素" << endl;
    for (int i = 1; i <= N - 1; i++)
    {
        cin >> c[i];
    }
    cout << "请输入f矩阵元素" << endl;
    for (int i = 1; i <= N; i++)
    {
        cin >> f[i];
    }  
    //以上为赋值模块
    Func_U(c, a, b, N); //调用求解U矩阵函数
    Func_L(b, a, U, N); //调用求解L矩阵函数
    Func_Y(f, a, L, N); //调用求解Y矩阵函数
    Func_x(U, Y, N); //调用求解方程组的解x函数
    for (int i = 1; i <= N; i++)
    {

        for (int j = 1; j <= N; j++)
        {
            X_x[i][j] = 0;
        }
    }//建立载体矩阵,由于LU分解矩阵的特点,先将矩阵的元素赋值为0,之后输出分解矩阵时只要将变化的元素重新赋值就行了
    for (int i = 2; i <= N; i++)
    {
        X_x[i][i - 1] = a[i];
    }//将载体矩阵的下对角线赋值为原始系数矩阵对角线a上的对应元素
    for (int i = 1; i <= N; i++)
    {
        X_x[i][i] = L[i];
    }//将载体矩阵的主对角线赋值为调用求解L矩阵函数求出的对应元素
    cout << "输出L矩阵" << endl;
    cout << setiosflags(ios::fixed) << setprecision(5);//输出结果保留5位小数
    for (int i = 1; i <= N; i++)
    {
        for (int j = 1; j <= N; j++)
        {
            cout << X_x[i][j] << " ";
        }
        cout << endl;
    }//输出L矩阵
    for (int i = 2; i <= N; i++)
    {
        X_x[i][i - 1] = 0;
    }//将载体矩阵的下对角线赋值为0
    for (int i = 1; i <= N; i++)
    {
        X_x[i][i] = 1;
    }//将载体矩阵的主对角线赋值为1
    for (int i = 2; i <= N; i++)
    {
        X_x[i - 1][i] = U[i - 1];
    }//将载体矩阵的上对角线赋值为调用求解U矩阵函数求出的对应元素
    cout << "输出U矩阵" << endl;
    for (int i = 1; i <= N; i++)
    {
        for (int j = 1; j <= N; j++)
        {
            cout << X_x[i][j] << " ";
        }
        cout << endl;
    }//输出U矩阵
    cout << "输出Y矩阵" << endl;
    for (int i = 1; i <= N; i++)
    {
        cout << Y[i] << " ";
    }//输出Y矩阵
    cout << endl;
    cout << "输出x矩阵" << endl;
    for (int i = 1; i <= N; i++)
    {
        cout << x[i] << endl;
    }//输出x矩阵
    return 0;
}

void Func_U(double c[], double a[], double b[], int n)
{
    U[1] = c[1] / b[1];
    for (int i = 2; i <= n - 1; i++)
    {
        U[i] = c[i] / (b[i] - a[i] * U[i - 1]);
    }
}

void Func_L(double b[], double a[], double U[], int n)
{
    L[1] = b[1];
    for (int i = 2; i <= n; i++)
    {
        L[i] = b[i] - a[i] * U[i - 1];
    }
}

void Func_Y(double f[], double a[], double L[], int n)
{
    Y[1] = f[1] / L[1];
    for (int i = 2; i <= n; i++)
    {
        Y[i] = (f[i] - a[i] * Y[i - 1]) / L[i];
    }
}

void Func_x(double U[], double Y[], int n)
{
    x[n] = Y[n];
    for (int i = 1; i <= n - 1; i++)
    {
        x[n - i] = Y[n - i] - U[n - i] * x[n + 1 -i ];
    }
}

例题:

追赶法解三对角方程组例题Python 追赶法求三对角方程_线性代数_15

 程序输出结果:

追赶法解三对角方程组例题Python 追赶法求三对角方程_线性代数_16

 P.S:在输出L和U矩阵的时候,没有考虑到负号的占位,所以导致了这种不对齐的输出结果。

 感谢

在输出矩阵的时候,小白的我不知道如何在里面的循环输出完成后,换行继续输出下一次循环,以此类推直到输出一个完整矩阵。