numeranal.c

二维热能方程的求解

#include <stdio.h>
#include <math.h>
typedef  double ATYPE;
typedef  double VTYPE;
typedef  float DTYPE;

void numeranal(nt,nx,ny,level,at,bt,ax,bx,ay,by,x,ierr)

int nt,nx,ny,level;
VTYPE *at,*bt,*ax,*bx,*ay,*by;
ATYPE x[][nx];
int *ierr;

{
   int ix,iy;
   VTYPE err,top,bot;
   VTYPE hx,hy,px,py,tsol,xsol;

   #if defined (_IBMR2) || defined (AIX) || defined (hpux)
       ATYPE truesol();
   #else
       ATYPE truesol_();
   #endif

   if(level==0) return;

   if(level>=1) printf("NUMERANAL: @ t=%.3g\n",*bt);
   if(level>=3) printf(" nt=%d nx=%d ny=%d\n",nt,nx,ny);

   err=top=bot=0.;
   hx=(*bx-*ax)/(nx-1);
   hy=(*by-*ay)/(ny-1);

   for(iy=0;iy<ny;iy++){
      py=hy*iy+*ay;
      for(ix=0;ix<nx;ix++){
         px=hx*ix+*ax;
         #if defined (_IBMR2) || defined (AIX) || defined (hpux)
            tsol=truesol(bt,&px,&py);
         #else
            tsol=truesol_(bt,&px,&py);
         #endif
         xsol=x[iy][ix];
         top+=(xsol-tsol)*(xsol-tsol);
         bot+=tsol*tsol;
       }
    }

   if( (sqrt(bot)/nx)>=1.e-8){
      err=sqrt(top/bot);
      if(err<0.99) printf(" Relative-L2-Error=%.3g\n",err);
    }
   err=sqrt(top)/nx;
   printf(" L2-Error=%.3g\n",err);
}