这个输出结果怎么不对啊？有大佬帮忙看一下吗感谢！！！

三支笔723 2023-03-23 10:39:54

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
int main()
{
 int N,M,i,j,k;
 double x,y,a,c;
 double dltx=1.5;//边界长=1.5
 double dlty=1.0;//边界宽=1
 printf("取横向上的网格数 N 为： ");
 scanf("%d",&N);
 printf("取纵向上的网格数 M 为： ");
 scanf("%d",&M);
 x=dltx/N;
 y=dlty/M;
 printf("所以横向上的步长为：%.2lf\n",x);
 printf("所以纵向上的步长为：%.2lf\n",y);
 printf("\n");
 N=N+1;
 M=M+1;
 double T[N][M],T0[N][M],r[N][M],b[N][M],Sp[N][M];
 double Tf,Lr,Imda1,Imda2,Imda3,h;
 double Eps,DT,DTmax;
 double aw,ae,an,as,ap;
 double q0=0,q1=10;
 Tf=300,h=10;
 Eps=1.e-5;
 Lr=x/y;
 for(k=0;k<4;k++)
 {
     for(i=0;i<N;i++)
     {
         for(j=0;j<M;j++)
         {
             T[i][j]=Tf;
             T0[i][j]=T[i][j];
         }
      }
 loop:
 for(i=0;i<N;i++)
 {
     for(j=0;j<M;j++)
     {
         if(k==0)
         {
            b[i][j]=2*T0[i][j]*T0[i][j]*x*y;
            r[i][j]=0.54+0.00058*(T[i][j]-273);
            Sp[i][j]=-T0[i][j];
            Imda1=2/(1/r[0][M-1]+1/(0.54+0.00058*(Tf-273)));

            Imda2=2/(1/r[i][j]+1/(0.54+0.00058*(Tf-273)));

            Imda3=2/(1/r[N-1][M-1]+1/(0.54+0.00058*(Tf-273)));
         }
          if(k==1)
         {
             b[i][j]=0;
             Sp[i][j]=0;
         }
          if(k==2)
          {
             b[i][j]=2*T0[i][j]*T0[i][j]*x*y;
             r[i][j]=401;
             Imda1=401;
             Imda2=401;
             Imda3=401;
             Sp[i][j]=-T0[i][j];
          }
         if(k==3)
         {
            b[i][j]=0;
            Sp[i][j]=0;
         }
     }
 }
 for(i=0;i<N;i++)//最左边温度为 T=400
{
    T[i][0]=400;
}
 for(j=1;j<M-1;j++)//求最上面的温度
{
    aw=(1/Lr)*2/(1/r[0][j]+1/r[0][j-1]);
    ae=(1/Lr)*2/(1/r[0][j]+1/r[0][j+1]);
    as=Lr*2/(1/r[0][j]+1/r[1][j]);
    ap=(ae+aw+as-Sp[0][j]*x*y);
    T[0][j]=(aw*T[0][j-1]+ae*T[0][j+1]+as*T[1][j]+b[0][j]+q1*x)/(ap);
}
 aw=(1/Lr)*2/(1/r[0][M-1]+1/r[0][M-2]);
 as=Lr*2/(1/r[0][M-1]+1/r[1][M-1]);
 ap=(aw+as-Sp[0][M-1]*x*y);
 T[0][M-1]=(aw*T[0][M-2]+as*T[1][M-1]+b[0][M-1]+q1*x+y*Tf/(1/h+0.5*x/Imda1))/(ap+y/(1/h+0.5*x/Imda1));
 for(j=1;j<M;j++)//求内部节点温度,横向
 {
     for(i=1;i<N-1;i++)
     {
     aw=(1/Lr)*2/(1/r[i][j]+1/r[i][j-1]);
     ae=(1/Lr)*2/(1/r[i][j]+1/r[i][j+1]);
     as=Lr*2/(1/r[i][j]+1/r[i+1][j]);
     an=Lr*2/(1/r[i][j]+1/r[i-1][j]);
     ap=(ae+aw+as+an-Sp[i][j]*x*y);

     T[i][j]=(aw*T[i][j-1]+ae*T[i][j+1]+as*T[i+1][j]+an*T[i-1][j]+b[i][j])/ap;
     if(j==(M-1))
     {
        T[i][j]=(aw*T[i][j-1]+as*T[i+1][j]+an*T[i-1][j]+b[i][j]+y*Tf/(1/h+0.5*x/Imda2))/(ap-ae+y/(1/h+0.5*x/Imda2));
     }
     }
 }
 for(j=1;j<M-1;j++)//求最下面的温度
 {
     aw=(1/Lr)*2/(1/r[N-1][j]+1/r[N-1][j-1]);
     ae=(1/Lr)*2/(1/r[N-1][j]+1/r[N-1][j+1]);
     an=Lr*2/(1/r[N-1][j]+1/r[N-2][j]);
     ap=(ae+aw+an-Sp[N-1][j]*x*y);
     T[N-1][j]=(aw*T[N-1][j-1]+ae*T[N-1][j+1]+an*T[N-2][j]+b[N-1][j]+q0*x)/(ap);
 }
 aw=(1/Lr)*2/(1/r[N-1][M-1]+1/r[N-1][j-1]);
 an=Lr*2/(1/r[N-1][M-1]+1/r[N-2][M-1]);
 ap=(aw+an-Sp[N-1][M-1]*x*y);
 T[N-1][M-1]=(aw*T[N-1][N-2]+an*T[N-2][M-1]+b[N-1][M-1]+q0*x+y*Tf/(1/h+0.5*x/Imda3))/(ap+y/(1/h+0.5*x/Imda3));
 DTmax=0.0;
 for(i=0;i<N;i++)
 {
    for(j=1;j<M;j++)
    {
       DT=fabs(T[i][j]-T0[i][j]);
       T0[i][j]=T[i][j];
       if(DT>DTmax)
       {
          DTmax=DT;
       }
     }
 }
 if(DTmax>Eps)goto loop;
 if(k==0) printf("保温材料，有内热源温度分布如下：\n");
 if(k==1) printf("保温材料，无内热源温度分布如下：\n");
 if(k==2) printf("铜体材料，有内热源温度分布如下：\n");
 if(k==3) printf("铜体材料，无内热源温度分布如下：\n");
 for(i=0;i<N;i++)
 {
    for(j=0;j<M;j++)
    {
        printf("%6.3f ",T[i][j]);
    }
     printf("\n");
 }
 printf("\n");
 }
 return 0;
}

输出结果：