关于雅克比迭代, 高斯塞德尔迭代以及SOR迭代的C++实现

算法不难,只要理解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]

10111101225x1x2x3=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(dij=1i1aijkj=i+1naijk)
需要先进行类似矩阵乘法的操作,将其求和并用

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(dij=1i1aijkj=i+1naijk+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大佬

    关于雅克比迭代, 高斯塞德尔迭代以及SOR迭代的C++实现

原创文章受到原创版权保护。转载请注明出处:https://www.ccppcoding.com/archives/333400

非原创文章文中已经注明原地址,如有侵权,联系删除

关注公众号【高性能架构探索】,第一时间获取最新文章

转载文章受原作者版权保护。转载请注明原作者出处!

(0)
上一篇 2023年3月1日 下午9:09
下一篇 2023年3月1日 下午9:09

相关推荐