p-r.for

fortran编写

      FUNCTION InitialF(alpha,beta,pi,x,y)
      implicit double precision(a-z)
      InitialF=sin(pi*alpha*x)*sin(beta*pi*y)
      END
      FUNCTION BoundX0(y,t)
      implicit double precision(a-z)
      BoundX0=0.0D0
      END
      FUNCTION BoundX1(y,t)
      implicit double precision(a-z)
      BoundX1=0.0D0
      END
      FUNCTION BoundY0(x,t)
      implicit double precision(a-z)
      BoundY0=0.0D0
      END
      FUNCTION BoundY1(x,t)
      implicit double precision(a-z)
      BoundY1=0.0D0
      END

      include 'fgraph.fi'
      implicit double precision (A-Z)
      integer j,l,DJ,DN
      parameter (DJ=20,T=1.0D0,DN=100,a=1.0D0,Width=1.0d0)
      parameter(alpha=1.0D0,beta=1.0D0,Pi=3.14159265358979D0)
      dimension u(0:DJ,0:DJ),u1(0:DJ,0:DJ),utrue(0:DJ,0:DJ)
      dimension a1(2:DJ-1),b1(1:DJ-1),c1(1:DJ-2),d1(1:DJ-1)
      !call clearscreen(2)
      h=Width/DJ
      Tao=T/DN
      r=a*Tao/h**2
      tt=0.0
      call Initial(alpha,beta,pi,u,DJ,h)
      do n=1,DN,1

         tt=(n-0.5)*Tao
         do l=1,DJ-1,1
             y=l*h
             tmp1=BoundX0(y,tt)
             tmp2=BoundX1(y,tt)
             call Initial1(a1,b1,c1,d1,DJ,r,u,l,tmp1,tmp2)
		
             call ZhuiGanFa(a1,b1,c1,d1,DJ)
             do j=1,DJ-1,1
                u1(j,l)=d1(j)
             end do
             u1(0,l)=tmp1
             u1(DJ,l)=tmp2
         end do
         do j=0,DJ,1
             x=j*h
             u1(j,0)=BoundY0(x,tt)
             u1(j,DJ)=BoundY1(x,tt)
         end do

         tt=n*Tao
         do j=1,DJ-1,1
             x=j*h
             tmp1=BoundY0(x,tt)
             tmp2=BoundY1(x,tt)

             call Initial2(a1,b1,c1,d1,DJ,r,u1,j,tmp1,tmp2)
             call ZHuiGanFa(a1,b1,c1,d1,DJ)
	     do l=1,DJ-1,1
		u(j,l)=d1(l)
             end do
             u(j,0)=tmp1
             u(j,DJ)=tmp2
         end do
         do l=0,DJ,1
             y=l*h
             u(0,l)=BoundX0(y,tt)
             u(DJ,l)=BoundX1(y,tt)
         end do
      call TrueValue(alpha,beta,pi,utrue,DJ,h,TT)
      call OutPut2(u,utrue,DJ)
        
      end do

	end

      subroutine Initial(alpha,beta,pi,u,DJ,h)
      implicit double precision (a-z)
      integer DJ,j,l
      dimension u(0:DJ,0:DJ)
      do l=0,DJ,1
          y=l*h
          do j=0,DJ,1
              x=j*h
              u(j,l)=InitialF(alpha,beta,pi,x,y)
          end do
      end do
      end

      subroutine Initial1(a,b,c,d,DJ,r,u,l,tmp1,tmp2)
      implicit double precision (a-z)
      integer DJ,j,l
      dimension a(2:DJ-1),b(1:DJ-1),c(1:DJ-2),d(1:DJ-1)
      dimension u(0:DJ,0:DJ)
      c(1)=-r/2.0
      b(1)=1+r
      j=1
      d(j)=(1.0-r)*u(j,l)-c(1)*(u(j,l+1)+u(j,l-1)+tmp1)
      do j=2,DJ-1,1
         a(j)=c(1)
         c(j-1)=c(1)
         b(j)=b(1)
         d(j)=(1.0-r)*u(j,l)-c(1)*(u(j,l+1)+u(j,l-1))
      end do
      j=DJ-1
cc    d(j)=(1.0-r)*u(j,l)-c(1)*(u(j+1,l)+u(j-1,l)+u1(j+1,l))
      d(j)=d(j)-c(1)*tmp2
      end

      subroutine Initial2(a,b,c,d,DJ,r,u1,j,tmp1,tmp2)
      implicit double precision (a-z)
      integer DJ,l,j
      dimension a(2:DJ-1),b(1:DJ-1),c(1:DJ-2),d(1:DJ-1)
      dimension u1(0:DJ,0:DJ)
      b(1)=1.0+r
      c(1)=-r/2.0
      l=1
      d(l)=(1.0-r)*u1(j,1)-c(1)*(u1(j+1,l)+u1(j-1,l)+tmp1)
      do l=2,DJ-1,1
         b(l)=b(1)
         c(l-1)=c(1)
         a(l)=c(1)
         d(l)=(1.0-r)*u1(j,l)-c(1)*(u1(j+1,l)+u1(j-1,l))
      end do
      l=DJ-1
      d(l)=d(l)-c(1)*tmp2
      end

      subroutine TrueValue(alpha,beta,pi,utrue,DJ,h,tt)
      implicit double precision (a-z)
      integer DJ,j,l
      dimension utrue(0:DJ,0:DJ)
      tmp1=exp(-pi**2*(alpha**2+beta**2)*tt)
      tmp2=alpha*pi
      tmp3=beta*pi
      do l=0,DJ,1
         y=l*h
         do j=0,DJ,1
            x=j*h
            utrue(j,l)=tmp1*sin(tmp2*x)*sin(tmp3*y)
         end do
      end do
      end

      subroutine Output2(u,utrue,DJ)
      integer DJ,j,l
      double precision u(0:DJ,0:DJ),utrue(0:DJ,0:DJ)
      do l=0,DJ,1
         do j=0,DJ,1
             write(*,111)u(j,l),utrue(j,l),abs(u(j,l)-utrue(j,l))
         end do
         print *,'Press Enter Key to continue...' ! Please notice here!
         read (*,*)
      end do
      write(*,*)char(7)
111   format(1x,3F25.10,F10.6)
      end

      SUBROUTINE ZhuiGanFa (A, B, C, D, DJ)
         IMPLICIT DOUBLE PRECISION (A-Z)
         INTEGER DJ
         DOUBLE PRECISION A(2:DJ-1),B(1:DJ-1),C(1:DJ-2),D(1:DJ-1),T
         N=DJ-1
         C(1) = C(1) / B(1)
         D(1) = D(1) / B(1)  !XiaBiao BiXu Shi 1
         DO I = 2 , N-1,1
            T = B(I) - C(I - 1) * A(I)
            C(I) = C(I) / T
            D(I) = (D(I) - D(I - 1) * A(I)) / T
         END DO
         D(N)= (D(N)-D(N-1)*A(N))/(B(N)-C(N-1)*A(N))
         DO I = N-1 , 1 , -1
            D(I) = D(I) - C(I) * D(I + 1)
         END DO
         !The D(I) return with results
         !DO I = 1 , DJ , 1
         !   X(I) = D(I)
         !END DO
      END