迭代法求解方程组，结果出不来
# include <stdio.h># include <stdlib.h>
# include <math.h>
int main()
{
int n=3,k;
double **a;
double f[n+1][n+1],x[n+1],y[n+1],u[n+1][n+1],g[n+1][n+1],R[n+1][n+1],b[(n+1)*(n+1)],ex[(n+1)*(n+1)],x0[(n+1)*(n+1)];
double PI=3.14159265359;
double hx,hy,max_error;
int maxn;
hx=PI/n;//x方向的剖分
hy=1.0/n;//y方向的剖分
int n7,n8;//n7-第一维长度,n8-第二维长度
int i,j;
n7=(n+1)*(n+1);
n8=(n+1)*(n+1);
for(i=0;i<=(n+1)*(n+1)-1;i++)//迭代初始值
{
x0[i]=0.0;
}
for(i=0;i<=n;i++)//X,Y 方向上的内节点
{
x[i]=i*hx;
y[i]=i*hy;
}
for(j=0;j<=n;j++)
{
for(i=0;i<=n;i++)
{
f[i][j]=cos(3.0*i*hx)*sin(PI*j*hy);//右边向量
R[i][j]=1.0/(9.0+PI*PI)*cos(3.0*i*hx)*sin(PI*j*hy);//真解
}
}
/*-----------------*/
//数组a[][]开辟空间
a=(double **)malloc(sizeof(double *)*n7);
for(i=0;i<=n7;++i)
{
a[i]=(double *)malloc(sizeof(double)*n8);
}
/*-----------------*/
//对a[][]初始化
for (j=0;j<=(n+1)*(n+1)-1;j++)
for (i=0;i<=(n+1)*(n+1)-1;i++)
a[i][j]=0.0;
/*-----------------*/
//对a[][]操作
for(i=0;i<=n;i++)
{
a[i][i]=1.0;
}
for(k=0;k<=n-1;k++)
{
a[n+1+k*(n+1)][n+1+k*(n+1)]=1;//主对角元中为1的元素
a[2*n+1+k*(n+1)][2*n+1+k*(n+1)]=1;//主对角元中为1的元素
}
for(k=1;k<n;k++)
{
a[k*n+k][k*n+k+1]=-1.0;
a[(k+1)*n+k][(k+1)*n+k-1]=-1.0;
}
for(k=1;k<n;k++)
{
for(i=k*n+k+1;i<=(k+1)*n+k-1;i++)
{
a[i][i-(n+1)]=-1.0/(hy*hy);//下对角上矩阵的元素
a[i][i+n+1]=-1.0/(hy*hy);//上对角上的矩阵的元素
a[i][i-1]=-1.0/(hx*hx);//主对角矩阵中的次对角元的元素
a[i][i+1]=-1.0/(hx*hx);//主对角矩阵中的次对角元的元素
a[i][i]=2.0*(1.0/(hx*hx)+1.0/(hy*hy));//主对角矩阵的主对角线上的元素
}
}
for(i=n*(n+1)-1;i<=(n+1)*(n+1)-1;i++)
{
a[i][i]=1;
}
/*-----------------*/
//输出a
printf("这个特殊矩阵是:\n");
for (i=0;i<=(n+1)*(n+1)-1;i++)
{
for (j=0;j<=(n+1)*(n+1)-1;j++)
printf("%.2f ",a[i][j]);
printf("\n");
}
/*------------*/
//输出右边向量
for(i=0;i<=n;i++)//g[][]初始化
{
for(j=0;j<=n;j++)
g[i][j]=0.0;
}
for(j=1;j<n;j++)
{
for(i=1;i<n;i++)
{
g[i][j]=f[i][j];
}
}
/*—————————*///打印右边向量g[][]
for(i=0;i<=n;i++)
{
for(j=0;j<=n;j++)
{
printf("%.2f",g[i][j]);
}
}
system("pause");
//将矩阵里g的值存放在向量b里;
k = 0;
for(j=0;j<=n;j++)
{
for(i=0;i<=n;i++)
{
b[k+i]=g[i][j];
}
k=k+n+1;
}
/*//将真解存放在向量ex中
k = 0;
for(j=0;j<=n;j++)
{
for(i=0;i<=n;i++)
{
ex[k+i]=R[i][j];
}
k=k+n+1;
}*/
void Jacobi(double **a, double g[][n+1], double x1[(n+1)*(n+1)],double x0[(n+1)*(n+1)], int maxn)//雅克比迭代法
{
int k,i,j;
float sum1,err=0.0;
for (k=0;k<maxn;k++)//k为迭代次数
{
max_error=0.0;
for (i=0;i<=(n+1)*(n+1)-1;i++)
{
for (j=0;j<=(n+1)*(n+1)-1;j++)
{
if(j!=i)
sum1+=a[i][j]*x0[j];
}
x1[i]=(b[i]-sum1)/a[i][i];
sum1=0.0;
}
for(i=0;i<=(n+1)*(n+1)-1;i++)
{
err=fabs(x1[i]-x0[i]);//绝对误差
if(err>max_error)
max_error=err;
x0[i]=x1[i];
}
if (err<pow(10,-6))
break;
}
printf("经过%d次雅克比迭代得到以下解向量\n", k+1);
//对矩阵u赋值0;
for(j=0;j<=n;j++)
for(i=0;i<=n;i++)
u[i][j]=0.0;
//将解向量x存放在矩阵u里;
k = 0;
for(i=0;i<=n;i++)
{
for(j=0;j<=n;j++)
{
u[i][j]=x0[k+j];
}
k=k+n+1;
}
k = n/4;//确定x方向每隔k个打印一下
for(i=k;i<n;i=i+k)
{
printf("(%.2f,0.25),R=%f,u=%f,e=%f\n",x[i],R[i][n/4],u[i][n/4],fabs(u[i][n/4]-R[i][n/4]));
printf("(%.2f,0.5),R=%f,u=%f,e=%f\n",x[i],R[i][n/2],u[i][n/2],fabs(u[i][n/2]-R[i][n/2]));
printf("(%.2f,0.75),R=%f,u=%f,e=%f\n",x[i],R[i][3*n/4],u[i][3*n/4],fabs(u[i][3*n/4]-R[i][3*n/4]));
}
}
void GaussSeidel(double **a, double g[][n+1],double x0[(n+1)*(n+1)],double x1[(n+1)*(n+1)], int maxn)//高斯-赛德尔迭代
{
int k,i,j;
float sum1=0.0,sum2=0.0,err=0, Temporary=0;
for (k=0;k<maxn;k++)//k为迭代次数
{
max_error=0.0;
for (i=0;i<=(n+1)*(n+1)-1;i++)
{
for (j=0;j<=(n+1)*(n+1)-1;j++)
{
for(j=0;j<i;j++)
{
sum1+=a[i][j]*x1[j];
}
for(j=i;j<(n+1)*(n+1);j++)
{
sum2=a[i][j]*x0[j];
}
}
x1[i]=(b[i]-sum1-sum2)/a[i][i];
}
for(i=0;i<=(n+1)*(n+1)-1;i++)
{
err=fabs(x1[i]-x0[i]);//绝对误差
if(err>max_error)
max_error=err;
x0[i]=x1[i];
}
if (err<pow(10,-6))
break;
}
printf("经过%d次高斯-赛德尔迭代得到以下解向量\n", k);
//对矩阵u赋值0;
for(i=0;i<=n;i++)
for(j=0;j<=n;j++)
u[i][j]=0.0;
//将解向量x存放在矩阵u里;
k = 0;
for(i=0;i<=n;i++)
{
for(j=0;j<=n;j++)
{
u[i][j]=x0[k+j];
}
k=k+n+1;
}
k=n/4;//确定x方向每隔k个打印一下
for(i=k;i<n;i=i+k)
{
printf("(%.2f,0.25),R=%f,u=%f,e=%f\n",x[i],R[i][n/4],u[i][n/4],fabs(u[i][n/4]-R[i][n/4]));
printf("(%.2f,0.5),R=%f,u=%f,e=%f\n",x[i],R[i][n/2],u[i][n/2],fabs(u[i][n/2]-R[i][n/2]));
printf("(%.2f,0.75),R=%f,u=%f,e=%f\n",x[i],R[i][3*n/4],u[i][3*n/4],fabs(u[i][3*n/4]-R[i][3*n/4]));
}
printf("请输入迭代最大次数\n");
scanf("%d",&maxn);
printf("请选择算法：\n1.雅克比迭代\n2.高斯赛德尔迭代\n");
scanf("%d",&k);
if(k==1)
Jacobi(a,g,x0,x1,maxn);//调用雅克比迭代法
else
GaussSeidel(a,g,x0,x1,maxn);//调用高斯赛德尔迭代法
}
}
