C FIGURE 4.8.1 IMPLICIT DOUBLE PRECISION(A-H,O-Z) C N = NUMBER OF POINTS IN X,Y DIRECTIONS C A = PDE PARAMETER PARAMETER (N=40, A=40.0) DIMENSION U(0:N,0:N),R(0:N,0:N),P(0:N,0:N), & AP(0:N,0:N),X(0:N),Y(0:N) DOUBLE PRECISION LAMBDA DX = 1.D0/N DY = DX C U = 0 DO 5 I=1,N-1 DO 5 J=1,N-1 U(I,J) = 0.0 5 CONTINUE C SET VALUES OF U,P ON BOUNDARY (THESE WILL C NEVER CHANGE) DO 10 L=0,N X(L) = L*DX Y(L) = L*DY U(L,0) = 0.0 P(L,0) = 0.0 U(L,N) = X(L) P(L,N) = 0.0 U(0,L) = 0.0 P(0,L) = 0.0 U(N,L) = Y(L)**3 P(N,L) = 0.0 10 CONTINUE C R = B - A*U C P = R DO 15 I=1,N-1 DO 15 J=1,N-1 R(I,J) = (U(I+1,J)-2*U(I,J)+U(I-1,J))/DX**2 & +(U(I,J+1)-2*U(I,J)+U(I,J-1))/DY**2 & - A*U(I,J) + X(I)*Y(J)*(A*Y(J)**2-6) P(I,J) = R(I,J) 15 CONTINUE C BEGIN CONJUGATE GRADIENT ITERATION DO 30 ITER=1,1000 C AP = A*P DO 20 I=1,N-1 DO 20 J=1,N-1 AP(I,J) = (-P(I+1,J)+2*P(I,J)-P(I-1,J))/DX**2 & +(-P(I,J+1)+2*P(I,J)-P(I,J-1))/DY**2 & + A*P(I,J) 20 CONTINUE C LAMBDA = R*P / P*AP LAMBDA = DOT(R,P,N)/DOT(P,AP,N) C U = U + LAMBDA*P CALL UPDATE(U,U,LAMBDA,P,N) C R = R - LAMBDA*AP CALL UPDATE(R,R,-LAMBDA,AP,N) C ALPHA = -R*AP / P*AP ALPHA = -DOT(R,AP,N)/DOT(P,AP,N) C P = R + ALPHA*P CALL UPDATE(P,R,ALPHA,P,N) C CALCULATE MAXIMUM ERROR ERMAX = 0.0 DO 25 I=1,N-1 DO 25 J=1,N-1 ERR = ABS(U(I,J) - X(I)*Y(J)**3) ERMAX = MAX(ERMAX,ERR) 25 CONTINUE MAXITR = ITER C STOP IF ERROR < 10**(-6) IF (ERMAX.LT.1.D-6) GO TO 35 30 CONTINUE 35 PRINT 40, MAXITR 40 FORMAT (I10) STOP END FUNCTION DOT(X,Y,N) IMPLICIT DOUBLE PRECISION(A-H,O-Z) DIMENSION X(0:N,0:N),Y(0:N,0:N) C CALCULATE DOT PRODUCT OF X AND Y DOT = 0.0 DO 5 I=1,N-1 DO 5 J=1,N-1 DOT = DOT + X(I,J)*Y(I,J) 5 CONTINUE RETURN END SUBROUTINE UPDATE(X,Y,S,Z,N) IMPLICIT DOUBLE PRECISION(A-H,O-Z) DIMENSION X(0:N,0:N),Y(0:N,0:N),Z(0:N,0:N) C CALCULATE X = Y + S*Z DO 5 I=1,N-1 DO 5 J=1,N-1 X(I,J) = Y(I,J) + S*Z(I,J) 5 CONTINUE RETURN END