AMI-5: An Adaptive
Multidimensional Integrator
AMI-5 is an adaptive multidimensional five-point integrator to approximate the iterated integral
with u not necessarily bounded.
Here, each lower-limit function fi and each upper-limit function gi
is a function of i variables, i = 0, ... , n - 1 (meaning, of
course, that f0 and g0 are constants).
You can download the zipped source with a couple of examples, or you can have fun
figuring out the simulation of recursive descent in a nonrecursive environment by looking at the source listing shown below.
<![CDATA[
C A M I - 5
C
C-----------------------------------------------------------------------
SUBROUTINE AMI5(N,A,B,F,G,U,BOUND,STEP,S,NS,V,IFLAG)
C-----------------------------------------------------------------------
C
C
C A D A P T I V E
C
C M U L T I D I M E N S I O N A L
C
C I N T E G R A T O R
C
C F I V E - P O I N T E X T R A P O L A T E D S I M P S O N
C
C
C-----------------------------------------------------------------------
C
C LUCIO TAVERNINI
C DIVISION OF MATHEMATICS, COMPUTER SCIENCE, AND STATISTICS
C UNIVERSITY OF TEXAS AT SAN ANTONIO
C SAN ANTONIO, TEXAS 78285
C
C AUGUST 15, 1983 -- REVISED JULY 12, 1988
C
C ADAPTED FROM THE ORIGINAL MDQUAD, WRITTEN FOR THE UNIVERSITY OF
C WISCONSIN COMPUTER CENTER LIBRARY IN 1968
C
C-----------------------------------------------------------------------
C
C --- ADAPTIVE APPROXIMATION TO AN N-FOLD ITERATED INTEGRAL ---
C
C THE ITERATED INTEGRALS ARE NUMBERED FROM 1 TO N STARTING
C WITH THE OUTERMOST INTEGRAL AND PROCEEDING TO THE INNERMOST.
C
C-----------------------------------------------------------------------
C
C
C INPUT PARAMETERS:
C
C
C - N -
C
C N IS THE DIMENSION OF THE INTEGRAL.
C
C DATATYPE: INTEGER.
C
C REQUIREMENT: 0 < N < 11.
C
C
C - A , B -
C
C A AND B ARE, RESPECTIVELY, THE LOWER AND UPPER LIMITS OF
C INTEGRATION FOR THE OUTERMOST INTEGRAL.
C
C DATATYPE: REAL*8.
C
C REQUIREMENT: A < B.
C
C
C - F , G -
C
C F AND G ARE ENTRY POINT NAMES OF USER-SUPPLIED FUNCTIONS
C WHICH DEFINE, RESPECTIVELY, THE LOWER AND UPPER LIMITS OF
C INTEGRATION. F MUST BE CODED BY THE USER TO BE EQUIVALENT
C TO THE FOLLOWING. (G MUST BE CODED SIMILARLY.)
C
C
C FUNCTION F(I,X)
C INTEGER I
C REAL*8 F,X(1)
C
C GO TO(1,2,...,N-1),I
C
C 1 F = VALUE OF LOWER LIMIT FOR LEVEL I+1
C AS A FUNCTION OF X(1)
C RETURN
C
C 2 F = VALUE OF LOWER LIMIT FOR LEVEL I+2
C AS A FUNCTION OF X(1) AND OF X(2)
C RETURN
C .
C .
C .
C
C N-1 F = VALUE OF LOWER LIMIT FOR LEVEL N
C AS A FUNCTION OF X(1),...,X(N-1)
C RETURN
C
C END
C
C
C REQUIREMENT: F MUST RETURN A VALUE NOT GREATER THAN THE
C VALUE RETURNED BY G OVER THE REGION OF
C INTEGRATION. IF THIS CONDITION FAILS TO BE
C SATISFIED, THE INTEGRAL TO WHICH THE LIMITS
C APPLY IS TAKEN TO BE ZERO. (THAT IS, THE
C LIMITS ARE TAKEN TO BE IDENTICAL.)
C
C
C - U -
C
C U IS THE ENTRY POINT NAME OF A USER-SUPPLIED REAL*8 FUNCTION
C WHICH EVALUATES THE INTEGRAND. U MUST BE CODED BY THE USER
C TO BE EQUIVALENT TO THE FOLLOWING.
C
C FUNCTION U(X)
C REAL*8 U,X(1)
C U = VALUE OF THE INTEGRAND AT X(1),...,X(N)
C RETURN
C END
C
C
C - BOUND -
C
C BOUND IS THE BOUND ON THE ESTIMATED GLOBAL DISCRETIZATION ERROR
C AMI USES SIMPSON'S RULE TO APPROXIMATE EACH ONE-
C DIMENSIONAL INTEGRAL. A NONUNIFORM GRID IS USED. THE
C SIMPSON APPROXIMATION COMPUTED ON EACH INTERVAL IS
C COMPARED TO THE RESULT OBTAINED BY APPLYING SIMPSON'S
C RULE TWICE, AFTER SPLITTING THE INTERVAL INTO TWO EQUAL
C HALVES. EXTRAPOLATION TO ZERO IS USED TO ESTIMATE THE
C LOCAL DISCRETIZATION ERROR.
C
C DATATYPE: REAL*8.
C
C REQUIREMENT: BOUND > 0.
C
C
C - STEP -
C
C STEP SPECIFIES THE SMALLEST INTEGRATION STEP WHICH CAN BE
C USED BY AMI. SHOULD THE INTEGRATION STEP ATTAIN A
C SMALLER VALUE DURING THE COMPUTATIONS, AMI ACCEPTS THE
C CURRENT ESTIMATE EVEN THOUGH THE ESTIMATED ERROR MAY
C EXCEED THE ERROR BOUND.
C
C DATATYPE: REAL*8.
C
C REQUIREMENT: STEP >= 0.
C
C
C - S , NS -
C
C S AND NS SPECIFY THE USER-SUPPLIED AMI WORKSPACE. S IS THE
C FIRST-ELEMENT NAME OF A FORTRAN ARRAY WHICH CAN BE USED BY
C AMI AS A SCRATCH WORKSPACE. NS SPECIFIES THE NUMBER
C OF ELEMENTS OF SIZE REAL*8 AVAILABLE TO AMI. AMI REVERTS
C TO A UNIFORM GRID IF INSUFFICIENT STORAGE IS AVAILABLE FOR
C THE USE OF ADAPTIVE INTEGRATION THROUGHOUT.
C
C DATATYPE: THE DATATYPE OF S IS NOT IMPORTANT,
C THE DATATYPE OF NS MUST BE INTEGER.
C
C REQUIREMENT: NS MUST INDICATE THE NUMBER OF ELEMENTS
C OF SIZE REAL*8 THAT ARE AVAILABLE. THE VALUE OF
C NS MUST BE AT LEAST 5*N. IF NS = 5*N,
C A UNIFORM GRID IS USED THROUGHOUT.
C
C
C-----------------------------------------------------------------------
C
C
C OUTPUT PARAMETERS:
C
C
C - V -
C
C V CONTAINS THE VALUE OF THE APPROXIMATE INTEGRAL.
C
C DATATYPE: REAL*8.
C
C
C - IFLAG -
C
C IFLAG RETURNS INFORMATION TO THE USER. THE MEANING OF
C IFLAG IS TO BE INTERPRETED AS FOLLOWS.
C
C IFLAG = 0 MEANS THAT THE COMPUTATION HAS BEEN COMPLETED
C NORMALLY. IF NOT, EACH DIGITS OF (THE BASE 10
C REPRESENTATION OF) IFLAG IS ASSIGNED THE VALUE 1 TO INDICATE
C THAT ONE OF THE CONDITIONS SHOWN IN THE TABLE GIVEN BELOW HAS
C OCCURRED. OTHERWISE, THE DIGIT HAS VALUE 0. THE DIGITS ARE
C NUMBERED FROM RIGHT TO LEFT STARTING WITH 1.
C
C
C DIGIT MEANING
C ----- -------
C
C 1 FUNCTION F HAS FAILED TO RETURN A VALUE LESS
C THAN THE VALUE RETURNED BY FUNCTION G (WHEN
C BOTH WERE EVALUATED AT THE SAME POINT) DURING
C THE COMPUTATIONS.
C
C 2 INSUFFICIENT STORAGE HAS BEEN ALLOCATED FOR
C THE USE OF ADAPTIVE INTEGRATION THROUGHOUT.
C A UNIFORM MESH HAS BEEN USED FOR PORTIONS
C OR FOR ALL OF THE COMPUTATIONS.
C
C 3 THE LOWER BOUND ON THE INTEGRATION STEP WAS
C ATTAINED DURING THE COMPUTATIONS.
C
C 4 A FATAL ERROR HAS OCCURRED. REFER TO THE
C PRINTER/TERMINAL OUTPUT FOR ADDITIONAL
C INFORMATION.
C
C DATATYPE: INTEGER.
C
C
C-----------------------------------------------------------------------
C
C
C COMMON:
C
C
C BEFORE CALLING AMI5 FOR THE FIRST TIME, THE USER MUST ASSIGN
C A FORTRAN CHANNEL FOR ANY ERROR MESSAGES THAT MAY BE GENERATED
C BY AMI5. THE COMMON BLOCK NAME IS OUTCHL. IT CONTAINS A
C SINGLE INTEGER VARIABLE, WHERE THE USER MUST STORE THE CHANNEL
C NUMBER:
C
C
C INTEGER OUTCHL
C COMMON/OUTCHL/OUTCHL
C
C
C-----------------------------------------------------------------------
C
C
C ----------
C PARAMETERS
C ----------
C
INTEGER IFLAG,N,NS
REAL*8 A,B,BOUND,STEP,V,U,F,G
C
C USER-ALLOCATED WORKSPACE
REAL*8 S(NS)
C
C ------
C LOCALS
C ------
C
INTEGER NC,NT
C
C FORTRAN CHANNEL FOR ERROR MESSAGES
INTEGER OUTCHL
COMMON/OUTCHL/OUTCHL
C
C ASSIGN WORKSPACE
NT = NS/5
NC = NT/N
C
C TEST PARAMETERS
C
IFLAG = 0
C
IF( N .LT. 11 .AND. N .GT. 0 ) GO TO 10
IFLAG = 1000
WRITE(OUTCHL,100)
WRITE(OUTCHL,101) N
10 CONTINUE
C
IF( B .GT. A ) GO TO 20
IF( IFLAG .EQ. 0 ) WRITE(OUTCHL,100)
IFLAG = 1000
WRITE(OUTCHL,102) A,B
20 CONTINUE
C
IF( BOUND .GT. 0 ) GO TO 30
IF( IFLAG .EQ. 0 ) WRITE(OUTCHL,100)
IFLAG = 1000
WRITE(OUTCHL,103) BOUND
30 CONTINUE
C
IF( STEP .GE. 0 ) GO TO 40
IF( FLAG .EQ. 0 ) WRITE(OUTCHL,100)
IFLAG = 1000
WRITE(OUTCHL,104) STEP
40 CONTINUE
C
IF( NC .GT. 0 ) GO TO 50
IF( IFLAG .EQ. 0 ) WRITE(OUTCHL,100)
IFLAG = 1000
WRITE(OUTCHL,105) NS
50 CONTINUE
C
60 IF( IFLAG .GT. 0 ) RETURN
C
CALL AMI5WK(U,F,G,N,A,B,BOUND,STEP,V,IFLAG,NC,
* S(1),S(NT+1),S(2*NT+1))
C ------
RETURN
C ------
C
100 FORMAT('0 * A M I 5 * INPUT PARAMETER ERROR(S):')
C
101 FORMAT(16X,'PARAMETER 1: N =',I10/
* 16X,'REQUIREMENT: 0 < N < 11.')
C
102 FORMAT(16X,'PARAMETER 2: A =',E15.6/
* 16X,'PARAMETER 3: B =',E15.6/
* 16X,'REQUIREMENT: A < B.')
C
103 FORMAT(16X,'PARAMETER 7: BOUND =',E15.6/
* 16X,'REQUIREMENT : BOUND > 0.')
C
104 FORMAT(16X,'PARAMETER 8: STEP =',E15.6/
* 16X,'REQUIREMENT: STEP >= 0.')
C
105 FORMAT(16X,'PARAMETER 10: NS =',I3/
* 16X,'REQUIREMENT: NS >= 5*N, WHERE N IS PARAMETER 1.')
C
END
C-----------------------------------------------------------------------
SUBROUTINE AMI5WK(U,F,G,N,AEP,BEP,EPSMAX,STEP,
* VALUE,IFLAG,NC,A,I,Q)
C-----------------------------------------------------------------------
C
C SUBSCRIPTING SCHEME USED FOR THE J-TH INTEGRATION INTERVAL
C AT RECURSION LEVEL L TO STORE VALUES OF THE INTEGRAND:
C
C
C
C
C Q(1,L,J) Q(2,L,J) Q(3,L,J)
C
C | | |
C | | |
C +-------------+-------------+-------------+-------------+
C | |
C | |
C
C P(1,L) P(2,L)
C
C
C ------------------
C VARIABLE WORKSPACE
C ------------------
REAL*8 A(N,NC),I(N,NC),Q(3,N,NC)
C
C ---------------
C FIXED WORKSPACE
C ---------------
INTEGER COUNTR(10),J(10),JEND(10),JTAB(10),
* K(10),M(10),MSTART(10),NSTEPS(10)
C
LOGICAL C001,C010,C100
C
REAL*8 AEP,BEP,EON,EPS(10),EPSMAX,U,F,G,H(10),HO3,HO6,
* I1,I2,LENGTH(10),P(2,10),
* STEP,SUM0(10),SUM1(10),SUM2(10),SUM3(10),V(11),VALUE,X(10)
C
C
C +--------------+
C | INITIALIZE |
C +--------------+
C
EON = EPSMAX/N
ASSIGN 500 TO J500
ASSIGN 1100 TO J1100
ASSIGN 2300 TO J2300
C001 = .FALSE.
C010 = .FALSE.
C100 = .FALSE.
A(1,1) = AEP
LENGTH(1) = BEP - AEP
L = 0
C
C +---------------------+
C | RECURSIVE DESCENT |
C +---------------------+
C
C --------------------------------------------------
C BEGIN: INITIALIZE TABLES FOR ADAPTIVE INTEGRATION.
C --------------------------------------------------
C
C INCREMENT LEVEL COUNTER.
100 L = L + 1
C
C SKIP EVALUATION OF F AND G FOR OUTER LIMITS OF INTEGRATION.
IF( L .EQ. 1 ) GO TO 200
A(L,1) = F(L-1,X)
LENGTH(L) = G(L-1,X) - A(L,1)
200 CONTINUE
C
C INITIALIZE APPROXIMATION AT THIS LEVEL TO ZERO
V(L) = 0
C
C TEST FOR INTEGRATION INTERVAL OF LENGTH ZERO
C001 = C001 .OR. LENGTH(L) .LT. 0
IF( LENGTH(L) .LE. 0 ) GO TO 5100
C
C H(L) IS THE DISTANCE BETWEEN GRID POINTS FOR LEVEL L.
H(L) = .5*LENGTH(L)
C
C EPS(L) IS USED TO BOUND THE ESTIMATED DISCRETIZATION
C ERROR AT LEVEL L.
EPS(L) = EON
C
C DECREASE ERROR BOUND TO SATISFY GLOBAL BOUND ON
C ASYMPTOTIC ERROR ESTIMATE.
LM1 = L - 1
IF( LM1 .LT. 1 ) GO TO 400
DO 300 INDEX = 1,LM1
300 EPS(L) = EPS(L)/LENGTH(INDEX)
400 CONTINUE
C
C INITIALIZE POINTERS FOR ADAPTIVE INTEGRATION WORKLIST.
K(L) = -1
J(L) = 1
MSTART(L) = NC
M(L) = NC
JEND(L) = 0
C
C -----------------------------------------
C COMPUTE POINTS FOR SIMPSON APPROXIMATION.
C -----------------------------------------
COUNTR(L) = 1
430 IF( COUNTR(L) .GT. 3 ) GO TO 600
X(L) = A(L,1) + ( COUNTR(L) - 1)*H(L)
C ---------------------------------------
C RECURSIVE DESCENT OR DIRECT EVALUATION?
IF( L .EQ. N) GO TO 450
JTAB(L) = J500
GO TO 100
450 V(L+1) = U(X)
C ---------------------------------------
C STORE POINT.
500 Q(COUNTR(L),L,1) = V(L+1)
COUNTR(L) = COUNTR(L) + 1
GO TO 430
600 CONTINUE
C
C --------------------------------
C EVALUATE SIMPSON APPROXIMATIONS.
C --------------------------------
I(L,1) = LENGTH(L)/6.*( Q(1,L,1) + 4.*Q(2,L,1) + Q(3,L,1) )
C
C -------------------------------------------------
C END: INITIALIZE TABLES FOR ADAPTIVE INTEGRATION.
C -------------------------------------------------
C
C +----------------------------------------+
C | VARIABLE STEP (ADAPTIVE) INTEGRATION |
C +----------------------------------------+
C
C -------------------------------------
C ARE WE THROUGH WITH CURRENT WORKLIST?
C -------------------------------------
800 IF( J(L) .EQ. JEND(L) ) GO TO 3300
C
C WE ARE NOT THROUGH: EVALUATE MIDPOINTS
C BETWEEN PREVIOUS POINTS.
COUNTR(L) = 1
1020 IF( COUNTR(L) .GT. 2 ) GO TO 1200
X(L) = A(L,J(L)) + ( COUNTR(L) - 0.5 )*H(L)
C ---------------------------------------
C RECURSIVE DESCENT OR DIRECT EVALUATION?
IF( L .EQ. N ) GO TO 1050
JTAB(L) = J1100
GO TO 100
1050 V(L+1) = U(X)
C ---------------------------------------
C STORE VALUE OF NEW MIDPOINT.
1100 P(COUNTR(L),L) = V(L+1)
COUNTR(L) = COUNTR(L) + 1
GO TO 1020
1200 CONTINUE
C
C -----------------------------------------------------------
C EVALUATE NEW SIMPSON APPROXIMATIONS OVER PREVIOUS INTERVAL.
C -----------------------------------------------------------
C
J1 = J(L)
HO6 = H(L)/6.
I1 = HO6*( Q(1,L,J1) + 4.*P(1,L) + Q(2,L,J1) )
I2 = HO6*( Q(2,L,J1) + 4.*P(2,L) + Q(3,L,J1) )
C
C ERR = ASYMPTOTIC ERROR ESTIMATE.
ERR = ( I1 + I2 - I(L,J1) )/15.
C
C TEST ASYMPTOTIC ESTIMATE AGAINST BOUND.
C100 = C100 .OR. H(L) .LT. STEP
IF( ABS(ERR) .LT. EPS(L) .OR. H(L) .LT. STEP ) GO TO 3200
C ------------------------------------
C ESTIMATED ERROR IS ABOVE UPPER BOUND
C ------------------------------------
C TEST FOR STORAGE OVERFLOW
IF( M(L) .NE. J1 ) GO TO 3100
C
C --------------------------------------
C BEGIN: RECOVER FROM STORAGE OVERFLOW.
C --------------------------------------
C
C CHANGE TO A UNIFORM INTEGRATION STEP.
C010 = .TRUE.
SUM1(L) = Q(1,L,J1) + Q(3,L,J1)
SUM2(L) = Q(2,L,J1) + P(1,L) + P(2,L)
C
C SAVE H(L)
Q(1,L,J1) = H(L)
H(L) = .5*H(L)
NSTEPS(L) = 4
SUM0(L) = I1 + I2
C
C REPEAT UNTIL ESTIMATED ERROR IS BELOW THE PRESCRIBED BOUND
2100 SUM3(L) = 0
C -----------------------------------
C BEGIN: LOOP TO SUM FUNCTION VALUES.
C -----------------------------------
COUNTR(L) = 1
2200 IF( COUNTR(L) .GT. NSTEPS(L) ) GO TO 2400
X(L) = A(L,J1) + ( COUNTR(L) - 0.5 )*H(L)
C ---------------------------------------
C RECURSIVE DESCENT OR DIRECT EVALUATION?
IF( L .EQ. N ) GO TO 2250
JTAB(L) = J2300
GO TO 100
2250 V(L+1) = U(X)
C ---------------------------------------
2300 SUM3(L) = SUM3(L) + V(L+1)
COUNTR(L) = COUNTR(L) + 1
GO TO 2200
2400 CONTINUE
C ---------------------------------
C END: LOOP TO SUM FUNCTION VALUES.
C ---------------------------------
C
C EVALUATE SIMPSON APPROXIMATION.
I2 = H(L)/6.*( SUM1(L) + 2.*SUM2(L) + 4.*SUM3(L) )
ERR = ( I2 - SUM0(L) )/15
C
C TEST ASYMPTOTIC ERROR ESTIMATE.
C100 = C100 .OR. H(L) .LT. STEP
IF(ABS( ERR) .LT. EPS(L) .OR. H(L) .LT. STEP ) GO TO 2500
SUM2(L) = SUM2(L) + SUM3(L)
NSTEPS(L) = 2*NSTEPS(L)
H(L) = .5*H(L)
SUM0(L) = I2
GO TO 2100
C
2500 V(L) = ERR + I2 + V(L)
C
C RESTORE H(L)
H(L) = Q(1,L,J(L))
J(L) = J(L) + K(L)
GO TO 800
C
C -----------------------------------
C END: RECOVER FROM STORAGE OVERFLOW.
C -----------------------------------
C
3100 CONTINUE
C
C MORE STORAGE IS AVAILABLE. SPLIT THE CURRENT
C INTERVAL INTO TWO SUBINTERVALS AND PLACE THEM
C IN THE NEXT WORKLIST
M1 = M(L)
M2 = M(L) + K(L)
A(L,M1) = A(L,J1)
A(L,M2) = A(L,J1) + H(L)
I(L,M1) = I1
I(L,M2) = I2
Q(1,L,M1) = Q(1,L,J1)
Q(2,L,M1) = P(1,L)
Q(3,L,M1) = Q(2,L,J1)
Q(1,L,M2) = Q(2,L,J1)
Q(2,L,M2) = P(2,L)
Q(3,L,M2) = Q(3,L,J1)
M(L) = M2 + K(L)
J(L) = J(L) + K(L)
GO TO 800
3200 CONTINUE
C
C ESTIMATED ERROR IS WITHIN BOUNDS. ACCEPT APPROXIMATION.
V(L) = ERR + I1 + I2 + V(L)
C UPDATE CURRENT WORKLIST POINTER FOR NEXT SUBINTERVAL.
J(L) = J(L) + K(L)
GO TO 800
3300 CONTINUE
C -----------------------------------------
C WE ARE THROUGH WITH THE CURRENT WORKLIST.
C -----------------------------------------
C
C --------------------------
C IS THE NEW WORKLIST EMPTY?
C --------------------------
IF( M(L) .EQ. MSTART(L) ) GO TO 5100
C
C NEW WORKLIST IS NOT EMPTY. TEST TO SEE IN WHICH DIRECTION
C POINTERS ARE TO BE MOVING.
IF(K(L).GT.0) GO TO 4100
C
C POINTERS ARE TO MOVE FROM LOW TO HIGH
3400 K(L) = 1
J(L) = M(L) + 1
JEND(L) = NC + 1
M(L) = 1
MSTART(L) = 1
EPS(L) = .5*EPS(L)
H(L) = .5*H(L)
GO TO 800
4100 CONTINUE
C
C POINTERS ARE TO MOVE FROM HIGH TO LOW
K(L) = -1
J(L) = M(L) - 1
JEND(L) = 0
M(L) = NC
MSTART(L) = NC
EPS(L) = .5*EPS(L)
H(L) = .5*H(L)
GO TO 800
5100 CONTINUE
C
C +---------------------------------+
C | RETURN FROM RECURSIVE DESCENT |
C +---------------------------------+
C
C BACK UP LEVEL COUNTER.
L = L - 1
C
C ARE WE THROUGH WITH LEVEL ONE?
IF( L .NE. 0 ) GO TO 5200
C
C ----------------------------
C WE ARE AT LEVEL ZERO: DONE.
C ----------------------------
VALUE = V(1)
IF( C001 ) IFLAG = IFLAG + 1
IF( C010 ) IFLAG = IFLAG + 10
IF( C100 ) IFLAG = IFLAG + 100
RETURN
5200 CONTINUE
C
C RESTORE RETURN ADDRESS.
JX = JTAB(L)
C RETURN.
GO TO JX, (500,1100,2300)
C
END
]]>