TWODQ

Computes a two-dimensional iterated integral.

Required Arguments

F — User-supplied FUNCTION to be integrated. The form is

F(X,Y), where

F(X,Y), where

X – First argument of F. (Input)

Y – Second argument of F. (Input)

F – The function value. (Output)

F must be declared EXTERNAL in the calling program.

A — Lower limit of outer integral. (Input)

B — Upper limit of outer integral. (Input)

G — User-supplied FUNCTION to evaluate the lower limits of the inner integral. The form is G(X), where

X – Only argument of G. (Input)

G – The function value. (Output)

G must be declared EXTERNAL in the calling program.

H — User-supplied FUNCTION to evaluate the upper limits of the inner integral. The form is H(X), where

X – Only argument of H. (Input)

H – The func`tion value. (Output)

H must be declared EXTERNAL in the calling program.

RESULT — Estimate of the integral from A to B of F. (Output)

Optional Arguments

ERRABS — Absolute accuracy desired. (Input)

Default: ERRABS = 1.e-3 for single precision and 1.d-8 for double precision.

Default: ERRABS = 1.e-3 for single precision and 1.d-8 for double precision.

ERRREL — Relative accuracy desired. (Input)

Default: ERRREL = 1.e-3 for single precision and 1.d-8 for double precision.

Default: ERRREL = 1.e-3 for single precision and 1.d-8 for double precision.

IRULE --- Choice of quadrature rule. (Input)

Default: IRULE = 2.

The Gauss-Kronrod rule is used with the following points:

Default: IRULE = 2.

The Gauss-Kronrod rule is used with the following points:

IRULE | Points |
---|---|

1 | 7-15 |

2 | 10-21 |

3 | 15-31 |

4 | 20-41 |

5 | 25-51 |

6 | 30-61 |

If the function has a peak singularity, use IRULE = 1. If the function is oscillatory, use IRULE = 6.

ERREST — Estimate of the absolute value of the error. (Output)

FORTRAN 90 Interface

Generic: CALL TWODQ (F, A, B, G, H, RESULT [, …])

Specific: The specific interface names are S_TWODQ and D_TWODQ.

FORTRAN 77 Interface

Single: CALL TWODQ (F, A, B, G, H, ERRABS, ERRREL, IRULE, RESULT, ERREST)

Double: The double precision name is DTWODQ.

Description

The routine TWODQ approximates the two-dimensional iterated integral

with the approximation returned in RESULT. An estimate of the error is returned in ERREST. The approximation is achieved by iterated calls to QDAG. Thus, this algorithm will share many of the characteristics of the routine QDAG. As in QDAG, several options are available. The absolute and relative error must be specified, and in addition, the Gauss-Kronrod pair must be specified (IRULE). The lower-numbered rules are used for less smooth integrands while the higher-order rules are more efficient for smooth (oscillatory) integrands.

Comments

1. Workspace may be explicitly provided, if desired, by use of T2ODQ/DT2ODQ. The reference is:

CALL T2ODQ (F, A, B, G, H, ERRABS, ERRREL, IRULE, RESULT, ERREST, MAXSUB, NEVAL, NSUBIN, ALIST, BLIST, RLIST, ELIST, IORD, WK, IWK)

The additional arguments are as follows:

MAXSUB — Number of subintervals allowed. (Input)

A value of 250 is used by TWODQ.

A value of 250 is used by TWODQ.

NEVAL — Number of evaluations of F. (Output)

NSUBIN — Number of subintervals generated in the outer integral. (Output)

ALIST — Array of length MAXSUB containing a list of the NSUBIN left endpoints for the outer integral. (Output)

BLIST — Array of length MAXSUB containing a list of the NSUBIN right endpoints for the outer integral. (Output)

RLIST — Array of length MAXSUB containing approximations to the NSUBIN integrals over the intervals defined by ALIST, BLIST, pertaining only to the outer integral. (Output)

ELIST — Array of length MAXSUB containing the error estimates of the NSUBIN values in RLIST. (Output)

IORD — Array of length MAXSUB. (Output)

Let K be NSUBIN if NSUBIN.LE.(MAXSUB/2 + 2), MAXSUB + 1 − NSUBIN otherwise. Then the first K locations contain pointers to the error estimates over the corresponding subintervals, such that ELIST(IORD(1)), …, ELIST(IORD(K)) form a decreasing sequence.

Let K be NSUBIN if NSUBIN.LE.(MAXSUB/2 + 2), MAXSUB + 1 − NSUBIN otherwise. Then the first K locations contain pointers to the error estimates over the corresponding subintervals, such that ELIST(IORD(1)), …, ELIST(IORD(K)) form a decreasing sequence.

WK — Work array of length 4 * MAXSUB, needed to evaluate the inner integral.

IWK — Work array of length MAXSUB, needed to evaluate the inner integral.

2. Informational errors

Type | Code | Description |
---|---|---|

4 | 1 | The maximum number of subintervals allowed has been reached. |

3 | 2 | Roundoff error, preventing the requested tolerance from being achieved, has been detected. |

3 | 3 | A degradation in precision has been detected. |

3. If EXACT is the exact value, TWODQ attempts to find RESULT such that ABS(EXACT − RESULT).LE.MAX(ERRABS, ERRREL * ABS(EXACT)). To specify only a relative error, set ERRABS to zero. Similarly, to specify only an absolute error, set ERRREL to zero.

Examples

Example 1

In this example, we approximate the integral

The value of the error estimate is machine dependent.

USE TWODQ_INT

USE UMACH_INT

IMPLICIT NONE

INTEGER IRULE, NOUT

REAL A, B, ERRABS, ERREST, ERRREL, F, G, H, RESULT

EXTERNAL F, G, H

! Get output unit number

CALL UMACH (2, NOUT)

! Set limits of integration

A = 0.0

B = 1.0

! Set error tolerances

ERRABS = 0.0

ERRREL = 0.01

! Parameter for oscillatory function

IRULE = 6

CALL TWODQ (F, A, B, G, H, RESULT, ERRABS, ERRREL, IRULE, errest)

! Print results

WRITE (NOUT,99999) RESULT, ERREST

99999 FORMAT (’ Result =’, F8.3, 13X, ’ Error estimate = ’, 1PE9.3)

END

!

REAL FUNCTION F (X, Y)

REAL X, Y

REAL COS

INTRINSIC COS

F = Y*COS(X+Y*Y)

RETURN

END

!

REAL FUNCTION G (X)

REAL X

G = 1.0

RETURN

END

!

REAL FUNCTION H (X)

REAL X

H = 3.0

RETURN

END

Output

Result = -0.514 Error estimate = 3.065E-06

Example 2

We modify the above example by assuming that the limits for the inner integral depend on x and, in particular, are g(x) = −2x and h(x) = 5x. The integral now becomes

The value of the error estimate is machine dependent.

USE TWODQ_INT

USE UMACH_INT

! Declare F, G, H

INTEGER IRULE, NOUT

REAL A, B, ERRABS, ERREST, ERRREL, F, G, H, RESULT

EXTERNAL F, G, H

!

CALL UMACH (2, NOUT)

! Set limits of integration

A = 0.0

B = 1.0

! Set error tolerances

ERRABS = 0.001

ERRREL = 0.0

! Parameter for oscillatory function

IRULE = 6

CALL TWODQ (F, A, B, G, H, RESULT, ERRABS, ERRREL, IRULE, ERREST)

! Print results

WRITE (NOUT,99999) RESULT, ERREST

99999 FORMAT (’ Computed =’, F8.3, 13X, ’ Error estimate = ’, 1PE9.3)

END

REAL FUNCTION F (X, Y)

REAL X, Y

!

REAL COS

INTRINSIC COS

!

F = Y*COS(X+Y*Y)

RETURN

END

REAL FUNCTION G (X)

REAL X

!

G = -2.0*X

RETURN

END

REAL FUNCTION H (X)

REAL X

!

H = 5.0*X

RETURN

END

Output

Computed = -0.083 Error estimate = 2.095E-06