雅克比和高斯赛德尔迭代求方程组的解

实验原理：

雅克比和高斯的不同之处是，k=1下面那个框框，高斯是a[i][j]*y[j]，雅克比是a[i][j]*x[j]，就是因为高斯用的是更新之后的x，雅克比没有用更新之后的x。

在代码编写过程中，我没有用到动态数组，感觉输入系数矩阵啥的用动态数组更好些。但是，我写了转化对角占优的部分，这部分也是非常重要的。

代码如下：

/**
 运用雅克比和高斯赛德尔公式求解方程组
 a 比较两种方法的收敛速度
 b 验证收敛条件的正确性
 **/
 #include<iostream>
 #include<math.h>
 #define N 3//4*4方程组
 double A[N][N];
 double C[N][N];
 double b[N];
 double d[N];
 using namespace std;
/**
调整系数矩阵为对角占优，并得到对角占优阵C
**/
void betterA( double A[N][N],double C[N][N],double b[N])
{
    double B[N][N];
    int pos[N];
    int i,j,pos_ele,k;
    double sum,max_ele;
    for (i=0;i<N;i++)
      {
           for(j=0;j<N;j++)
        {
            B[i][j]=fabs(A[i][j]);
        }

      }
        for(i=0;i<N;i++)
        {
            max_ele=B[i][0];
            pos_ele=0;
            sum=0;
            for(j=0;j<N;j++)
            //找这一行最大的元素
            {
                if(B[i][j] > max_ele)
               {
                   max_ele=B[i][j];
                   pos_ele=j;//第i行最大元素在第j列，要标记这个j
               }
               sum=sum+B[i][j];//这一行的元素之和
            }
            sum=sum-max_ele;
            if(max_ele<=sum)//这一行最大元素仍然小于所在行的其他元素的绝对值之和
                {
                    cout<<"不能转化为对角占优矩阵！因为第"<<i<<"行的最大元素"<<max_ele<<"小于同行其他元素之和:"<<sum<<endl;
                    //cout<<sum<<endl<<max_ele<<endl<<pos_ele<<endl<<j;
                    return;
                }
            else if(max_ele>sum)
            {
                pos[i]=pos_ele;//记录第i行最大元素所在的行的下标
                cout<<"第"<<i<<"行最大元素在第"<<pos_ele<<"列"<<endl;
            }
        }

   /**
    set the row pos of matrix C with A(i,:)
    **/
    for(i=0;i<N;i++)
    {
        for(j=0;j<N;j++)
        {
            C[pos[i]][j]=A[i][j];
            d[pos[i]]=b[i];
        }
    }
    return;
}
void jakbi(double C[N][N], double b[N], double x[N], double y[N], double eps, int n)
{
    int k=1,i,j;
    while(1)//小于最大迭代次数
    {
        for(i=0;i<N;i++)
        {
            x[i]=y[i];
            double s=0;
            for(j=0;j<N;j++)
            {
                if(j-i)
                    s=s+C[i][j]*x[j];
                y[i]=(b[i]-s)/C[i][i];
            }
        }
        double max=0;
        for(i=0;i<N;i++)
        {
            if(max<fabs(y[i]-x[i]))//无穷范数
                max=fabs(y[i]-x[i]);
        }
        if(max<eps)//
            {
                cout<<"雅克比迭代成功！迭代次数是"<<k<<endl;
                cout<<"迭代结果是";
           for(int m=0;m<N;m++)
                  cout<<"x["<<m<<"]="<<y[m]<<endl;
                break;
            }
        else if(max>=eps && k<n)
        {
            k++;
            continue;//不满足精度要求，继续迭代
        }
        if(k==n)
      {
        cout<<"迭代失败"<<endl;
        break;
       }
    }
}
void Gauss(double C[N][N], double b[N], double x[N], double y[N], double eps, int n)
{
    int k=1,i,j;
    while(1)//小于最大迭代次数
    {
        for(i=0;i<N;i++)
        {
            x[i]=y[i];//y[0]初值
            double s=0;
            for(j=0;j<N;j++)
            {
                if(j-i)
                    s=s+C[i][j]*y[j];
                y[i]=(b[i]-s)/C[i][i];
            }
        }
        double max=0;
        for(i=0;i<N;i++)
        {
            if(max<fabs(y[i]-x[i]))
                max=fabs(y[i]-x[i]);
        }
        if(max<eps)//
            {
                cout<<"高斯-赛德尔迭代成功！迭代次数是"<<k<<endl;
                cout<<"迭代结果是";
                for(int m=0;m<N;m++)
                  cout<<"x["<<m<<"]="<<y[m]<<endl;
                break;
            }
        else if(max>=eps && k<n)
        {
            k++;
            continue;//不满足精度要求，继续迭代
        }
        if(k==n)
      {
        cout<<"迭代失败"<<endl;
        break;
       }
    }
}
int main()
{
    /**
 输入系数矩阵
 **/
    int i,j,n;
    double eps;
    cout<<"请输入系数矩阵"<<endl;
    for(i=0;i<N;i++)
        for(j=0;j<N;j++)
        cin>>A[i][j];
    cout<<"请输入b["<<N<<"]"<<endl;
    for(i=0;i<N;i++)
        cin>>b[i];
    betterA(A,C,b);
    cout<<"转化为对角占优阵是"<<endl;
    for(i=0;i<N;i++)
    {
        for(j=0;j<N;j++)
            cout<<C[i][j]<<" ";
        cout<<d[i]<<endl;//转化后的b
    }
    cout<<"请输入精度和最大迭代次数";
    cin>>eps>>n;
    double x[N];
    double y[N];
    for(i=0;i<N;i++)
        {
            x[i]=0;
            y[i]=0;//赋初值
        }
    jakbi(C,d,x,y,eps,n);
    for(i=0;i<N;i++)
    {
        x[i]=0;
        y[i]=0;
    }
    Gauss(C,d,x,y,eps,n);
    return 0;
}