C FIGURE 4.7.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),X(0:N),Y(0:N) DX = 1.D0/N DY = DX C TRY VARIOUS VALUES OF THE PARAMETER OMEGA DO 35 IW=0,19 W = 1.0 + IW*0.05 C SET UNKNOWNS TO ZERO STARTING VALUES DO 5 I=1,N-1 DO 5 J=1,N-1 U(I,J) = 0.0 5 CONTINUE C SET VALUES OF U 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 U(L,N) = X(L) U(0,L) = 0.0 U(N,L) = Y(L)**3 10 CONTINUE C BEGIN SOR ITERATION DO 20 ITER=1,10000 ERMAX = 0.0 C UPDATE UNKNOWNS USING SOR FORMULA DO 15 I=1,N-1 DO 15 J=1,N-1 U(I,J) = U(I,J) - W*( (-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) ) / (2/DX**2+2/DY**2+A) ERR = ABS(U(I,J) - X(I)*Y(J)**3) ERMAX = MAX(ERMAX,ERR) 15 CONTINUE MAXITR = ITER C STOP IF ERROR < 10**(-6) IF (ERMAX.LT.1.D-6) GO TO 25 20 CONTINUE 25 PRINT 30, W,MAXITR 30 FORMAT (5X,F8.2,I10) 35 CONTINUE STOP END