SUBROUTINE Q1DAX(F,A,B,EPS,R,E,NINT,RST,W,NMAX,FMIN,FMAX,KF,IFLAG) C***BEGIN PROLOGUE Q1DAX C***DATE WRITTEN 821118 C***REVISION DATE 821118 C***CATEGORY NO. D1B1A1 C***KEYWORDS ADAPTIVE QUADRATURE, AUTOMATIC QUADRATURE C***AUTHOR KAHANER, DAVID K., SCIENTIFIC COMPUING DIVISION, NBS. C***PURPOSE APPROXIMATES ONE DIMENSIONAL INTEGRALS OF USER DEFINED C FUNCTIONS, FLEXIBLE USAGE. C C***DESCRIPTION C C Q1DAX IS A FLEXIBLE SUBROUTINE FOR THE AUTOMATIC EVALUATION C OF DEFINITE INTEGRALS OF A USER DEFINED FUNCTION C OF ONE VARIABLE. C C FOR EASIER TO USE ROUTINES SEE Q1DA OR Q1DB. C C CAPABILITIES OF Q1DAX (IN ADDITION TO THOSE OF Q1DA, Q1DB) C INCLUDE: C ABILITY TO RESTART A CALCULATION TO GREATER C ACCURACY WITHOUT PENALTY... C ABILITY TO SPECIFY AN INITIAL PARTITION OF C THE INTEGRATION INTERVAL... C ABILITY TO INCREASE THE WORK SPACE TO HANDLE C MORE DIFFICULT PROBLEMS... C OUTPUT OF LARGEST/SMALLEST INTEGRAND VALUE FOR C APPLICATIONS SUCH AS SCALING GRAPHS... C C A R G U M E N T S I N T H E C A L L S E Q U E N C E C C F (INPUT) THE NAME OF YOUR INTEGRAND FUNCTION. C THIS NAME MUST APPEAR IN AN EXTERNAL STATEMENT C IN ANY PROGRAM WHICH CALLS Q1DAX. C YOU MUST WRITE F IN THE FORM C FUNCTION F(X) C F=(EVALUATE INTEGRAND AT THE POINT X) C RETURN C END C A C B (INPUT) ENDPOINTS OF INTEGRATION INTERVAL C EPS (INPUT) ACCURACY TO WHICH THE INTEGRAL IS TO BE CALCULATED. C Q1DAX WILL TRY TO ACHIEVE RELATIVE ACCURACY, C E.G. SET EPS=.01 FOR 2 DIGITS, .001 FOR 3, ETC. C R (OUTPUT) THE ESTIMATE OF THE INTEGRAL C E (OUTPUT) THE ESTIMATE OF THE ABSOLUTE ERROR IN R. C NINT (INPUT C OUTPUT) C AS AN INPUT QUANTITY, NINT MUST BE SET TO C THE NUMBER OF SUBINTERVALS IN THE INITIAL C PARTITION OF [A,B]. FOR MOST PROBLEMS C THIS IS JUST 1, THE INTERVAL [A,B] ITSELF. C NINT MUST BE LESS THAN NMAX, SEE BELOW. C NINT IS USEFUL IF YOU WOULD LIKE TO HELP C Q1DAX LOCATE A DIFFICULT SPOT ON [A,B]. C IN THIS REGARD NINT IS USED ALONG C WITH THE ARRAY W (SEE BELOW). IF YOU SET C NINT=1 IT IS NOT NECESSARY TO BE CONCERNED C WITH W, EXCEPT THAT IT MUST BE DIMENSIONED... C AS AN EXAMPLE OF MORE GENERAL APPLICATIONS, C IF [A,B]=[0,1] BUT THE INTEGRAND JUMPS AT 0.3, C IT WOULD BE WISE TO SET NINT=2 AND THEN SET C W(1,1)=0.0 (LEFT ENDPOINT) C W(2,1)=0.3 (SINGULAR POINT) C W(3,1)=1.0 (RIGHT ENDPOINT) C IF YOU SET NINT GREATER THAN 1, BE SURE TO C CHECK THAT YOU HAVE ALSO SET C W(1,1)=A AND W(NINT+1,1)=B C AS AN OUTPUT QUANTITY, NINT GIVES THE C NUMBER OF SUBINTERVALS IN THE FINAL C PARTITION OF [A,B]. C RST (INPUT) A LOGICAL VARIABLE (E.G. TRUE OR FALSE) C SET RST=.FALSE. FOR INITIAL CALL TO Q1DAX C SET RST=.TRUE. FOR A SUBSEQUENT CALL, C E.G. ONE FOR WHICH MORE ACCURACY IS C DESIRED (SMALLER EPS). A RESTART ONLY C MAKES SENSE IF THE PRECEDING CALL RETURNED C WITH A VALUE OF IFLAG (SEE BELOW) LESS THAN 3. C ON A RESTART YOU MAY NOT CHANGE THE VALUES OF ANY C OTHER ARGUMENTS IN THE CALL SEQUENCE, EXCEPT EPS. C W(NMAX,6) (INPUT & C SCRATCH) C W IS AN ARRAY USED BY Q1DAX. C YOU M U S T INCLUDE A DIMENSION STATEMENT IN C YOUR CALLING PROGRAM TO ALLOCATE THIS STORAGE. C THIS SHOULD BE OF THE FORM C DIMENSION W(NMAX,6) C WHERE NMAX IS AN INTEGER. AN ADEQUATE VALUE OF C NMAX IS 50. IF YOU SET NINT>1 YOU MUST ALSO C INITIALIZE W, SEE NINT ABOVE. C NMAX (INPUT) AN INTEGER EQUAL TO THE FIRST SUBSCRIPT IN THE C DIMENSION STATEMENT FOR THE ARRAY W. THIS IS C ALSO EQUAL TO THE MAXIMUM NUMBER OF SUBINTERVALS C PERMITTED IN THE INTERNAL PARTITION OF [A,B]. C A VALUE OF 50 IS AMPLE FOR MOST PROBLEMS. C FMIN C FMAX (OUTPUT) THE SMALLEST AND LARGEST VALUES OF THE INTEGRAND C WHICH OCCURRED DURING THE CALCULATION. THE C ACTUAL INTEGRAND RANGE ON [A,B] MAY, OF COURSE, C BE GREATER BUT PROBABLY NOT BY MORE THAN 10%. C KF (OUTPUT) THE ACTUAL NUMBER OF INTEGRAND EVALUATIONS USED C BY Q1DAX TO APPROXIMATE THIS INTEGRAL. KF C WILL ALWAYS BE AT LEAST 30. C IFLAG (OUTPUT) TERMINATION FLAG...POSSIBLE VALUES ARE C 0 NORMAL COMPLETION, E SATISFIES C EEPS*ABS(R) C 2 NORMAL COMPLETION, E SATISFIES C EEPS C 3 NORMAL COMPLETION BUT EPS WAS TOO SMALL TO C SATISFY ABSOLUTE OR RELATIVE ERROR REQUEST. C 4 ABORTED CALCULATION BECAUSE OF SERIOUS ROUNDING C ERROR. PROBABLY E AND R ARE CONSISTENT. C 5 ABORTED CALCULATION BECAUSE OF INSUFFICIENT STORAGE. C R AND E ARE CONSISTENT. PERHAPS INCREASING NMAX C WILL PRODUCE BETTER RESULTS. C 6 ABORTED CALCULATION BECAUSE OF SERIOUS DIFFICULTIES C MEETING YOUR ERROR REQUEST. C 7 ABORTED CALCULATION BECAUSE EITHER EPS, NINT OR NMAX C HAS BEEN SET TO AN ILLEGAL VALUE. C 8 ABORTED CALCULATION BECAUSE YOU SET NINT>1 BUT FORGOT C TO SET W(1,1)=A AND W(NINT+1,1)=B C C T Y P I C A L P R O B L E M S E T U P C C DIMENSION W(50,6) C LOGICAL RST C EXTERNAL F C A=0.0 C B=1.0 C W(1,1)=A C W(2,1)=.3 [SET INTERNAL PARTITION POINT AT .3] C W(3,1)=B C NINT=2 [INITIAL PARTITION HAS 2 INTERVALS] C RST=.FALSE. C EPS=.001 C NMAX=50 C C 1 CALL Q1DAX(F,A,B,EPS,R,E,NINT,RST,W,NMAX,FMIN,FMAX,KF,IFLAG) C C IF(EPS.EQ. .0001 .OR. IFLAG.GE.3)STOP C RST=.TRUE C EPS=.0001 [ASK FOR ANOTHER DIGIT] C GO TO 1 C END C FUNCTION F(X) C IF(X.LT. .3) C 1 THEN C F=X**(0.2)*ALOG(X) C ELSE C F=SIN(X) C ENDIF C RETURN C END C C C R E M A R K C C WHEN YOU USE Q1ADX WITH NINT=1, WE HAVE BUILT A SMALL C AMOUNT OF RANDOMIZATION INTO IT. REPEATED CALLS DURING C THE SAME RUN WILL PRODUCE DIFFERENT, BUT HOPEFULLY C CONSISTENT, RESULTS. C C E N D O F D O C U M E N T A T I O N C***REFERENCES (NONE) C***ROUTINES CALLED (R1MACH,UNI,GL15T) C***END PROLOGUE Q1DAX