线性方程组数值解法（3）

 

雅可比迭代法：

#include<stdio.h>

#include<stdlib.h>

#include<math.h>

#define eps 1e-4//题目要求精度

#define MaxSize 100//最大迭代次数

#define n 6//矩阵的维数

double a[n][n] = {

             {4,-1,0,-1,0,0},

             {-1,4,-1,0,-1,0},

             {0,-1,4,-1,0,-1},

             {-1,0,-1,4,-1,0},

             {0,-1,0,-1,4,-1},

             {0,0,-1,0,-1,4} };

//如果是一列的向量可以用一维数组存储

double b[n]={0,5,-2,5,-2,6};//b

double X[n]={0,0,0,0,0,0};//解向量

double f[n]={0,0,0,0,0,0};//Invdiag*b

double Invdiag[n][n]={0};//D的逆

double Mat_Jacobi[n][n] = { 0 };//B=Invdiag*(L+U)

double temp[n]={0,0,0,0,0,0};//保存旧的X

double add1[n]={0,0,0,0,0,0};//保存B*X

double d;

//计算系数矩阵的Jacobi迭代矩阵B

void Jacobi_Mat()

{

      int i,j;

      for(j=0;j<n;j++)

      {

             if(a[j][j] != 0)

             {

                    d=a[j][j];

                    Invdiag[j][j] = 1/a[j][j];//D的逆矩阵

                    for(i=0;i<n;i++)

                    {

                           if(i!=j)

                           {

                                  Mat_Jacobi[i][j]=(-a[i][j])/d;//B矩阵

                           }

                           else

                           {

                                  Mat_Jacobi[i][j] = 0;

                           }

                    }

             }

             else

             {

                    printf("对角元素存在零元,无法求取Jacobi迭代矩阵!\n");

                    return;

             }

      }

}

//由于列向量的特殊性，二范数相当于各分量的绝对值的平方和的开方

double GetNorm()

{

      int i;

      double sum=0;

      double e;

      for (i=0;i<n;i++)

      {

             sum=sum+fabs(X[i]-temp[i])*fabs(X[i]-temp[i]);//计算X1-X0

             e=sqrt(sum);//计算X1-X0的二范数

      }

      return e;

}

//雅可比迭代公式

void Jacobi()

{

      //i表示行，j表示列，k用来计算迭代次数

      int i,j,k;

      double max,sum;

      //判断系数矩阵是否严格对角占优，雅可比迭代法是否收敛

      for(i=0;i<n;i++)

      {

             sum=0;

             max=fabs(a[i][i]);

             for(j=0;j<n;j++)

             {

                    if(i!=j)

                    {

                           sum=sum+fabs(a[i][j]);

                    }

             }

             if(sum>=max)

             {

                    printf("雅可比迭代法发散！\n");

             }

      }

      //计算f=(D逆)*b

      for(i=0;i<n;i++)

      {

             f[i]=Invdiag[i][i]*b[i];

      }

      for(k=0;k<MaxSize;k++)

      {

             //计算add1=B*x

             for(i=0;i<n;i++)

             {

                    sum=0;

                    for(j=0;j<n;j++)

                    {

                           sum=sum+Mat_Jacobi[i][j]*X[j];

                    }

                    add1[i]=sum;

             }

             //保存前一次的X解向量到temp

             for(i=0;i<n;i++)

             {

                    temp[i]=X[i];

             }

             //计算新的X=add1+f

             for(i=0;i<n;i++)

             {

                    X[i]=add1[i]+f[i];

             }

             //打印第k+1次迭代的解向量

             if(GetNorm()<eps)