算法不难,只要理解Jacobi的算法另外两种就很简单了,因为他们都是建立在Jacobi迭代的基础上.
首先讲一下程序,主要是由输入功能,和迭代方式调用功能实现块构成的.
点击这里下载源码.
提取码:mlpb
程序变量
int m; //m*n的系数矩阵
int max_times;// z最大迭代次数
int times = 0;//当前迭代次数
//int mx_mis = 99999999;
string choice;//功能选择
double mis;//误差范围
double mtrx[sz][sz];//系数矩阵
double b[sz];//常数矩阵
double init[sz];//初始向量
double ans[sz];//答案
//double max_mis[sz];//
主程序
int main()
{
input();
func();
return 0;
}
输入功能
这里就是输入计算需要的系数矩阵,常数矩阵等基本条件.
void input()
{
cout<<"请输入迭代方式(J,GS,SOR)"<<endl;
cin>>choice;
cout << "请输入m" << endl;
cin >> m;
cout << "请输入系数矩阵m*n(m)" << endl;
for ( int i = 0; i < m; i++ )
{
for ( int j = 0; j < m; j++ )
{
cin >> mtrx[i][j];
}
}
cout << "请输入常数矩阵b:" << endl;
for ( int i = 0; i < m; i++ )
{
cin >> b[i];
}
cout << "请输入初始向量" << endl;
for ( int i = 0; i < m; i++ )
{
cin >> init[i];
ans[i] = init[i];
}
cout << "请输入最大迭代次数:" << endl;
cin >> max_times;
// cout << endl << "请输入最大误差:"<<endl;
// cin>>mis;
}
输出功能
迭代结束后进行输出
void output()
{
cout << "迭代" << times+1 << "次,答案是" << endl;
for ( int i = 0; i < m; i++ )
{
// cout << init[i] << endl;
printf ( "%.4f\t", init[i] );
}
printf ( "\n" );
}
迭代调用功能
为了main函数里干净点,把功能函数统统放进func里
void func()
{
if(choice=="J")
Jacobi();
else if(choice=="GS")
GS();
else
SOR();
}
Jacobi迭代
首先确立迭代截止条件,就是迭代次数达到所设(或及误差满足条件,但我这里么有实现).然后在while里进行具体迭代.
形如:
[
10
−
1
−
2
−
1
10
−
2
−
1
−
1
5
]
[
x
1
x
2
x
3
]
=
[
7.2
8.3
4.2
]
\left[ \begin{matrix} 10 & -1 & -2 \\ -1 & 10 & -2 \\ -1 & -1 & 5 \end{matrix} \right] \left[ \begin{matrix} x_1\\ x_2 \\ x_3 \end{matrix} \right] = \left[ \begin{matrix} 7.2\\ 8.3\\ 4.2 \end{matrix} \right]
⎣⎡10−1−1−110−1−2−25⎦⎤⎣⎡x1x2x3⎦⎤=⎣⎡7.28.34.2⎦⎤
这个矩阵.进行Jacobi迭代时,按照迭代格式
x
i
k
+
1
=
1
a
i
i
(
d
i
−
∑
j
=
1
i
−
1
a
i
j
k
−
∑
j
=
i
+
1
n
a
i
j
k
)
x_i^{k+1}= \frac {1} {a_{ii}} \left(d_i- \sum_{j=1}^{i-1}a_{ij}^{k}-\sum_{j=i+1}^{n}a_{ij}^k \right)
xik+1=aii1(di−j=1∑i−1aijk−j=i+1∑naijk)
需要先进行类似矩阵乘法的操作,将其求和并用
d
i
d_i
di减去,然后除
a
i
i
a_{ii}
aii,得到的值就是
x
i
k
+
1
x_i^{k+1}
xik+1.
void Jacobi()
{
do
{
mx_mis = 0;
//根据初始向量进行向量乘法.求出每一行对应的解
//把结果放在ans中
for ( int i = 0; i < m; i++ )
{
double sum = 0;
for ( int j = 0; j < m; j++ )
{
if ( j != i )
{
sum += init[j] * mtrx[i][j];
}
}
ans[i] = ( b[i] - sum ) / mtrx[i][i];
// if ( i > 0 )
if ( fabs ( init[i] - ans[i] ) > mx_mis )
mx_mis = fabs ( init[i] - ans[i] );
}
//迭代完成一次, 就更新一次init数组,以便下一次迭代
for ( int i = 0; i < m; i++ )
{
init[i] = ans[i];
}
// mx_mis=max_mis[0];
// for(int i=0;i<m;i++)
// {
// if(max_mis[i]>mx_mis)
// mx_mis= max_mis[i];
// }
output();
}
while ( ( times++ < max_times ) && ( mis_range < mx_mis ) );
}
GS迭代
迭代格式是
x
i
k
+
1
=
1
a
i
i
(
d
i
−
∑
j
=
1
i
−
1
a
i
j
k
−
∑
j
=
i
+
1
n
a
i
j
k
+
1
)
x_i^{k+1}= \frac {1} {a_{ii}} \left(d_i- \sum_{j=1}^{i-1}a_{ij}^{k}-\sum_{j=i+1}^{n}a_{ij}^{k+1} \right)
xik+1=aii1(di−j=1∑i−1aijk−j=i+1∑naijk+1)
对比Jacobi迭代,在代码上的实现区别就是不需要ans数组暂存,直接将新的
x
i
k
+
1
x_i^{k+1}
xik+1赋给init[i].自然也不需要init和ans的同步.
void GS()
{
double temp;
do
{
mx_mis = 0;
//根据初始向量进行向量乘法.求出每一行对应的解
for ( int i = 0; i < m; i++ )
{
double sum = 0;
for ( int j = 0; j < m; j++ )
{
if ( j != i )
{
sum += init[j] * mtrx[i][j];
}
}
temp = init[i];
init[i] = ( b[i] - sum ) / mtrx[i][i];
if ( fabs ( init[i] - temp ) > mx_mis )
mx_mis = fabs ( init[i] - temp );
}
output();
}
while ( ( times++ < max_times ) && ( mx_mis > mis_range ) );
}
SOR迭代
void SOR()
{
double w;
cout << "请输入松弛因子w:" << endl;
cin >> w;
double temp;
do
{
mx_mis = 0;
//根据初始向量进行向量乘法.求出每一行对应的解
for ( int i = 0; i < m; i++ )
{
double sum = 0;
for ( int j = 0; j < m; j++ )
{
if ( j != i )
{
sum += init[j] * mtrx[i][j];
}
}
temp = init[i];
init[i] = init[i] + w * ( ( b[i] - sum ) / mtrx[i][i] - init[i] );
if ( fabs ( init[i] - temp ) > mx_mis )
mx_mis = fabs ( init[i] - temp );
}
output();
}while ( ( times++ < max_times ) && ( mx_mis > mis_range ) );
}
上面是注意分析,下面是源代码.
如有错误请多指正.
#include<iostream>
#include<math.h>
#define sz 10
using namespace std;
/*****************************
@Jacobi迭代解线性方程组
@author:jawi
@date:2020/03/03
@tip:本程序默认输入矩阵均为理想矩阵
*****************************/
void input();
void Jacobi();
void GS();
void SOR();
void func();
void output();
int m; //m*n的系数矩阵
int max_times;// z最大迭代次数
int times = 0;//当前迭代次数
double mx_mis = 0;//当前最大误差
string choice;//功能选择
double mis;//误差范围
double mtrx[sz][sz];//系数矩阵
double b[sz];//常数矩阵
double init[sz];//初始向量
double ans[sz];//答案
double mis_range;//最大误差范围
int main()
{
input();
func();
return 0;
}
void func()
{
if ( choice == "J" )
Jacobi();
else if ( choice == "GS" )
GS();
else
SOR();
}
void input()
{
cout << "请输入迭代方式(J,GS,SOR)" << endl;
cin >> choice;
cout << "请输入m" << endl;
cin >> m;
cout << "请输入系数矩阵m*n(m)" << endl;
for ( int i = 0; i < m; i++ )
{
for ( int j = 0; j < m; j++ )
{
cin >> mtrx[i][j];
}
}
cout << "请输入常数矩阵b:" << endl;
for ( int i = 0; i < m; i++ )
{
cin >> b[i];
}
cout << "请输入初始向量" << endl;
for ( int i = 0; i < m; i++ )
{
cin >> init[i];
ans[i] = init[i];
}
cout << "请输入最大迭代次数:" << endl;
cin >> max_times;
cout << "请输入最大误差范围:" << endl;
cin >> mis_range;
// cout << endl << "请输入最大误差:"<<endl;
// cin>>mis;
}
void Jacobi()
{
do
{
mx_mis = 0;
//根据初始向量进行向量乘法.求出每一行对应的解
//把结果放在ans中
for ( int i = 0; i < m; i++ )
{
double sum = 0;
for ( int j = 0; j < m; j++ )
{
if ( j != i )
{
sum += init[j] * mtrx[i][j];
}
}
ans[i] = ( b[i] - sum ) / mtrx[i][i];
// if ( i > 0 )
if ( fabs ( init[i] - ans[i] ) > mx_mis )
mx_mis = fabs ( init[i] - ans[i] );
}
//迭代完成一次, 就更新一次init数组,以便下一次迭代
for ( int i = 0; i < m; i++ )
{
init[i] = ans[i];
}
// mx_mis=max_mis[0];
// for(int i=0;i<m;i++)
// {
// if(max_mis[i]>mx_mis)
// mx_mis= max_mis[i];
// }
output();
}
while ( ( times++ < max_times ) && ( mis_range < mx_mis ) );
}
void GS()
{
double temp;
do
{
mx_mis = 0;
//根据初始向量进行向量乘法.求出每一行对应的解
for ( int i = 0; i < m; i++ )
{
double sum = 0;
for ( int j = 0; j < m; j++ )
{
if ( j != i )
{
sum += init[j] * mtrx[i][j];
}
}
temp = init[i];
init[i] = ( b[i] - sum ) / mtrx[i][i];
if ( fabs ( init[i] - temp ) > mx_mis )
mx_mis = fabs ( init[i] - temp );
}
output();
}
while ( ( times++ < max_times ) && ( mx_mis > mis_range ) );
}
void SOR()
{
double w;
cout << "请输入松弛因子w:" << endl;
cin >> w;
double temp;
do
{
mx_mis = 0;
//根据初始向量进行向量乘法.求出每一行对应的解
for ( int i = 0; i < m; i++ )
{
double sum = 0;
for ( int j = 0; j < m; j++ )
{
if ( j != i )
{
sum += init[j] * mtrx[i][j];
}
}
temp = init[i];
init[i] = init[i] + w * ( ( b[i] - sum ) / mtrx[i][i] - init[i] );
if ( fabs ( init[i] - temp ) > mx_mis )
mx_mis = fabs ( init[i] - temp );
}
output();
}while ( ( times++ < max_times ) && ( mx_mis > mis_range ) );
}
void output()
{
cout << "迭代" << times+1 << "次,答案是" << endl;
for ( int i = 0; i < m; i++ )
{
// cout << init[i] << endl;
printf ( "%.4f\t", init[i] );
}
printf ( "\n" );
}
原文链接: https://www.cnblogs.com/klaus08/p/15105036.html
欢迎关注
微信关注下方公众号,第一时间获取干货硬货;公众号内回复【pdf】免费获取数百本计算机经典书籍;
也有高质量的技术群,里面有嵌入式、搜广推等BAT大佬
原创文章受到原创版权保护。转载请注明出处:https://www.ccppcoding.com/archives/333400
非原创文章文中已经注明原地址,如有侵权,联系删除
关注公众号【高性能架构探索】,第一时间获取最新文章
转载文章受原作者版权保护。转载请注明原作者出处!