gk.for

有限元求解柏松方程。本文采用FORTRAN语言编制程序。程序中大部分变量采用有名公共区存储方式存储

	PROGRAM MAIN
	INTEGER N,M,NG,NP,ND,NE,II(3,2*80*80),LEE,NS,
	+        KD(81*81+1),IBE(81*81)
	DOUBLE PRECISION A(1000000),F(81*81),XY(2,81*81),
	+                 AE(81*81),DIF(2,81*81)
	COMMON /C1/N,M,NP,ND/C2/NS,KD,NG,LEE/C3/A,F/C4/DIF/C5/NE,IBE,AE
	COMMON /C6/XY,II 
	N=4
	M=N
	NG=(N+1)*(M+1)
	NP=1
	ND=3
	LEE=2*N*M
	OPEN(2,FILE='RESULTS.DAT')
	CALL IFG 
	CALL FIXP 
      CALL GKD 
	WRITE(*,*) 'NS=',NS
	WRITE(2,*) 'NS=',NS
	IF(NS.GT.1000000) THEN
	  PRINT*,'NS OVERFLOAT'
	  STOP
	ENDIF
	CALL QEB 
	CALL BDE
	CALL SOLGS 
	WRITE(2,100) (F(I),I=1,NG)
	WRITE(2,*)'DIFX====================================='
	CALL DIFF 
	WRITE(2,100) (DIF(1,I),I=1,NG)
	WRITE(2,*)'DIFY====================================='
	WRITE(2,100) (DIF(2,I),I=1,NG)
  100	FORMAT(1X,5F12.4)
	CLOSE(2)
	OPEN(9,FILE='FIG.DAT')
	WRITE(9,*)' VARIABLES= "X" "Y" "F" "DIFX" "DIFY" '
	WRITE(9,*)'ZONE I=',N+1,'J= ',M+1
	DO 110 I=1,NG
	  WRITE(9,200) XY(1,I),XY(2,I),F(I),DIF(1,I),DIF(2,I)
  110	CONTINUE
      CLOSE(9)
	OPEN(9,FILE='LINEX.DAT')
	WRITE(9,*)' VARIABLES= "X","U(Y=0)","DU/DY(Y=3.14)",
	+"DU/DY(Y=1.57)","DU/DY(Y=0.785)","DU/DY(Y=0)"  '
	DO 111 I=1,N+1
	  WRITE(9,200) XY(1,I),F(I),DIF(2,M*(N+1)+I),DIF(2,M/2*(N+1)+I),
	+               DIF(2,M/4*(N+1)+I),DIF(2,I)
  111	CONTINUE
      CLOSE(9)

	OPEN(9,FILE='LINEY.DAT')
	WRITE(9,*)' VARIABLES= "Y" "DU/DX(X=3.14)","U(X=0)","U(X=1.57)",
	+ "U(X=2.355)","U(X=3.14)"'
      DO 140 I=1,M+1
	  WRITE(9,200) XY(2,I*(N+1)),DIF(1,I*(N+1)),F((I-1)*(N+1)+1),
	+               F((I-1)*(N+1)+N/2),
     +               F((I-1)*(N+1)+3*N/4),F(I*(N+1))
  140 CONTINUE
	CLOSE(9)

  200 FORMAT(1X,6F12.4)
      END


	SUBROUTINE FIXP 

	INTEGER N,M,NG,NP,ND,NE,II(3,2*80*80),LEE,
	+        IBE(81*81),NS,KD(81*81+1)
	DOUBLE PRECISION XY(2,81*81),
	+                 AE(81*81)
	COMMON /C1/N,M,NP,ND/C2/NS,KD,NG,LEE/C5/NE,IBE,AE
	COMMON /C6/XY,II
	PRINT*,'FIXP:','N=',N,'    M=',M

	NE=M+1+N
	DO 10 I=1,N+1
	  IBE(I)=I
	  AE(I)=SIN(XY(1,IBE(I)))
   10 CONTINUE
      DO 20 I=1,M
	  J=I+N+1
	  IBE(J)=I*(N+1)+1
	  AE(J)=0.0
   20 CONTINUE

      
      END

      SUBROUTINE GK(NI,NJ,ADJ,AIJ)
	INTEGER NI,NJ
	DOUBLE PRECISION AIJ,ADJ
	IF((NI.EQ.1).AND.(NJ.EQ.1)) AIJ=1.+ADJ/6.
	IF(((NI.EQ.2).AND.(NJ.EQ.2)).OR.((NI.EQ.3).AND.(NJ.EQ.3))) 
	+   AIJ=0.5+ADJ/6.
	IF(((NI.EQ.1).AND.(NJ.NE.1)).OR.((NJ.EQ.1).AND.(NI.NE.1))) 
	+   AIJ=-0.5+ADJ/12.
	IF(((NI.EQ.2).AND.(NJ.EQ.3)).OR.((NJ.EQ.2).AND.(NI.EQ.3))) 
	+   AIJ=0.+ADJ/12.
	END

	SUBROUTINE GF(NI,N,M,LE,YI,FE)
	INTEGER NI,LE,N,M
	DOUBLE PRECISION FE,YI
	PARAMETER(PI=3.1415926)
	HY=PI/M
	FE=0.
	DO 10 J=1,M
	  IF(LE.EQ.2*N*J) THEN
	    IF(NI.EQ.3) FE=SIN(YI+HY)+(COS(YI+HY)-COS(YI))/HY+SIN(YI)-
     +                   SIN(YI+HY)
	    IF(NI.EQ.1) FE=-SIN(YI)-(COS(YI)-COS(YI-HY))/HY
	  END IF
   10	CONTINUE
      END
	
      SUBROUTINE AK(I,J,AIJ)
	INTEGER I,J
	DOUBLE PRECISION AIJ

	INTEGER NS,KD(81*81+1)
	DOUBLE PRECISION A(1000000),F(81*81)
	COMMON /C2/NS,KD,NG,LEE/C3/A,F

	IF(I.GE.J) THEN
	  IGP=KD(I+1)-I+J   
	ELSE
	  IGP=KD(J+1)-J+I
	END IF
	A(IGP)=A(IGP)+AIJ
	END

	SUBROUTINE QEB 
 	DOUBLE PRECISION ADJ,AIJ,FE

	INTEGER N,M,NG,NP,ND,II(3,2*80*80),LEE,NS,
	+        KD(81*81+1)
	DOUBLE PRECISION A(1000000),F(81*81),XY(2,81*81)
	COMMON /C1/N,M,NP,ND/C2/NS,KD,NG,LEE/C3/A,F
	COMMON /C6/XY,II
	PARAMETER(PI=3.1415926)
	PRINT*,'QEB:','N=',N,'    M=',M
	ADJ=PI/N*PI/M
	DO 10 LP=1,NG
   10	  F(LP)=0.
 	DO 20 LP=1,NS
   20	  A(LP)=0.
      DO 30 LE=1,LEE
	  DO 40 NI=1,ND
	    I=II(NI,LE)
	    CALL GF(NI,N,M,LE,XY(2,I),FE)
	    F(I)=F(I)+FE
	    DO 50 NJ=1,NI
	      J=II(NJ,LE)
	      CALL GK(NI,NJ,ADJ,AIJ)
	      CALL AK(I,J,AIJ)
   50	    CONTINUE
   40   CONTINUE
   30 CONTINUE
      END


	SUBROUTINE SOLGS 
	DOUBLE PRECISION RES0(81*81),RES(81*81)

	INTEGER NG,LEE,NS,KD(81*81+1)
	DOUBLE PRECISION A(1000000),F(81*81)
	COMMON /C2/NS,KD,NG,LEE/C3/A,F

	NTIMES=0
	DO 11 I=1,NG
	  RES0(I)=0.
	  RES(I)=0.
   11 CONTINUE
	DO 100 K=1,5000 
	NTIMES=NTIMES+1
      DO 10 I=1,NG
	  SUM=0.
	  IIG=KD(I+1)-I
	  MI=KD(I)-IIG+1
	  DO 20 J=MI,I-1
	    IGP=IIG+J
	  	SUM=SUM+A(IGP)*RES(J)
   20   CONTINUE
        DO 30 J=I+1,NG  
		JJG=KD(J+1)-J
	    IGP=JJG+I
	    MJ=KD(J)-JJG+1
	    IF(I.GE.MJ) SUM=SUM+A(IGP)*RES(J)
   30   CONTINUE
        RES(I)=(F(I)-SUM)/A(KD(I+1))
	  IF(ERR.LT.ABS(RES(I)-RES0(I))) ERR=ABS(RES(I)-RES0(I))
	  RES0(I)=RES(I)
   10   CONTINUE
	  WRITE(*,*) 'TIMES=',NTIMES,'   ERR=',ERR
        IF(ERR.LT.1.E-6) THEN
	    DO 40 I=1,NG
	      F(I)=RES(I)
   40     CONTINUE
          PRINT*,'NTIMES=',NTIMES
          RETURN
        ENDIF
	  ERR=0.
  100	CONTINUE
      DO 110 I=1,NG
	      F(I)=RES(I)
  110 CONTINUE

	END

	SUBROUTINE BDE 

	INTEGER NG,NE,LEE,NS,
	+        KD(81*81+1),IBE(81*81)
	DOUBLE PRECISION A(1000000),F(81*81),
	+                 AE(81*81)
	COMMON /C2/NS,KD,NG,LEE/C3/A,F/C5/NE,IBE,AE
	DO 10 L=1,NE
	  I=IBE(L)
	  X0=AE(L)
	  IIG=KD(I+1)-I
	  MI=KD(I)-IIG+1
	  DO 20 J=MI,I-1
	    IGP=IIG+J
	    F(J)=F(J)-X0*A(IGP)
	    A(IGP)=0.
   20   CONTINUE
        A(KD(I+1))=1.
	  F(I)=X0
	  DO 30 J=I+1,NG
	  JJG=KD(J+1)-J
	  IGP=I+JJG
	  MJ=KD(J)-JJG+1
	  IF(I.GE.MJ) THEN
	    F(J)=F(J)-A(IGP)*X0
	    A(IGP)=0.
	  END IF
   30   CONTINUE
   10 CONTINUE
      END



	SUBROUTINE GKD 
	INTEGER N,M,NG,NP,ND,II(3,2*80*80),LEE,NS,
	+        KD(81*81+1)
	DOUBLE PRECISION XY(2,81*81)
	COMMON /C1/N,M,NP,ND/C2/NS,KD,NG,LEE
	COMMON /C6/XY,II
	KD(1)=0
	PRINT*,'NG=',NG
	DO 10, L=1,NG
	   MIN=NG
	   DO 20, LE=1,LEE
	     DO 30, J=1,ND
	       IF(L.EQ.II(J,LE)) GOTO 40
   30	     CONTINUE
           GOTO 70
   40      DO 50, K=1,ND
             IF(II(K,LE).LT.MIN) MIN=II(K,LE)
   50      CONTINUE
   70		 CONTINUE
   20   CONTINUE
        LI=NP*(L-MIN)
	  I=NP*(L-1)
	  DO 60, NN=1,NP
	    KD(I+NN+1)=KD(I+NN)+LI+NN
   60   CONTINUE
   10 CONTINUE
      NS=KD(NG*NP+1)
	END


	SUBROUTINE UDIFF(NI,NFLAG,UDIF,LE,ADJ)
	INTEGER NI,NFLAG,LE
	DOUBLE PRECISION UDIF,ADJ
	INTEGER II(3,2*80*80)
	DOUBLE PRECISION XY(2,81*81)
	COMMON /C6/XY,II 

	IF(NFLAG.EQ.1) THEN
	  IF(NI.EQ.1) UDIF=(XY(2,II(2,LE))-XY(2,II(3,LE)))/ADJ
	  IF(NI.EQ.2) UDIF=(XY(2,II(3,LE))-XY(2,II(1,LE)))/ADJ
	  IF(NI.EQ.3) UDIF=(XY(2,II(1,LE))-XY(2,II(2,LE)))/ADJ
	ENDIF
	IF(NFLAG.EQ.2) THEN
	  IF(NI.EQ.1) UDIF=(XY(1,II(3,LE))-XY(1,II(2,LE)))/ADJ
	  IF(NI.EQ.2) UDIF=(XY(1,II(1,LE))-XY(1,II(3,LE)))/ADJ
	  IF(NI.EQ.3) UDIF=(XY(1,II(2,LE))-XY(1,II(1,LE)))/ADJ
	ENDIF
	END

      SUBROUTINE DIFF 
      INTEGER NFLAG
	DOUBLE PRECISION TADJ(81*81),UDIF,ADJ

	INTEGER N,M,NG,NP,ND,II(3,2*80*80),LEE,NS,
	+        KD(81*81+1)
	DOUBLE PRECISION A(1000000),F(81*81),XY(2,81*81),
	+                 DIF(2,81*81)
	COMMON /C1/N,M,NP,ND/C2/NS,KD,NG,LEE/C3/A,F/C4/DIF
	COMMON /C6/XY,II 

	PARAMETER(PI=3.1415926)
	PRINT*,'DIFF:','N=',N,'    M=',M

	ADJ=PI/N*PI/M
	DO 10, I=1,NG
	  DIF(1,I)=0.
	  DIF(2,I)=0.
	  TADJ(I)=0.
   10 CONTINUE
      DO 20, LE=1,LEE
	  DO 30, NI=1,ND
	    DO 40, NFLAG=1,2
	      CALL UDIFF(NI,NFLAG,UDIF,LE,ADJ)
	      DIF(NFLAG,II(1,LE))=DIF(NFLAG,II(1,LE))+
     +                           UDIF*F(II(NI,LE))*ADJ
	      DIF(NFLAG,II(2,LE))=DIF(NFLAG,II(2,LE))+
     +                           UDIF*F(II(NI,LE))*ADJ
	      DIF(NFLAG,II(3,LE))=DIF(NFLAG,II(3,LE))+
     +                           UDIF*F(II(NI,LE))*ADJ

   40     CONTINUE
	    TADJ(II(NI,LE))=TADJ(II(NI,LE))+ADJ

   30   CONTINUE
   20 CONTINUE
      
	DO 50, I=1,NG
	   IF(TADJ(I).GT.0.) THEN
	     DIF(1,I)=DIF(1,I)/TADJ(I)
           DIF(2,I)=DIF(2,I)/TADJ(I)
	   END IF
   50 CONTINUE
      END

      SUBROUTINE IFG 
	PARAMETER (PI=3.1415926)
	DOUBLE PRECISION HX,HY

	INTEGER N,M,NG,NP,ND,II(3,2*80*80),LEE,NS,KD(81*81+1)
	DOUBLE PRECISION XY(2,81*81)
	COMMON /C1/N,M,NP,ND/C2/NS,KD,NG,LEE
	COMMON /C6/XY,II 
	PRINT*,'IFG:','N=',N,'    M=',M
	HX=PI/N
	HY=PI/M
	DO 10, J=1,M+1
	  DO 20, I=1,N+1
	    NG=(J-1)*(N+1)+I
	    XY(1,NG)=(I-1)*HX
	    XY(2,NG)=(J-1)*HY
   20   CONTINUE
   10 CONTINUE
	LE=0
	DO 30, J=1,M
	  NJ=(J-1)*(N+1)
	  DO 40, I=1,N
	    LE=LE+1
	    II(1,LE)=NJ+I
	    II(2,LE)=NJ+I+1
	    II(3,LE)=NJ+N+1+I
	    LE=LE+1
	    II(1,LE)=NJ+N+1+I+1
	    II(2,LE)=NJ+N+1+I
	    II(3,LE)=NJ+I+1
   40   CONTINUE
   30 CONTINUE
      print*,'LEE=',LEE
	END