在一些实际问题中,例如解常微分方程边值问题,解热传导方程以及船体数学放样中建立三次样条函数等,都会要求解系数矩阵为对角占优的三对角线方程组;——《数值分析(第五版)》
即为Ax=f;
其中,当
>1时,
=0时,且满足下面三个条件时,即可对A矩阵进行LU分解,求x矩阵;1.
\left |c _{1} \right |>0" title="\left |b _{1}\right |>\left |c _{1} \right |>0" style="width: 114px; visibility: visible;" data-type="block">
;2.
\left |a _{n} \right |>0" title="\left |b _{n}\right |>\left |a _{n} \right |>0" style="width: 118px; visibility: visible;" data-type="block">
;3.
,
,
LU分解:
其中
;
三对角线方程组的追赶法公式:
1.计算{
}的递推公式;
2.解Ly=f;
3.解Ux=y;
程序主要思路:
赋值模块;
在程序中使用的赋值方式是按照对角线的方向进行赋值,由于只有三对角线,其他数值为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 ];
}
}
例题:
程序输出结果:
P.S:在输出L和U矩阵的时候,没有考虑到负号的占位,所以导致了这种不对齐的输出结果。
感谢
在输出矩阵的时候,小白的我不知道如何在里面的循环输出完成后,换行继续输出下一次循环,以此类推直到输出一个完整矩阵。