解线性方程组的Jacobi迭代法和Gauss-Seidel迭代法

最新推荐文章于 2023-11-09 22:05:40 发布

ceamky

最新推荐文章于 2023-11-09 22:05:40 发布

阅读量2.8k

点赞数

分类专栏： C/C++ 文章标签： ini math.h include

本文链接：https://blog.csdn.net/ceamky/article/details/1644708

版权

C/C++ 专栏收录该内容

18 篇文章 0 订阅

订阅专栏

#include <math.h>
#include <stdio.h>
#define eps 0.0001
#define max 100
#define N 10

double norm_inf(double x[],int n)
{
double norm;
int i;
norm=fabs(x[0]);
for(i=1;i<n;i++)
{
if(fabs(x[i])>norm)
norm=fabs(x[i]);
}
return norm;
}

void jacobi(double a[N][N],double g[N],int n)
{
double b[N][N]={0},x0[N],x1[N],x1_x0[N],norm,temp;
int i,j,k;
for(i=0;i<n;i++)
{
  g[i]=g[i]/a[i][i];
  for(j=0;j<n;j++)
  {
   if(i==j)
    continue;
   b[i][j]=-a[i][j]/a[i][i];
  }
}
for(i=0;i<n;i++)
{
  x0[i]=0;
  x1[i]=1;
  x1_x0[i]=x1[i]-x0[i];
}
k=0;
norm=norm_inf(x1_x0,n);
while((norm>=eps)&&(k<max))
{
  for(i=0;i<n;i++)
   x0[i]=x1[i];
  for(i=0;i<n;i++)
  {
   temp=0;
   for(j=0;j<n;j++)
    temp=temp+b[i][j]*x0[j];
   x1[i]=temp+g[i];
   x1_x0[i]=x1[i]-x0[i];
  }
  norm=norm_inf(x1_x0,n);
  k++;
}
for(i=0;i<n;i++)
  printf("x[%d]=%lf/n",i,x1[i]);
printf("%d times iteration./n",k);
}

void seidel(double a[N][N],double g[N],int n)
{
double b[N][N]={0},x0[N],x1[N],x1_x0[N],norm,temp;
int i,j,k;
for(i=0;i<n;i++)
{
  g[i]=g[i]/a[i][i];
  for(j=0;j<n;j++)
  {
   if(i==j)
    continue;
   b[i][j]=-a[i][j]/a[i][i];
  }
}
for(i=0;i<n;i++)
{
  x0[i]=0;
  x1[i]=1;
  x1_x0[i]=x1[i]-x0[i];
}
k=0;
norm=norm_inf(x1_x0,n);
while((norm>=eps)&&(k<max))
{
  for(i=0;i<n;i++)
   x0[i]=x1[i];
  for(i=0;i<n;i++)
  {
   temp=0;
   for(j=0;j<=i-1;j++)
    temp=temp+b[i][j]*x1[j];
   for(j=i+1;j<n;j++)
    temp=temp+b[i][j]*x0[j];
   x1[i]=temp+g[i];
   x1_x0[i]=x1[i]-x0[i];
  }
  norm=norm_inf(x1_x0,n);
  k++;
}
for(i=0;i<n;i++)
  printf("x[%d]=%lf/n",i,x1[i]);
printf("%d times iteration./n",k);
}

int main()
{
double a[N][N]={{5,2,3,2},{2,4,1,-2},{1,-3,4,3},{3,2,2,8}},g[N]={-1,5,4,-6};
double b[N][N]={{5,2,3,2},{2,4,1,-2},{1,-3,4,3},{3,2,2,8}},f[N]={-1,5,4,-6};
jacobi(a,g,4);
printf("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~/n");
seidel(b,f,4);
return 0;
}