一、实验目的
1、熟悉求解线性方程组的有关理论和方法;
2、会编制LU 分解法、雅可比及高斯—塞德尔迭代法德程序;
3、通过实际计算,进一步了解各种方法的优缺点,选择合适的数值方法。
二、算法描述
3.1 矩阵直接三角分解法
算法
将方程组Ax=b 中的A分解为A=LU,其中L为单位下三角矩阵,U为上三角矩阵,则方程组Ax=b化为解2个方程组Ly=b,Ux=y,具体算法如下:
3.2 迭代法
3.2.1 雅可比迭代法
算法:设方程组Ax=b系数矩阵的对角线元素,M为迭代次数容许的最大值,ε为容许误差。
3.2.2 高斯-塞尔德迭代法
算法:设方程组Ax=b的系数矩阵的对角线元素,,M为迭代次数容许的最大值,ε为容许误差
三、源程序
1、矩阵直接三角分解法
#include <iostream>
#include <math.h>
#include <iomanip>
using namespace std;
/*
1 2 -12 8
5 4 7 -2
-3 7 9 5
6 -12 -8 3
27 4 11 49
*/
int main()
{
cout << "矩阵直接三角分解法\n\n请输入矩阵阶数:" << endl;
int n;
cin >> n;
//数组a[n+1][n+1];
int** a = new int* [n+1];
for (int i = 0; i <= n; i++)
{
a[i] = new int[n+1];
}
//上三角矩阵
double** u = new double* [n + 1];
for (int i = 0; i <= n; i++)
{
u[i] = new double[n + 1];
}
for (int i = 1; i <= n; i++)
{
for (int j = 1; j <= n; j++)
{
u[i][j]=0;
}
}
//单位下三角矩阵
double** l = new double* [n + 1];
for (int i = 0; i <= n; i++)
{
l[i] = new double[n + 1];
}
for (int i = 1; i <= n; i++)
{
for (int j = 1; j <= n; j++)
{
if(i==j)l[i][j] = 1;
else l[i][j] = 0;
}
}
double* b = new double[n+1];
double* y = new double[n + 1];
double* x = new double[n + 1];
cout << "请输入矩阵A的值:" << endl;
for (int i = 1; i<=n; i++)
{
for (int j =1 ; j<=n; j++)
{
cin >> a[i][j];
}
}
cout << "请输入矩阵b的值:" << endl;
for (int i = 1; i<=n; i++)
{
cin >> b[i];
}
for (int j = 1; j <= n; j++)
{
u[1][j] = a[1][j];
}
for (int i = 2; i <= n; i++)
{
l[i][1] = a[i][1]*1.0 / a[1][1];
}
for (int k = 1; k <= n; k++)
{
for (int j = k; j <= n; j++)
{
u[k][j] = a[k][j];
for (int q = 1; q <= k - 1; q++)
{
u[k][j] -= l[k][q] * u[q][j];
}
}
for (int i = k + 1; i <= n; i++)
{
l[i][k] = a[i][k];
for (int q = 1; q <= k - 1; q++)
{
l[i][k] -= l[i][q] * u[q][k];
}
l[i][k] = l[i][k]*1.0/u[k][k];
}
}
y[1] = b[1];
for (int k = 2; k <= n; k++)
{
y[k] = b[k];
for (int q = 1; q <= k - 1; q++)
{
y[k] -= l[k][q] * y[q];
}
}
x[n] = y[n]*1.0 / u[n][n];
for (int k = n - 1; k >= 1; k--)
{
x[k] = y[k];
for (int q = k + 1; q <= n; q++)
{
x[k] -= u[k][q] * x[q];
}
x[k] = x[k] *1.0/u[k][k];
}
cout << "结果为:" << endl;
for (int i = 1; i <= n; i++)
{
cout << "X[" << i-1 << "]=" << setiosflags(ios::fixed) << setprecision(6) << x[i]<<endl;
}
return 0;
}
2、雅可比迭代法
#include <iostream>
#include <math.h>
#include <iomanip>
using namespace std;
/*
5 2 1
2 8 -3
1 -3 -6
8 21 1
*/
int main()
{
cout << "雅可比迭代法\n\n请输入矩阵阶数:" << endl;
int n,count,e;
cin >> n;
//数组a[n+1][n+1];
int** a = new int* [n + 1];
for (int i = 0; i <= n; i++)
{
a[i] = new int[n + 1];
}
double* b = new double[n + 1];
cout << "请输入矩阵A的值:" << endl;
for (int i = 1; i <= n; i++)
{
for (int j = 1; j <= n; j++)
{
cin >> a[i][j];
}
}
for (int i = 1; i <= n; i++)
{
for (int j = 1; j <= n; j++)
{
if (i == j && a[i][j] == 0)
{
cout << "对角线元素不能为0,请重新输入!" << endl;
return 0;
}
}
}
cout << "请输入矩阵b的值:" << endl;
for (int i = 1; i <= n; i++)
{
cin >> b[i];
}
cout << "请输入迭代次数允许的最大值:" << endl;
cin >> count;
cout << "请输入容许误差:" << endl;
cin >> e;
double** x = new double* [n + 1];
for (int i = 0; i <= n; i++)
{
x[i] = new double[count+1];
}
double* x1 = new double[n + 1];//最终结果
//设置初始向量
for (int i = 1; i <= n; i++)
{
x[i][0]=0;
}
int k,num=0;
for (k = 0;;)
{
int sum = 0;
for (int i = 1; i <= n; i++)
{
x[i][k + 1] = b[i];
for (int j = 1; j <= n; j++)
{
if (j != i)
{
x[i][k + 1] -= a[i][j] * x[j][k];
}
}
x[i][k + 1] /= a[i][i];
x1[i-1] = x[i][k + 1];
sum += abs(x[i][k + 1] - x[i][k]);
}
if (sum >= e)
{
if (k >= count)
{
//cout << "超出迭代次数允许的最大值,程序结束";
break;
}
else
{
k = k + 1;
}
}
}
for (int i = 0; i <n; i++)
{
cout << "X[" << i << "]=" << setiosflags(ios::fixed) << setprecision(6) << x1[i] << endl;
}
return 0;
}
3、高斯-赛德尔迭代法
#include <iostream>
#include <math.h>
#include <iomanip>
using namespace std;
/*
8 -3 2
4 11 -1
6 3 12
20 33 36
*/
int main()
{
cout << "高斯-赛德尔迭代法\n\n请输入矩阵阶数:" << endl;
int n, count, e;
cin >> n;
//数组a[n+1][n+1];
int** a = new int* [n + 1];
for (int i = 0; i <= n; i++)
{
a[i] = new int[n + 1];
}
double* b = new double[n + 1];
cout << "请输入矩阵A的值:" << endl;
for (int i = 1; i <= n; i++)
{
for (int j = 1; j <= n; j++)
{
cin >> a[i][j];
}
}
for (int i = 1; i <= n; i++)
{
for (int j = 1; j <= n; j++)
{
if (i == j && a[i][j] == 0)
{
cout << "对角线元素不能为0,请重新输入!" << endl;
return 0;
}
}
}
cout << "请输入矩阵b的值:" << endl;
for (int i = 1; i <= n; i++)
{
cin >> b[i];
}
cout << "请输入迭代次数允许的最大值:" << endl;
cin >> count;
cout << "请输入容许误差:" << endl;
cin >> e;
double** x = new double* [n + 1];
for (int i = 0; i <= n; i++)
{
x[i] = new double[count + 1];
}
double* x1 = new double[n + 1];//最终结果
//设置初始向量
for (int i = 1; i <= n; i++)
{
x[i][0] = 0;
}
int k;
for (k = 0;;)
{
int sum = 0;
for (int i = 1; i <= n; i++)
{
x[i][k + 1] = b[i];
for (int j = 1; j <= i-1; j++)
{
x[i][k + 1] -= a[i][j] * x[j][k+1];
}
for (int j = i+1; j <= n; j++)
{
x[i][k + 1] -= a[i][j] * x[j][k];
}
x[i][k + 1] /= a[i][i];
x1[i - 1] = x[i][k + 1];
sum += abs(x[i][k + 1] - x[i][k]);
}
if (sum >= e)
{
if (k >= count)
{
//cout << "超出迭代次数允许的最大值,程序结束";
break;
}
else
{
k = k + 1;
}
}
}
for (int i = 0; i < n; i++)
{
cout << "X[" << i << "]=" << setiosflags(ios::fixed) << setprecision(6) << x1[i] << endl;
}
return 0;
}