C Solution of Poisson's equation -△u = f in the rectangle
C [XMIN,XMAX]×[YMIN,YMAX]  by Fourier series
C 近似精度は Fourier 展開係数の項数 KOSU, および, 
C 長方形上の定積分を４点公式で計算するときの辺分割数 IIB に依存する
C 章末課題 12.1 に対応するためグラフィック表示をやめ, 
C 積分の辺分割数は 320 に固定して真値に替え，これから得られた
C Fourier 係数を用いて N=10, 20, 40, 80, 160 に対する部分和の値に
C Richardson 補外法を適用して，無限和の真値を予測し, かつ
C 部分和の誤差の 1/N に関するオーダーを見るように変更したもの
C 値の比較は中央値だけで行っている
      PROGRAM POISN1
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C      PARAMETER(KOSUMX=160,IIBCT=6,INTMES=320)
      PARAMETER(KOSUMX=320,IIBCT=5,INTMES=640)
C 時間が掛かりすぎると思う人は, 上の行の代わりに次を用いよ
C     PARAMETER(KOSUMX=80,IIBCT=4,INTMES=160)
      PARAMETER(XMIN=0,XMAX=1,YMIN=0,YMAX=1)
      REAL*8 CF(KOSUMX,KOSUMX),T(0:IIBCT-1,IIBCT)
      CALL FRCOEF(CF,KOSUMX,XMIN,XMAX,YMIN,YMAX,INTMES)
      KOSU=10
      DO 100 I=1,IIBCT
          T(0,I)=FFOUR(0.5D0,0.5D0,CF,KOSU,KOSUMX)
          KOSU=KOSU*2
  100 CONTINUE
      M=2
      DO 300 I=1,IIBCT-1
         DO 200 I2=1,IIBCT-I
            T(I,I2)=(M*T(I-1,I2+1)-T(I-1,I2))/(M-1)
  200    CONTINUE
         M=M*2
  300 CONTINUE
      WRITE(*,*)'Richardson extrapolation of 1st order:'
      DO 500 I=0,IIBCT-1
          WRITE(*,'(5F18.15)')(T(I-I2+1,I2),I2=1,I+1)
  500 CONTINUE
      WRITE(*,*)'Error:'
      WRITE(*,'(4F19.15)')((T(0,I)-T(IIBCT-1,1)),I=1,IIBCT)
      WRITE(*,*)'Ratio:'
      WRITE(*,'(4F19.15)')((T(0,I+1)-T(IIBCT-1,1))/
     +            (T(0,I)-T(IIBCT-1,1)),I=1,IIBCT-1)
      M=4
      DO 800 I=1,IIBCT-1
         DO 700 I2=1,IIBCT-I
            T(I,I2)=(M*T(I-1,I2+1)-T(I-1,I2))/(M-1)
  700    CONTINUE
         M=M*4
  800 CONTINUE
      WRITE(*,*)'Richardson extrapolation of 2nd order:'
      DO 1000 I=0,IIBCT-1
          WRITE(*,'(5F18.15)')(T(I-I2+1,I2),I2=1,I+1)
 1000 CONTINUE
      WRITE(*,*)'Error:'
      WRITE(*,'(4F19.15)')((T(0,I)-T(IIBCT-1,1)),I=1,IIBCT)
      WRITE(*,*)'Ratio:'
      WRITE(*,'(4F19.15)')((T(0,I+1)-T(IIBCT-1,1))/
     +              (T(0,I)-T(IIBCT-1,1)),I=1,IIBCT-1)
      M=8
      DO 3000 I=1,IIBCT-1
         DO 2000 I2=1,IIBCT-I
            T(I,I2)=(M*T(I-1,I2+1)-T(I-1,I2))/(M-1)
 2000    CONTINUE
         M=M*8
 3000 CONTINUE
      WRITE(*,*)'Richardson extrapolation of 3rd order:'
      DO 4000 I=0,IIBCT-1
          WRITE(*,'(5F18.15)')(T(I-I2+1,I2),I2=1,I+1)
 4000 CONTINUE
      WRITE(*,*)'Error:'
      WRITE(*,'(4F19.15)')((T(0,I)-T(IIBCT-1,1)),I=1,IIBCT)
      WRITE(*,*)'Ratio:'
      WRITE(*,'(4F19.15)')((T(0,I+1)-T(IIBCT-1,1))/
     +               (T(0,I)-T(IIBCT-1,1)),I=1,IIBCT-1)
      M=16
      DO 8000 I=1,IIBCT-1
         DO 7000 I2=1,IIBCT-I
            T(I,I2)=(M*T(I-1,I2+1)-T(I-1,I2))/(M-1)
 7000    CONTINUE
         M=M*16
 8000 CONTINUE
      WRITE(*,*)'Richardson extrapolation of 4th order:'
      DO 9000 I=0,IIBCT-1
          WRITE(*,'(5F18.15)')(T(I-I2+1,I2),I2=1,I+1)
 9000 CONTINUE
      WRITE(*,*)'Error:'
      WRITE(*,'(4F19.15)')((T(0,I)-T(IIBCT-1,1)),I=1,IIBCT)
      WRITE(*,*)'Ratio:'
      WRITE(*,'(4F19.15)')((T(0,I+1)-T(IIBCT-1,1))/
     +               (T(0,I)-T(IIBCT-1,1)),I=1,IIBCT-1)
      END
C
C 方程式の右辺の函数のフーリエ展開の係数を計算し配列CFに格納
      SUBROUTINE FRCOEF(CF,KOSUMX,XMIN,XMAX,YMIN,YMAX,IIB)
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER(PI=3.14159265358979D0)
      REAL*8 CF(KOSUMX,KOSUMX)
      EXTERNAL G
      CC=4/(XMAX-XMIN)/(YMAX-YMIN)/PI/PI
      XX=(XMAX-XMIN)**2
      YY=(YMAX-YMIN)**2
      DO 200 I=1,KOSUMX
         DO 100 J=1,KOSUMX+1-I
            CF(I,J)=SQINT(G,I,J,XMIN,XMAX,YMIN,YMAX,IIB)
            CF(I,J)=CF(I,J)*CC/(I*I/XX+J*J/YY)
C            IF (DABS(CF(I,J)).LT.1.0D-8) CF(I,J)=0.0D0
  100    CONTINUE
  200 CONTINUE
      RETURN
      END
C
C 配列CFに格納されたデータを係数とするフーリエ級数の和を計算する
      FUNCTION FFOUR(X,Y,CF,KOSU,KOSUMX)
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER(PI=3.1415926535897932384D0)
      REAL*8 CF(KOSUMX,KOSUMX)
      FFOUR=0
      DO 200 I=1,KOSU
         DO 100 J=1,KOSU+1-I
C            IF (CF(I,J).NE.0.0D0) THEN
               FFOUR=FFOUR+CF(I,J)*DSIN(I*PI*X)*DSIN(J*PI*Y)
C            END IF
  100    CONTINUE
  200 CONTINUE
      RETURN
      END
C
C 函数F(X,Y)の長方形[XMIN,XMAX]×[YMIN,YMAX]上の二次元定積分を
C 四点公式で計算
      FUNCTION SQINT(F,I,J,XMIN,XMAX,YMIN,YMAX,IIB)
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER(GP=1.732050807568877/3)
      EXTERNAL F
      HX=(XMAX-XMIN)/IIB
      HY=(YMAX-YMIN)/IIB
      GX1=HX*(1-GP)/2
      GX2=HX*(1+GP)/2
      GY1=HY*(1-GP)/2
      GY2=HY*(1+GP)/2
      SQINT=0
      X=XMIN
      DO 70 K=1,IIB
         Y=YMIN
         DO 50 L=1,IIB
            F1=F(X+GX1,Y+GY1,I,J)
            F2=F(X+GX1,Y+GY2,I,J)
            F3=F(X+GX2,Y+GY1,I,J)
            F4=F(X+GX2,Y+GY2,I,J)
            SQINT=SQINT+F1+F2+F3+F4
            Y=Y+HY
   50    CONTINUE
         X=X+HX
   70 CONTINUE
      SQINT=SQINT*HX*HY/4
      RETURN
      END
C
C 方程式の右辺の函数のフーリエ展開の係数を計算するための被積分函数
      FUNCTION G(X,Y,I,J)
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER(XMIN=0,XMAX=1,YMIN=0,YMAX=1,PI=3.14159265358979323D0,
     +   C1=PI/(XMAX-XMIN),C2=PI/(YMAX-YMIN))
      G=F(X,Y)*DSIN(C1*I*(X-XMIN))*DSIN(C2*J*(Y-YMIN))
      RETURN
      END
C
C GIVEN RIGHT-HAND SIDE
      FUNCTION F(X,Y)
      REAL*8 F,X,Y
C      IF ((DABS(X-0.5).LT.0.2D0).AND.(DABS(Y-0.5).LT.0.2D0)) THEN
C       F=5
C      ELSE
C       F=0
C      END IF
      F=5
      RETURN
      END