numeranalf.f

二维热能方程的求解

c=======================================================================
      subroutine numeranalf(nt,nx,ny,level,at,bt,ax,bx,ay,by,x,ierr)
c=======================================================================
      implicit none
      integer nt,nx,ny,level,ierr
      real*8 at,bt,ax,bx,ay,by
      real*8 x(nx,ny)

      integer ix,iy
      real*8 err,top,bot,top8,bot8,diff
      real*8 hx,hy,px,py,tsol,xsol
      real*8 truesol

      if(level.eq.0) return

      if(level.ge.1) print'("NUMERANALF: @ t=",f7.3)',bt
      if(level.ge.3) print*,"nx=",nx," ny=",ny

      err=0.
      top=0.
      bot=0.
      top8=0.
      bot8=0.

      hx=(bx-ax)/(nx-1)
      hy=(by-ay)/(ny-1)

      do iy=1,ny
         py=hy*(iy-1)+ay
      do ix=1,nx
         px=hx*(ix-1)+ax
         tsol=truesol(bt,px,py)
         diff=(tsol-x(ix,iy))
         top=top+diff**2
         bot=bot+tsol**2
         top8=max(top8,dabs(diff))
         bot8=max(bot8,dabs(tsol))
      enddo
      enddo

      if( (sqrt(bot)/nx).ge.1.e-8) then
         err=dsqrt(top/bot)
         print'(" Relative-L2-Error=",1pe8.2)',err
         print'(" Relative-L8-Error=",1pe8.2)',top8/bot8
      endif
      err=sqrt(top)/nx
      print'(" L2-Error=",1pe8.2,2x,"L8-Error=",1pe8.2)',err,top8

      return
      end