C*********************************************************************** C * C EXTRAPOLATION METHOD * C * C*********************************************************************** C C C C TO APPROXIMTE THE SOLUTION OF THE INITIAL VALUE PROBLEM: C Y'=F(T,Y), A<=T<=B, Y(A)=ALPHA, C WITH LOCAL ERROR WITHIN A GIVEN TOLERANCE: C C INPUT: ENDPOINTS A,B; INITIAL CONDITION ALPHA; C TOLERANCE TOL; MAXIMUM STEPSIZE HMAX; C MINIMUM STEPSIZE HIMIN. C C OUTPUT: T,W,H WHERE W APPROXIMATES Y(T) AND STEPSIZE H WAS C USED OR A MESSAGE THAT MINIMUM STEPSIZE EXCEEDED. C DIMENSION NK(8),Y(8),Q(7,7) CHARACTER NAME*14,NAME1*14,AA*1 INTEGER INP,OUP,FLAG LOGICAL OK F(TZ,WZ)= WZ-TZ*TZ+1 C STEP 1 DATA NK/2,3,4,6,8,12,16,24/ OPEN(UNIT=5,FILE='CON',ACCESS='SEQUENTIAL') OPEN(UNIT=6,FILE='CON',ACCESS='SEQUENTIAL') WRITE(6,*) 'This is the Extrapolation Method.' WRITE(6,*) 'Has the function F been created in the program? ' WRITE(6,*) 'Enter Y or N ' WRITE(6,*) ' ' READ(5,*) AA IF(( AA .EQ. 'Y' ) .OR. ( AA .EQ. 'y' )) THEN OK = .FALSE. 10 IF (OK) GOTO 11 WRITE(6,*) 'Input left and right endpoints separated by' WRITE(6,*) 'blank - include decimal point' WRITE(6,*) ' ' READ(5,*) A, B IF (A.GE.B) THEN WRITE(6,*) 'Left endpoint must be less' WRITE(6,*) 'than right endpoint' ELSE OK = .TRUE. ENDIF GOTO 10 11 OK = .FALSE. 13 WRITE(6,*) 'Input the initial condition.' WRITE(6,*) ' ' READ(5,*) ALPHA OK = .FALSE. 14 IF (OK) GOTO 15 WRITE(6,*) 'Input tolerance ' WRITE(6,*) ' ' READ(5,*) TOL IF (TOL.LE.0.0) THEN WRITE(6,*) 'Tolerance must be positive ' ELSE OK = .TRUE. ENDIF GOTO 14 15 OK = .FALSE. 16 IF (OK) GOTO 17 WRITE(6,*) 'Input minimum and maximum mesh spacing' WRITE(6,*) 'separated by blank.' WRITE(6,*) ' ' READ(5,*) HMIN, HMAX IF ((HMIN.GE.HMAX) .OR. (HMIN.LE.0.0)) THEN WRITE(6,*) 'Minimum mesh spacing must be positive ' WRITE(6,*) 'and less than maximum mesh spacing.' ELSE OK = .TRUE. ENDIF GOTO 16 17 CONTINUE ELSE WRITE(6,*) 'The program will end so that the function F ' WRITE(6,*) 'can be created ' OK = .FALSE. ENDIF IF(.NOT.OK) GOTO 700 WRITE(6,*) 'Select output destinations: ' WRITE(6,*) '1. Screen ' WRITE(6,*) '2. Text file ' WRITE(6,*) 'Enter 1 or 2 ' WRITE(6,*) ' ' READ(5,*) FLAG IF ( FLAG .EQ. 2 ) THEN WRITE(6,*) 'Input the file name in the form - ' WRITE(6,*) 'drive:name.ext' WRITE(6,*) 'with the name contained within quotes' WRITE(6,*) 'as example: ''A:OUTPUT.DTA'' ' WRITE(6,*) ' ' READ(5,*) NAME1 OUP = 3 OPEN(UNIT=OUP,FILE=NAME1,STATUS='NEW') ELSE OUP = 6 ENDIF WRITE(OUP,*) 'GRAGG EXTRAPOLATION METHOD' WRITE(OUP,6) 6 FORMAT(12X,'t(i)',12X,'w(i)',15X,'h',2X,'k') C STEP 2 TO=A WO=ALPHA H=HMAX C IFLAG IS USED TO EXIT THE LOOP IN STEP 4 IFLAG=1 C STEP 3 DO 110 I=1,7 DO 110 J=1,I C ( Q(I,J) = H(J)**2/H(I+1)**2 ) 110 Q(I,J)=(FLOAT(NK(I+1))/FLOAT(NK(J)))**2 C STEP 4 100 IF (IFLAG.NE.1) GOTO 200 C STEP 5 K=1 C WHEN DESIRED ACCURACY ACHIEVED, NFLAG IS SET TO 1 NFLAG=0 C STEP 6 300 IF (K.GT.8 .OR. NFLAG.NE.0) GOTO 400 C STEP 7 HK=H/NK(K) T=TO W2=WO C EULER FIRST STEP W3=W2+HK*F(T,W2) T=TO+HK C STEP 8 M=NK(K)-1 DO 20 J=1,M W1=W2 W2=W3 C MIDPOINT METHOD W3=W1+2*HK*F(T,W2) 20 T=TO+(J+1)*HK C STEP 9 C SMOOTHING TO COMPUTE Y(K,1) Y(K)=(W3+W2+HK*F(T,W3))/2 C STEP 10 C NOTE: Y(K-1)=Y(K-1,1), Y(K-2)=Y(K-1,2),..., Y(1)=Y(K-1,K-1) SINCE C ONLY PREVIOUS ROW OF TABLE REQUIRED IF(K.GE.2) THEN C STEP 11 J=K C SAVE Y(K-1,K-1) V=Y(1) C STEP 12 500 IF (J.LT.2) GOTO 600 C EXTRAPOLATION TO COMPUTE Y(J-1)=Y(K,K-J+2) C NOTE: Y(J-1) = (H(J-1)**2*Y(J)-H(K)**2*Y(J-1))/(H(J-1)**2-H(K)**2) Y(J-1)=Y(J)+(Y(J)-Y(J-1))/(Q(K-1,J-1)-1) J=J-1 GOTO 500 C STEP 13 600 IF(ABS(Y(1)-V).LE.TOL) NFLAG=1 C Y(1) ACCEPTED AS NEW W END IF C STEP 14 K=K+1 GOTO 300 C STEP 15 400 K=K-1 C STEP 16 IF(NFLAG.EQ.0) THEN C STEP 17 C NEW VALUE FOR W REJECTED, DECREASE H H=H/2 C STEP 18 IF(H.LT.HMIN) THEN WRITE(OUP,1) IFLAG=0 END IF ELSE C STEP 19 C NEW VALUE FOR W ACCEPTED WO=Y(1) TO=TO+H C STEP 20 WRITE(OUP,2) TO,WO,H,K IF(ABS(TO-B).LT.HMIN/2)THEN IFLAG=0 ELSE IF(TO+H.GT.B)THEN H=B-TO ELSE IF(K.LE.3.AND.H.LT.HMAX/2) THEN H=2*H ENDIF ENDIF ENDIF END IF GOTO 100 C STEP 21 200 CONTINUE 700 CLOSE(UNIT=5) CLOSE(UNIT=OUP) IF(OUP.NE.6) CLOSE(UNIT=6) STOP 1 FORMAT(1X,'FAILURE') 2 FORMAT(1X,E15.8,1X,E15.8,1X,E15.8,1X,I3) END