\begin{page}{manpageXXd01}{NAG On-line Documentation: d01} \beginscroll \begin{verbatim} D01(3NAG) Foundation Library (12/10/92) D01(3NAG) D01 -- Quadrature Introduction -- D01 Chapter D01 Quadrature 1. Scope of the Chapter This chapter provides routines for the numerical evaluation of definite integrals in one or more dimensions and for evaluating weights and abscissae of integration rules. 2. Background to the Problems The routines in this chapter are designed to estimate: (a) the value of a one-dimensional definite integral of the form: b / |f(x)dx (1) / a where f(x) is defined by the user, either at a set of points (x ,f(x )), for i=1,2,...,n where a=x w f(x ). (5) / -- i i a i=1 The points x within the interval [a,b] are known as the i abscissae, and the w are known as the weights. i More generally, if the integrand has the form (2), the corresponding formula is b N / -- |w(x)g(x)dx~= > w g(x ). (6) / -- i i a i=1 If the integrand is known only at a fixed set of points, these points must be used as the abscissae, and the weighted sum is calculated using finite-difference methods. However, if the functional form of the integrand is known, so that its value at any abscissa is easily obtained, then a wide variety of quadrature rules are available, each characterised by its choice of abscissae and the corresponding weights. The appropriate rule to use will depend on the interval [a,b] - whether finite or otherwise - and on the form of any w(x) factor in the integrand. A suitable value of N depends on the general behaviour of f(x); or of g(x), if there is a w(x) factor present. Among possible rules, we mention particularly the Gaussian formulae, which employ a distribution of abscissae which is optimal for f(x) or g(x) of polynomial form. The choice of basic rules constitutes one of the principles on which methods for one-dimensional integrals may be classified. The other major basis of classification is the implementation strategy, of which some types are now presented. (a) Single rule evaluation procedures A fixed number of abscissae, N, is used. This number and the particular rule chosen uniquely determine the weights and abscissae. No estimate is made of the accuracy of the result. (b) Automatic procedures The number of abscissae, N, within [a,b] is gradually increased until consistency is achieved to within a level of accuracy (absolute or relative) requested by the user. There are essentially two ways of doing this; hybrid forms of these two methods are also possible: (i) whole interval procedures (non-adaptive) A series of rules using increasing values of N are successively applied over the whole interval [a,b]. It is clearly more economical if abscissae already used for a lower value of N can be used again as part of a higher-order formula. This principle is known as optimal extension. There is no overlap between the abscissae used in Gaussian formulae of different orders. However, the Kronrod formulae are designed to give an optimal (2N+1)-point formula by adding (N+1) points to an N-point Gauss formula. Further extensions have been developed by Patterson. (ii) adaptive procedures The interval [a,b] is repeatedly divided into a number of sub-intervals, and integration rules are applied separately to each sub-interval. Typically, the subdivision process will be carried further in the neighbourhood of a sharp peak in the integrand, than where the curve is smooth. Thus, the distribution of abscissae is adapted to the shape of the integrand. Subdivision raises the problem of what constitutes an acceptable accuracy in each sub-interval. The usual global acceptability criterion demands that the sum of the absolute values of the error estimates in the sub-intervals should meet the conditions required of the error over the whole interval. Automatic extrapolation over several levels of subdivision may eliminate the effects of some types of singularities. An ideal general-purpose method would be an automatic method which could be used for a wide variety of integrands, was efficient (i.e., required the use of as few abscissae as possible), and was reliable (i.e., always gave results within the requested accuracy). Complete reliability is unobtainable, and generally higher reliability is obtained at the expense of efficiency, and vice versa. It must therefore be emphasised that the automatic routines in this chapter cannot be assumed to be 100% reliable. In general, however, the reliability is very high. 2.2. Multi-dimensional Integrals A distinction must be made between cases of moderately low dimensionality (say, up to 4 or 5 dimensions), and those of higher dimensionality. Where the number of dimensions is limited, a one-dimensional method may be applied to each dimension, according to some suitable strategy, and high accuracy may be obtainable (using product rules). However, the number of integrand evaluations rises very rapidly with the number of dimensions, so that the accuracy obtainable with an acceptable amount of computational labour is limited; for example a product 9 of 3-point rules in 20 dimensions would require more than 10 integrand evaluations. Special techniques such as the Monte Carlo methods can be used to deal with high dimensions. (a) Products of one-dimensional rules Using a two-dimensional integral as an example, we have b b [ b ] 1 2 N [ 2 ] / / -- [ / ] | | f(x,y)dydx~= > w [ | f(x ,y)dy] (7) / / -- i[ / i ] a a i=1 [ a ] 1 2 [ 2 ] b b 1 2 N N / / -- -- | | f(x,y)dydx~= > > w v f(x ,y ) (8) / / -- -- i j i j a a i=1 j=1 1 2 where (w ,x ) and (v ,y ) are the weights and abscissae of i i i i the rules used in the respective dimensions. A different one-dimensional rule may be used for each dimension, as appropriate to the range and any weight function present, and a different strategy may be used, as appropriate to the integrand behaviour as a function of each independent variable. For a rule-evaluation strategy in all dimensions, the formula (8) is applied in a straightforward manner. For automatic strategies (i.e., attempting to attain a requested accuracy), there is a problem in deciding what accuracy must be requested in the inner integral(s). Reference to formula (7) shows that the presence of a limited but random error in the y-integration for different values of x can produce a 'jagged' function of x, which i may be difficult to integrate to the desired accuracy and for this reason products of automatic one-dimensional routines should be used with caution (see also Lyness [3]). (b) Monte Carlo methods These are based on estimating the mean value of the integrand sampled at points chosen from an appropriate statistical distribution function. Usually a variance reducing procedure is incorporated to combat the fundamentally slow rate of convergence of the rudimentary form of the technique. These methods can be effective by comparison with alternative methods when the integrand contains singularities or is erratic in some way, but they are of quite limited accuracy. (c) Number theoretic methods These are based on the work of Korobov and Conroy and operate by exploiting implicitly the properties of the Fourier expansion of the integrand. Special rules, constructed from so-called optimal coefficients, give a particularly uniform distribution of the points throughout n-dimensional space and from their number theoretic properties minimize the error on a prescribed class of integrals. The method can be combined with the Monte Carlo procedure. (d) Sag-Szekeres method By transformation this method seeks to induce properties into the integrand which make it accurately integrable by the trapezoidal rule. The transformation also allows effective control over the number of integrand evaluations. (e) Automatic adaptive procedures An automatic adaptive strategy in several dimensions normally involves division of the region into subregions, concentrating the divisions in those parts of the region where the integrand is worst behaved. It is difficult to arrange with any generality for variable limits in the inner integral(s). For this reason, some methods use a region where all the limits are constants; this is called a hyper-rectangle. Integrals over regions defined by variable or infinite limits may be handled by transformation to a hyper-rectangle. Integrals over regions so irregular that such a transformation is not feasible may be handled by surrounding the region by an appropriate hyper-rectangle and defining the integrand to be zero outside the desired region. Such a technique should always be followed by a Monte Carlo method for integration. The method used locally in each subregion produced by the adaptive subdivision process is usually one of three types: Monte Carlo, number theoretic or deterministic. Deterministic methods are usually the most rapidly convergent but are often expensive to use for high dimensionality and not as robust as the other techniques. 2.3. References Comprehensive reference: [1] Davis P J and Rabinowitz P (1975) Methods of Numerical Integration. Academic Press. Special topics: [2] Gladwell I (1986) Vectorisation of one dimensional quadrature codes. Techincal Report. TR7/86 NAG. [3] Lyness J N (1983) When not to use an automatic quadrature routine. SIAM Review. 25 63--87. [4] Piessens R, De Doncker-Kapenga E, Uberhuber C and Kahaner D (1983) QUADPACK, A Subroutine Package for Automatic Integration. Springer-Verlag. [5] Sobol I M (1974) The Monte Carlo Method. The University of Chicago Press. [6] Stroud A H (1971) Approximate Calculation of Multiple Integrals. Prentice-Hall. 3. Recommendations on Choice and Use of Routines The following three sub-sections consider in turn routines for: one-dimensional integrals over a finite interval, and over a semi-infinite or an infinite interval; and multi-dimensional integrals. Within each sub-section, routines are classified by the type of method, which ranges from simple rule evaluation to automatic adaptive algorithms. The recommendations apply particularly when the primary objective is simply to compute the value of one or more integrals, and in these cases the automatic adaptive routines are generally the most convenient and reliable, although also the most expensive in computing time. Note however that in some circumstances it may be counter- productive to use an automatic routine. If the results of the quadrature are to be used in turn as input to a further computation (e.g. an 'outer' quadrature or an optimization problem), then this further computation may be adversely affected by the 'jagged performance profile' of an automatic routine; a simple rule-evaluation routine may provide much better overall performance. For further guidance, the article Lyness [3] is recommended. 3.1. One-dimensional Integrals over a Finite Interval (a) Integrand defined as a set of points If f(x) is defined numerically at four or more points, then the Gill-Miller finite difference method (D01GAF) should be used. The interval of integration is taken to coincide with the range of x-values of the points supplied. It is in the nature of this problem that any routine may be unreliable. In order to check results independently and so as to provide an alternative technique the user may fit the integrand by Chebyshev series using E02ADF and then use routines E02AJF and E02AKF to evaluate its integral (which need not be restricted to the range of the integration points, as is the case for D01GAF). A further alternative is to fit a cubic spline to the data using E02BAF and then to evaluate its integral using E02BDF. (b) Integrand defined as a function If the functional form of f(x) is known, then one of the following approaches should be taken. They are arranged in the order from most specific to most general, hence the first applicable procedure in the list will be the most efficient. However, if the user does not wish to make any assumptions about the integrand, the most reliable routine to use will be D01AJF, although this will in general be less efficient for simple integrals. (i) Rule-evaluation routines If f(x) is known to be sufficiently well-behaved (more precisely, can be closely approximated by a polynomial of moderate degree), a Gaussian routine with a suitable number of abscissae may be used. D01BBF may be used if it is required to examine the weights and abscissae. In this case, the user should write the code for the evaluation of quadrature summation (6). (ii) Automatic adaptive routines Firstly, several routines are available for integrands of the form w(x)g(x) where g(x) is a ' smooth' function (i.e., has no singularities, sharp peaks or violent oscillations in the interval of integration) and w(x) is a weight function of one of the following forms: (alpha) (beta) k l if w(x)=(b-x) (x-a) (log(b-x)) (log(x-a)) where k,l=0 or 1, (alpha),(beta)>-1: use D01APF; if w(x)=1/(x-c): use D01AQF (this integral is called the Hilbert transform of g); if w(x)=cos((omega)x) or sin((omega)x): use D01ANF (this routine can also handle certain types of singularities in g(x)). Secondly, there are some routines for general f(x). If f(x) is known to be free of singularities, though it may be oscillatory, D01AKF may be used. The most powerful of the finite interval integration routine is D01AJF (which can cope with singularities of several types). It may be used if none of the more specific situations described above applies. D01AJF is very reliable, particularly where the integrand has singularities other than at an end-point, or has discontinuities or cusps, and is therefore recommended where the integrand is known to be badly- behaved, or where its nature is completely unknown. Most of the routines in this chapter require the user to supply a function or subroutine to evaluate the integrand at a single point. If f(x) has singularities of certain types, discontinuities or sharp peaks occurring at known points, the integral should be evaluated separately over each of the subranges or D01ALF may be used. 3.2. One-dimensional Integrals over a Semi-infinite or Infinite Interval (a) Integrand defined as a set of points If f(x) is defined numerically at four or more points, and the portion of the integral lying outside the range of the points supplied may be neglected, then the Gill-Miller finite difference method, D01GAF, should be used. (b) Integrand defined as a function (i) Rule evaluation routines If f(x) behaves approximately like a polynomial in x, apart from a weight function of the form -(beta)x e (beta)>0 (semi-infinite interval, lower limit finite); -(beta)x or e (beta)<0 (semi-infinite interval, upper limit finite); 2 -(beta)(x-(alpha)) or e (beta)>0 (infinite interval); or if f(x) behaves approximately like a -1 polynomial in (x+B) (semi-infinite range); then the Gaussian routines may be used. D01BBF may be used if it is required to examine the weights and abscissae. In this case, the user should write the code for the evaluation of quadrature summation (6). (ii) Automatic adaptive routines D01AMF may be used, except for integrands which decay slowly towards an infinite end-point, and oscillate in sign over the entire range. For this class, it may be possible to calculate the integral by integrating between the zeros and invoking some extrapolation process. D01ASF may be used for integrals involving weight functions of the form cos((omega)x) and sin((omega)x) over a semi-infinite interval (lower limit finite). The following alternative procedures are mentioned for completeness, though their use will rarely be necessary: (1) If the integrand decays rapidly towards an infinite end-point, a finite cut-off may be chosen, and the finite range methods applied. (2) If the only irregularities occur in the finite part (apart from a singularity at the finite limit, with which D01AMF can cope), the range may be divided, with D01AMF used on the infinite part. (3) A transformation to finite range may be employed, e.g. 1-t x= --- or x=-log t t e will transform (0,infty) to (1,0) while for infinite ranges we have +infty infty / / | f(x)dx= | [f(x)+f(-x)]dx. / / -infty 0 If the integrand behaves badly on (-infty,0) and well on (0,infty) or vice versa it is better to compute it as 0 infty / / | f(x)dx+ | f(x)dx. / / -infty 0 This saves computing unnecessary function values in the semi-infinite range where the function is well behaved. 3.3. Multi-dimensional Integrals A number of techniques are available in this area and the choice depends to a large extent on the dimension and the required accuracy. It can be advantageous to use more than one technique as a confirmation of accuracy particularly for high dimensional integrations. Many of the routines incorporate the transformation procedure REGION which allows general product regions to be easily dealt with in terms of conversion to the standard n-cube region. (a) Products of one-dimensional rules (suitable for up to about 5 dimensions) If f(x ,x ,...,x ) is known to be a sufficiently well- 1 2 n behaved function of each variable x , apart possibly from i weight functions of the types provided, a product of Gaussian rules may be used. These are provided by D01BBF. In this case, the user should write the code for the evaluation of quadrature summation (6). Rules for finite, semi-infinite and infinite ranges are included. The one-dimensional routines may also be used recursively. For example, the two-dimensional integral b b 1 2 / / I= | | f(x,y)dydx / / a a 1 2 may be expressed as b 1 / I= | F(x)dx, / a 1 where b 2 / F(x)= | f(x,y)dy. / a 2 The user segment to evaluate F(x) will call the integration routine for the y-integration, which will call another user segment for f(x,y) as a function of y (x being effectively a constant). Note that, as Fortran is not a recursive language, a different library integration routine must be used for each dimension. Apart from this restriction, the full range of one-dimensional routines are available, for finite/infinite intervals, constant/variable limits, rule evaluation/automatic strategies etc. (b) Automatic routines (D01GBF and D01FCF) Both routines are for integrals of the form b b b 1 2 n / / / | | ... | f(x ,x ,...,x )dx dx ...dx . / / / 1 2 n n n-1 1 a a a 1 2 n D01GBF is an adaptive Monte Carlo routine. This routine is usually slow and not recommended for high accuracy work. It is a robust routine that can often be used for low accuracy results with highly irregular integrands or when n is large. D01FCF is an adaptive deterministic routine. Convergence is fast for well-behaved integrands. Highly accurate results can often be obtained for n between 2 and 5, using significantly fewer integrand evaluations than would be required by D01GBF. The routine will usually work when the integrand is mildly singular and for n<=10 should be used before D01GBF. If it is known in advance that the integrand is highly irregular, it is best to compare results from at least two different routines. There are many problems for which one or both of the routines will require large amounts of computing time to obtain even moderately accurate results. The amount of computing time is controlled by the number of integrand evaluations allowed by the user, and users should set this parameter carefully, with reference to the time available and the accuracy desired. 3.4. Decision Trees (i) One-dimensional integrals over a finite interval. (If in doubt, follow the downward branch.) Please see figure in printed Reference Manual (ii) One-dimensional integrals over a semi-infinite or infinite interval. (If in doubt, follow the downward branch.) Please see figure in printed Reference Manual D01 -- Quadrature Contents -- D01 Chapter D01 Quadrature D01AJF 1-D quadrature, adaptive, finite interval, strategy due to Piessens and de Doncker, allowing for badly-behaved integrands D01AKF 1-D quadrature, adaptive, finite interval, method suitable for oscillating functions D01ALF 1-D quadrature, adaptive, finite interval, allowing for singularities at user-specified break-points D01AMF 1-D quadrature, adaptive, infinite or semi-infinite interval D01ANF 1-D quadrature, adaptive, finite interval, weight function cos((omega)x) or sin((omega)x) D01APF 1-D quadrature, adaptive, finite interval, weight function with end-point singularities of algebraico- logarithmic type D01AQF 1-D quadrature, adaptive, finite interval, weight function 1/(x-c), Cauchy principal value (Hilbert transform) D01ASF 1-D quadrature, adaptive, semi-infinite interval, weight function cos((omega)x) or sin((omega)x) D01BBF Weights and abscissae for Gaussian quadrature rules D01FCF Multi-dimensional adaptive quadrature over hyper- rectangle D01GAF 1-D quadrature, integration of function defined by data values, Gill-Miller method D01GBF Multi-dimensional quadrature over hyper-rectangle, Monte Carlo method \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXd01ajf}{NAG On-line Documentation: d01ajf} \beginscroll \begin{verbatim} D01AJF(3NAG) Foundation Library (12/10/92) D01AJF(3NAG) D01 -- Quadrature D01AJF D01AJF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose D01AJF is a general-purpose integrator which calculates an approximation to the integral of a function f(x) over a finite interval [a,b]: b / I= |f(x)dx. / a 2. Specification SUBROUTINE D01AJF (F, A, B, EPSABS, EPSREL, RESULT, 1 ABSERR, W, LW, IW, LIW, IFAIL) INTEGER LW, IW(LIW), LIW, IFAIL DOUBLE PRECISION F, A, B, EPSABS, EPSREL, RESULT, ABSERR, W 1 (LW) EXTERNAL F 3. Description D01AJF is based upon the QUADPACK routine QAGS (Piessens et al [3]). It is an adaptive routine, using the Gauss 10-point and Kronrod 21-point rules. The algorithm, described by de Doncker [1], incorporates a global acceptance criterion (as defined by Malcolm and Simpson [2]) together with the (epsilon)-algorithm (Wynn [4]) to perform extrapolation. The local error estimation is described by Piessens et al [3]. The routine is suitable as a general purpose integrator, and can be used when the integrand has singularities, especially when these are of algebraic or logarithmic type. D01AJF requires the user to supply a function to evaluate the integrand at a single point. The routine D01ATF(*) uses an identical algorithm but requires the user to supply a subroutine to evaluate the integrand at an array of points. Therefore D01ATF(*) will be more efficient if the evaluation can be performed in vector mode on a vector- processing machine. 4. References [1] De Doncker E (1978) An Adaptive Extrapolation Algorithm for Automatic Integration. Signum Newsletter. 13 (2) 12--18. [2] Malcolm M A and Simpson R B (1976) Local Versus Global Strategies for Adaptive Quadrature. ACM Trans. Math. Softw. 1 129--146. [3] Piessens R, De Doncker-Kapenga E, Uberhuber C and Kahaner D (1983) QUADPACK, A Subroutine Package for Automatic Integration. Springer-Verlag. [4] Wynn P (1956) On a Device for Computing the e (S ) m n Transformation. Math. Tables Aids Comput. 10 91--96. 5. Parameters 1: F -- DOUBLE PRECISION FUNCTION, supplied by the user. External Procedure F must return the value of the integrand f at a given point. Its specification is: DOUBLE PRECISION FUNCTION F (X) DOUBLE PRECISION X 1: X -- DOUBLE PRECISION Input On entry: the point at which the integrand f must be evaluated. F must be declared as EXTERNAL in the (sub)program from which D01AJF is called. Parameters denoted as Input must not be changed by this procedure. 2: A -- DOUBLE PRECISION Input On entry: the lower limit of integration, a. 3: B -- DOUBLE PRECISION Input On entry: the upper limit of integration, b. It is not necessary that a= 4. 10: IW(LIW) -- INTEGER array Output On exit: IW(1) contains the actual number of sub-intervals used. The rest of the array is used as workspace. 11: LIW -- INTEGER Input On entry: the dimension of the array IW as declared in the (sub)program from which D01AJF is called. The number of sub-intervals into which the interval of integration may be divided cannot exceed LIW. Suggested value: LIW = LW/4. Constraint: LIW >= 1. 12: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. Users who are unfamiliar with this parameter should refer to the Essential Introduction for details. On exit: IFAIL = 0 unless the routine detects an error or gives a warning (see Section 6). For this routine, because the values of output parameters may be useful even if IFAIL /=0 on exit, users are recommended to set IFAIL to -1 before entry. It is then essential to test the value of IFAIL on exit. 6. Error Indicators and Warnings Errors or warnings specified by the routine: If on entry IFAIL = 0 or -1, explanatory error messages are output on the current error message unit (as defined by X04AAF). IFAIL= 1 The maximum number of subdivisions allowed with the given workspace has been reached without the accuracy requirements being achieved. Look at the integrand in order to determine the integration difficulties. If the position of a local difficulty within the interval can be determined (e.g. a singularity of the integrand or its derivative, a peak, a discontinuity, etc) you will probably gain from splitting up the interval at this point and calling the integrator on the subranges. If necessary, another integrator, which is designed for handling the type of difficulty involved, must be used. Alternatively, consider relaxing the accuracy requirements specified by EPSABS and EPSREL, or increasing the amount of workspace. IFAIL= 2 Round-off error prevents the requested tolerance from being achieved. The error may be under-estimated. Consider requesting less accuracy. IFAIL= 3 Extremely bad local integrand behaviour causes a very strong subdivision around one (or more) points of the interval. IFAIL= 4 The requested tolerance cannot be achieved, because the extrapolation does not increase the accuracy satisfactorily; the returned result is the best which can be obtained. The same advice applies as in the case of IFAIL = 1. IFAIL= 5 The integral is probably divergent, or slowly convergent. Please note that divergence can occur with any non-zero value of IFAIL. IFAIL= 6 On entry LW < 4, or LIW < 1. 7. Accuracy The routine cannot guarantee, but in practice usually achieves, the following accuracy: |I-RESULT|<=tol, where tol=max{|EPSABS|,|EPSREL|*|I|}, and EPSABS and EPSREL are user-specified absolute and relative error tolerance. Moreover it returns the quantity ABSERR which, in normal circumstances, satisfies |I-RESULT|<=ABSERR<=tol. 8. Further Comments The time taken by the routine depends on the integrand and the accuracy required. If IFAIL /= 0 on exit, then the user may wish to examine the contents of the array W, which contains the end-points of the sub-intervals used by D01AJF along with the integral contributions and error estimates over the sub-intervals. Specifically, for i=1,2,...,n, let r denote the approximation to i the value of the integral over the sub-interval [a ,b ] in the i i partition of [a,b] and e be the corresponding absolute error i estimate. b i n / -- Then, | f(x)dx~=r and RESULT= > r , unless D01AJF terminates / i -- i a i=1 i while testing for divergence of the integral (see Piessens et al [3], Section 3.4.3). In this case, RESULT (and ABSERR) are taken to be the values returned from the extrapolation process. The value of n is returned in IW(1), and the values a , b , e and r i i i i are stored consecutively in the array W, that is: a =W(i), i b =W(n+i), i e =W(2n+i) and i r =W(3n+i). i 9. Example To compute 2(pi) / x sin(30x) | ----------------dx. / 0 /( ( x )2) / (1-(-----) ) \/ ( (2(pi)) ) The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXd01akf}{NAG On-line Documentation: d01akf} \beginscroll \begin{verbatim} D01AKF(3NAG) Foundation Library (12/10/92) D01AKF(3NAG) D01 -- Quadrature D01AKF D01AKF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose D01AKF is an adaptive integrator, especially suited to oscillating, non-singular integrands, which calculates an approximation to the integral of a function f(x) over a finite interval [a,b]: b / I= |f(x)dx. / a 2. Specification SUBROUTINE D01AKF (F, A, B, EPSABS, EPSREL, RESULT, 1 ABSERR, W, LW, IW, LIW, IFAIL) INTEGER LW, IW(LIW), LIW, IFAIL DOUBLE PRECISION F, A, B, EPSABS, EPSREL, RESULT, ABSERR, W 1 (LW) EXTERNAL F 3. Description D01AKF is based upon the QUADPACK routine QAG (Piessens et al [3] ). It is an adaptive routine, using the Gauss 30-point and Kronrod 61-point rules. A 'global' acceptance criterion (as defined by Malcolm and Simpson [1]) is used. The local error estimation is described in Piessens et al [3]. Because this routine is based on integration rules of high order, it is especially suitable for non-singular oscillating integrands. D01AKF requires the user to supply a function to evaluate the integrand at a single point. The routine D01AUF(*) uses an identical algorithm but requires the user to supply a subroutine to evaluate the integrand at an array of points. Therefore D01AUF(*) will be more efficient if the evaluation can be performed in vector mode on a vector- processing machine. D01AUF(*) also has an additional parameter KEY which allows the user to select from six different Gauss-Kronrod rules. 4. References [1] Malcolm M A and Simpson R B (1976) Local Versus Global Strategies for Adaptive Quadrature. ACM Trans. Math. Softw. 1 129--146. [2] Piessens R (1973) An Algorithm for Automatic Integration. Angewandte Informatik. 15 399--401. [3] Piessens R, De Doncker-Kapenga E, Uberhuber C and Kahaner D (1983) QUADPACK, A Subroutine Package for Automatic Integration. Springer-Verlag. 5. Parameters 1: F -- DOUBLE PRECISION FUNCTION, supplied by the user. External Procedure F must return the value of the integrand f at a given point. Its specification is: DOUBLE PRECISION FUNCTION F (X) DOUBLE PRECISION X 1: X -- DOUBLE PRECISION Input On entry: the point at which the integrand f must be evaluated. F must be declared as EXTERNAL in the (sub)program from which D01AKF is called. Parameters denoted as Input must not be changed by this procedure. 2: A -- DOUBLE PRECISION Input On entry: the lower limit of integration, a. 3: B -- DOUBLE PRECISION Input On entry: the upper limit of integration, b. It is not necessary that a= 4. See IW below. 10: IW(LIW) -- INTEGER array Output On exit: IW(1) contains the actual number of sub-intervals used. The rest of the array is used as workspace. 11: LIW -- INTEGER Input On entry: the dimension of the array IW as declared in the (sub)program from which D01AKF is called. The number of sub-intervals into which the interval of integration may be divided cannot exceed LIW. Suggested value: LIW = LW/4. Constraint: LIW >= 1. 12: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. Users who are unfamiliar with this parameter should refer to the Essential Introduction for details. On exit: IFAIL = 0 unless the routine detects an error or gives a warning (see Section 6). For this routine, because the values of output parameters may be useful even if IFAIL /=0 on exit, users are recommended to set IFAIL to -1 before entry. It is then essential to test the value of IFAIL on exit. 6. Error Indicators and Warnings Errors or warnings specified by the routine: If on entry IFAIL = 0 or -1, explanatory error messages are output on the current error message unit (as defined by X04AAF). IFAIL= 1 The maximum number of subdivisions allowed with the given workspace has been reached without the accuracy requirements being achieved. Look at the integrand in order to determine the integration difficulties. Probably another integrator which is designed for handling the type of difficulty involved must be used. Alternatively, consider relaxing the accuracy requirements specified by EPSABS and EPSREL, or increasing the amount of workspace. IFAIL= 2 Round-off error prevents the requested tolerance from being achieved. Consider requesting less accuracy. IFAIL= 3 Extremely bad local integrand behaviour causes a very strong subdivision around one (or more) points of the interval. The same advice applies as in the case of IFAIL = 1. IFAIL= 4 On entry LW < 4, or LIW < 1. 7. Accuracy The routine cannot guarantee, but in practice usually achieves, the following accuracy: |I-RESULT|<=tol, where tol=max{|EPSABS|,|EPSREL|*|I|}, and EPSABS and EPSREL are user-specified absolute and relative error tolerances. Moreover it returns the quantity ABSERR which, in normal circumstances satisfies |I-RESULT|<=ABSERR<=tol. 8. Further Comments The time taken by the routine depends on the integrand and the accuracy required. If IFAIL /= 0 on exit, then the user may wish to examine the contents of the array W, which contains the end-points of the sub-intervals used by D01AKF along with the integral contributions and error estimates over these sub-intervals. Specifically, for i=1,2,...,n, let r denote the approximation to i the value of the integral over the sub-interval [a ,b ] in the i i partition of [a,b] and e be the corresponding absolute error i b i n / -- estimate. Then, | f(x)dx~=r and RESULT= > r . The value of n / i -- i a i=1 i is returned in IW(1), and the values a , b , e and r are stored i i i i consecutively in the array W, that is: a =W(i), i b =W(n+i), i e =W(2n+i) and i r =W(3n+i). i 9. Example To compute 2(pi) / | x sin(30x) cosx dx. / 0 The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXd01alf}{NAG On-line Documentation: d01alf} \beginscroll \begin{verbatim} D01ALF(3NAG) Foundation Library (12/10/92) D01ALF(3NAG) D01 -- Quadrature D01ALF D01ALF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose D01ALF is a general purpose integrator which calculates an approximation to the integral of a function f(x) over a finite interval [a,b]: b / I= |f(x)dx / a where the integrand may have local singular behaviour at a finite number of points within the integration interval. 2. Specification SUBROUTINE D01ALF (F, A, B, NPTS, POINTS, EPSABS, EPSREL, 1 RESULT, ABSERR, W, LW, IW, LIW, IFAIL) INTEGER NPTS, LW, IW(LIW), LIW, IFAIL DOUBLE PRECISION F, A, B, POINTS(*), EPSABS, EPSREL, 1 RESULT, ABSERR, W(LW) EXTERNAL F 3. Description D01ALF is based upon the QUADPACK routine QAGP (Piessens et al [3]). It is very similar to D01AJF, but allows the user to supply difficult. It is an adaptive routine, using the Gauss 10-point and Kronrod 21-point rules. The algorithm described by de Doncker [1], incorporates a global acceptance criterion (as defined by Malcolm and Simpson [2]) together with the (epsilon)-algorithm (Wynn [4]) to perform extrapolation. The user-supplied 'break- points'always occur as the end-points of some sub-interval during the adaptive process. The local error estimation is described by Piessens et al [3]. 4. References [1] De Doncker E (1978) An Adaptive Extrapolation Algorithm for Automatic Integration. Signum Newsletter. 13 (2) 12--18. [2] Malcolm M A and Simpson R B (1976) Local Versus Global Strategies for Adaptive Quadrature. ACM Trans. Math. Softw. 1 129--146. [3] Piessens R, De Doncker-Kapenga E, Uberhuber C and Kahaner D (1983) QUADPACK, A Subroutine Package for Automatic Integration. Springer-Verlag. [4] Wynn P (1956) On a Device for Computing the e (S ) m n Transformation. Math. Tables Aids Comput. 10 91--96. 5. Parameters 1: F -- DOUBLE PRECISION FUNCTION, supplied by the user. External Procedure F must return the value of the integrand f at a given point. Its specification is: DOUBLE PRECISION FUNCTION F (X) DOUBLE PRECISION X 1: X -- DOUBLE PRECISION Input On entry: the point at which the integrand f must be evaluated. F must be declared as EXTERNAL in the (sub)program from which D01ALF is called. Parameters denoted as Input must not be changed by this procedure. 2: A -- DOUBLE PRECISION Input On entry: the lower limit of integration, a. 3: B -- DOUBLE PRECISION Input On entry: the upper limit of integration, b. It is not necessary that a= 0. 5: POINTS(NPTS) -- DOUBLE PRECISION array Input On entry: the user-specified break-points. Constraint: the break-points must all lie within the interval of integration (but may be supplied in any order). 6: EPSABS -- DOUBLE PRECISION Input On entry: the absolute accuracy required. If EPSABS is negative, the absolute value is used. See Section 7. 7: EPSREL -- DOUBLE PRECISION Input On entry: the relative accuracy required. If EPSREL is negative, the absolute value is used. See Section 7. 8: RESULT -- DOUBLE PRECISION Input On entry: the approximation to the integral I. 9: ABSERR -- DOUBLE PRECISION Output On exit: an estimate of the modulus of the absolute error, which should be an upper bound for |I-RESULT|. 10: W(LW) -- DOUBLE PRECISION array Output On exit: details of the computation, as described in Section 8. 11: LW -- INTEGER Input On entry: the dimension of the array W as declared in the (sub)program from which D01ALF is called. The value of LW (together with that of LIW below) imposes a bound on the number of sub-intervals into which the interval of integration may be divided by the routine. The number of sub-intervals cannot exceed (LW-2*NPTS-4)/4. The more difficult the integrand, the larger LW should be. Suggested value: a value in the range 800 to 2000 is adequate for most problems. Constraint: LW>=2*NPTS+8. 12: IW(LIW) -- INTEGER array Output On exit: IW(1) contains the actual number of sub-intervals used. The rest of the array is used as workspace. 13: LIW -- INTEGER Input On entry: the dimension of the array IW as declared in the (sub)program from which D01ALF is called. The number of sub-intervals into which the interval of integration may be divided cannot exceed (LIW-NPTS-2)/2. Suggested value: LIW = LW/2. Constraint: LIW >= NPTS + 4. 14: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. Users who are unfamiliar with this parameter should refer to the Essential Introduction for details. On exit: IFAIL = 0 unless the routine detects an error or gives a warning (see Section 6). For this routine, because the values of output parameters may be useful even if IFAIL /=0 on exit, users are recommended to set IFAIL to -1 before entry. It is then essential to test the value of IFAIL on exit. 6. Error Indicators and Warnings Errors or warnings specified by the routine: If on entry IFAIL = 0 or -1, explanatory error messages are output on the current error message unit (as defined by X04AAF). IFAIL= 1 The maximum number of subdivisions allowed with the given workspace has been reached, without the accuracy requirements being achieved. Look at the integrand in order to determine the integration difficulties. If the position of a local difficulty within the interval can be determined (e.g. a singularity of the integrand or its derivative, a peak, a discontinuity, etc) it should be supplied to the routine as an element of the vector POINTS. If necessary, another integrator should be used, which is designed for handling the type of difficulty involved. Alternatively, consider relaxing the accuracy requirements specified by EPSABS and EPSREL, or increasing the amount of workspace. IFAIL= 2 Round-off error prevents the requested tolerance from being achieved. The error may be under-estimated. Consider requesting less accuracy. IFAIL= 3 Extremely bad local integrand behaviour causes a very strong subdivision around one (or more) points of the interval. The same advice applies as in the case of IFAIL = 1. IFAIL= 4 The requested tolerance cannot be achieved, because the extrapolation does not increase the accuracy satisfactorily; the result returned is the best which can be obtained. The same advice applies as in the case IFAIL = 1. IFAIL= 5 The integral is probably divergent, or slowly convergent. Please note that divergence can also occur with any other non-zero value of IFAIL. IFAIL= 6 The input is invalid: break-points are specified outside the integration range, NPTS > LIMIT or NPTS < 0. RESULT and ABSERR are set to zero. IFAIL= 7 On entry LW<2*NPTS+8, or LIW r unless D01ALF / i -- i a i=1 i terminates while testing for divergence of the integral (see Piessens et al [3] Section 3.4.3). In this case, RESULT (and ABSERR) are taken to be the values returned from the extrapolation process. The value of n is returned in IW(1), and the values a , b , e and r are stored consecutively in the i i i i array W, that is: a =W(i), i b =W(n+i), i e =W(2n+i) and i r =W(3n+i). i 9. Example To compute 1 / 1 | ---------dx. / 0 \/|x-1/7| A break-point is specified at x=1/7, at which point the integrand is infinite. (For definiteness the function FST returns the value 0.0 at this point.) The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXd01amf}{NAG On-line Documentation: d01amf} \beginscroll \begin{verbatim} D01AMF(3NAG) Foundation Library (12/10/92) D01AMF(3NAG) D01 -- Quadrature D01AMF D01AMF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose D01AMF calculates an approximation to the integral of a function f(x) over an infinite or semi-infinite interval [a,b]: b / I= |f(x)dx / a 2. Specification SUBROUTINE D01AMF (F, BOUND, INF, EPSABS, EPSREL, RESULT, 1 ABSERR, W, LW, IW, LIW, IFAIL) INTEGER INF, LW, IW(LIW), LIW, IFAIL DOUBLE PRECISION F, BOUND, EPSABS, EPSREL, RESULT, ABSERR, 1 W(LW) EXTERNAL F 3. Description D01AMF is based on the QUADPACK routine QAGI (Piessens et al [3]) [0,1] using one of the identities: a 1 / / ( 1-t) 1 | f(x)dx= |f(a- ---) --dt / / ( t ) 2 -infty 0 t infty 1 / / ( 1-t) 1 | f(x)dx= |f(a+ ---) --dt / / ( t ) 2 a 0 t infty infty 1 / / /[ ( 1-t) ( -1+t)] 1 | f(x)dx= | (f(x)+f(-x))dx= |[f( ---)+f( ----)] --dt / / /[ ( t ) ( t )] 2 -infty 0 0 t where a represents a finite integration limit. An adaptive procedure, based on the Gauss seven-point and Kronrod 15-point rules, is then employed on the transformed integral. The algorithm, described by de Doncker [1], incorporates a global acceptance criterion (as defined by Malcolm and Simpson [2]) together with the (epsilon)-algorithm (Wynn [4]) to perform extrapolation. The local error estimation is described by Piessens et al [3]. 4. References [1] De Doncker E (1978) An Adaptive Extrapolation Algorithm for Automatic Integration. Signum Newsletter. 13 (2) 12--18. [2] Malcolm M A and Simpson R B (1976) Local Versus Global Strategies for Adaptive Quadrature. ACM Trans. Math. Softw. 1 129--146. [3] Piessens R, De Doncker-Kapenga E, Uberhuber C and Kahaner D (1983) QUADPACK, A Subroutine Package for Automatic Integration. Springer-Verlag. [4] Wynn P (1956) On a Device for Computing the e (S ) m n Transformation. Math. Tables Aids Comput. 10 91--96. 5. Parameters 1: F -- DOUBLE PRECISION FUNCTION, supplied by the user. External Procedure On entry: the point at which the integrand f must be evaluated. Its specification is: DOUBLE PRECISION FUNCTION F (X) DOUBLE PRECISION X 1: X -- DOUBLE PRECISION Input On entry: the point at which the integrand f must be evaluated. F must be declared as EXTERNAL in the (sub)program from which D01AMF is called. Parameters denoted as Input must not be changed by this procedure. 2: BOUND -- DOUBLE PRECISION Input On entry: the finite limit of the integration range (if present). BOUND is not used if the interval is doubly infinite. 3: INF -- INTEGER Input On entry: indicates the kind of integration range: if INF =1, the range is [BOUND, +infty) if INF =-1, the range is (-infty, BOUND] if INF =+2, the range is (-infty, +infty). Constraint: INF =-1, 1 or 2. 4: EPSABS -- DOUBLE PRECISION Input On entry: the absolute accuracy required. If EPSABS is negative, the absolute value is used. See Section 7. 5: EPSREL -- DOUBLE PRECISION Input On entry: the relative accuracy required. If EPSREL is negative, the absolute value is used. See Section 7. 6: RESULT -- DOUBLE PRECISION Output On exit: the approximation to the integral I. 7: ABSERR -- DOUBLE PRECISION Output On exit: an estimate of the modulus of the absolute error, which should be an upper bound for |I-RESULT|. 8: W(LW) -- DOUBLE PRECISION array Output On exit: details of the computation, as described in Section 8. 9: LW -- INTEGER Input On entry: the dimension of the array W as declared in the (sub)program from which D01AMF is called. The value of LW (together with that of LIW below) imposes a bound on the number of sub-intervals into which the interval of integration may be divided by the routine. The number of sub-intervals cannot exceed LW/4. The more difficult the integrand, the larger LW should be. Suggested value: a value in the range 800 to 2000 is adequate for most problems. Constraint: LW >= 4. 10: IW(LIW) -- INTEGER array Output On exit: IW(1) contains the actual number of sub-intervals used. The rest of the array is used as workspace. 11: LIW -- INTEGER Input On entry: the dimension of the array IW as declared in the (sub)program from which D01AMF is called. The number of sub-intervals into which the interval of integration may be divided cannot exceed LIW. Suggested value: LIW = LW/4. Constraint: LIW >= 1. 12: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. Users who are unfamiliar with this parameter should refer to the Essential Introduction for details. On exit: IFAIL = 0 unless the routine detects an error or gives a warning (see Section 6). For this routine, because the values of output parameters may be useful even if IFAIL /=0 on exit, users are recommended to set IFAIL to -1 before entry. It is then essential to test the value of IFAIL on exit. 6. Error Indicators and Warnings Errors or warnings specified by the routine: If on entry IFAIL = 0 or -1, explanatory error messages are output on the current error message unit (as defined by X04AAF). IFAIL= 1 The maximum number of subdivisions allowed with the given workspace has been reached without the requested accuracy requirements being achieved. Look at the integrand in order to determine the integration difficulties. If the position of a local difficulty within the interval can be determined (e.g. a singularity of the integrand or its derivative, a peak, a discontinuity, etc) you will probably gain from splitting up the interval at this point and calling D01AMF on the infinite subrange and an appropriate integrator on the finite subrange. Alternatively, consider relaxing the accuracy requirements specified by EPSABS and EPSREL, or increasing the amount of workspace. IFAIL= 2 Round-off error prevents the requested tolerance from being achieved. The error may be underestimated. Consider requesting less accuracy. IFAIL= 3 Extremely bad local integrand behaviour causes a very strong subdivision around one (or more) points of the interval. The same advice applies as in the case of IFAIL = 1. IFAIL= 4 The requested tolerance cannot be achieved, because the extrapolation does not increase the accuracy satisfactorily; the returned result is the best which can be obtained. The same advice applies as in the case of IFAIL = 1. IFAIL= 5 The integral is probably divergent, or slowly convergent. It must be noted that divergence can also occur with any other non-zero value of IFAIL. IFAIL= 6 On entry LW < 4, or LIW < 1, or INF /= -1, 1 or 2. Please note that divergence can occur with any non-zero value of IFAIL. 7. Accuracy The routine cannot guarantee, but in practice usually achieves, the following accuracy: |I-RESULT|<=tol, where tol=max{|EPSABS|,|EPSREL|*|I|}, and EPSABS and EPSREL are user-specified absolute and relative error tolerances. Moreover it returns the quantity ABSERR, which, in normal circumstances, satisfies |I-RESULT|<=ABSERR<=tol. 8. Further Comments The time taken by the routine depends on the integrand and the accuracy required. If IFAIL /= 0 on exit, then the user may wish to examine the contents of the array W, which contains the end-points of the sub-intervals used by D01AMF along with the integral contributions and error estimates over these sub-intervals. Specifically, for i=1,2,...,n, let r denote the approximation to i the value of the integral over the sub-interval [a ,b ] in the i i partition of [a,b] and e be the corresponding absolute error i b i n / -- estimate. Then, | f(x)dx~=r and RESULT= > r unless D01AMF / i -- i a i=1 i terminates while testing for divergence of the integral (see Piessens et al [3] Section 3.4.3). In this case, RESULT (and ABSERR) are taken to be the values returned from the extrapolation process. The value of n is returned in IW(1), and the values a , b , e and r are stored consecutively in the i i i i array W, that is: a =W(i),b =W(n+i),e =W(2n+i) and r =W(3n+i). i i i i Note: that this information applies to the integral transformed to (0,1) as described in Section 3, not to the original integral. 9. Example To compute infty / 1 | --------dx. / 0 (x+1)\/x The exact answer is (pi). The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXd01anf}{NAG On-line Documentation: d01anf} \beginscroll \begin{verbatim} D01ANF(3NAG) Foundation Library (12/10/92) D01ANF(3NAG) D01 -- Quadrature D01ANF D01ANF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose D01ANF calculates an approximation to the sine or the cosine transform of a function g over [a,b]: b b / / I= |g(x)||,sin((omega)x)||,dx or I= |g(x)||,cos((omega)x)||,dx / / a a (for a user-specified value of (omega)). 2. Specification SUBROUTINE D01ANF (G, A, B, OMEGA, KEY, EPSABS, EPSREL, 1 RESULT, ABSERR, W, LW, IW, LIW, IFAIL) INTEGER KEY, LW, IW(LIW), LIW, IFAIL DOUBLE PRECISION G, A, B, OMEGA, EPSABS, EPSREL, RESULT, 1 ABSERR, W(LW) EXTERNAL G 3. Description D01ANF is based upon the QUADPACK routine QFOUR (Piessens et al [3]). It is an adaptive routine, designed to integrate a function of the form g(x)w(x), where w(x) is either sin((omega)x) or cos((omega)x). If a sub-interval has length -l L=|b-a|2 then the integration over this sub-interval is performed by means of a modified Clenshaw-Curtis procedure (Piessens and Branders [2]) if L(omega)>4 and l<=20. In this case a Chebyshev-series approximation of degree 24 is used to approximate g(x), while an error estimate is computed from this approximation together with that obtained using Chebyshev-series of degree 12. If the above conditions do not hold then Gauss 7-point and Kronrod 15-point rules are used. The algorithm, described in [3], incorporates a global acceptance criterion (as defined in Malcolm and Simpson [1]) together with the (epsilon)-algorithm Wynn [4] to perform extrapolation. The local error estimation is described in [3]. 4. References [1] Malcolm M A and Simpson R B (1976) Local Versus Global Strategies for Adaptive Quadrature. ACM Trans. Math. Softw. 1 129--146. [2] Piessens R and Branders M (1975) Algorithm 002. Computation of Oscillating Integrals. J. Comput. Appl. Math. 1 153--164. [3] Piessens R, De Doncker-Kapenga E, Uberhuber C and Kahaner D (1983) QUADPACK, A Subroutine Package for Automatic Integration. Springer-Verlag. [4] Wynn P (1956) On a Device for Computing the e (S ) m n Transformation. Math. Tables Aids Comput. 10 91--96. 5. Parameters 1: G -- DOUBLE PRECISION FUNCTION, supplied by the user. External Procedure G must return the value of the function g at a given point. Its specification is: DOUBLE PRECISION FUNCTION G (X) DOUBLE PRECISION X 1: X -- DOUBLE PRECISION Input On entry: the point at which the function g must be evaluated. G must be declared as EXTERNAL in the (sub)program from which D01ANF is called. Parameters denoted as Input must not be changed by this procedure. 2: A -- DOUBLE PRECISION Input On entry: the lower limit of integration, a. 3: B -- DOUBLE PRECISION Input On entry: the upper limit of integration, b. It is not necessary that a= 4. 12: IW(LIW) -- INTEGER array Output On exit: IW(1) contains the actual number of sub-intervals used. The rest of the array is used as workspace. 13: LIW -- INTEGER Input On entry: the dimension of the array IW as declared in the (sub)program from which D01ANF is called. The number of sub-intervals into which the interval of integration may be divided cannot exceed LIW/2. Suggested value: LIW = LW/2. Constraint: LIW >= 2. 14: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. Users who are unfamiliar with this parameter should refer to the Essential Introduction for details. On exit: IFAIL = 0 unless the routine detects an error or gives a warning (see Section 6). For this routine, because the values of output parameters may be useful even if IFAIL /=0 on exit, users are recommended to set IFAIL to -1 before entry. It is then essential to test the value of IFAIL on exit. 6. Error Indicators and Warnings Errors or warnings specified by the routine: If on entry IFAIL = 0 or -1, explanatory error messages are output on the current error message unit (as defined by X04AAF). IFAIL= 1 The maximum number of subdivisions allowed with the given workspace has been reached without the accuracy requested being achieved. Look at the integrand in order to determine the integration difficulties. If the position of a local difficulty within the interval can be determined (e.g. a singularity of the integrand or its derivative, a peak, a discontinuity, etc) you will probably gain from splitting up the interval at this point and calling the integrator on the subranges. If necessary, another integrator, which is designed for handling the type of difficulty involved, must be used. Alternatively consider relaxing the accuracy requirements specified by EPSABS and EPSREL, or increasing amount of workspace. IFAIL= 2 Round-off error prevents the requested tolerance from being achieved. The error may be underestimated. Consider requesting less accuracy. IFAIL= 3 Extremely bad local behaviour of g(x) causes a very strong subdivision around one (or more) points of the interval. The same advice applies as in the case of IFAIL = 1. IFAIL= 4 The requested tolerance cannot be achieved because the extrapolation does not increase the accuracy satisfactorily; the returned result is the best which can be obtained. The same advice applies as in the case of IFAIL = 1. IFAIL= 5 The integral is probably divergent, or slowly convergent. It must be noted that divergence can occur with any non-zero value of IFAIL. IFAIL= 6 On entry KEY < 1, or KEY > 2. IFAIL= 7 On entry LW < 4, or LIW < 2. 7. Accuracy The routine cannot guarantee, but in practice usually achieves, the following accuracy: |I-RESULT|<=tol, where tol=max{|EPSABS|,|EPSREL|*|I|}, and EPSABS and EPSREL are user-specified absolute and relative tolerances. Moreover it returns the quantity ABSERR, which, in normal circumstances, satisfies |I-RESULT|<=ABSERR<=tol. 8. Further Comments The time taken by the routine depends on the integrand and on the accuracy required. If IFAIL /= 0 on exit, then the user may wish to examine the contents of the array W, which contains the end-points of the sub-intervals used by D01ANF along with the integral contributions and error estimates over these sub-intervals. Specifically, for i=1,2,...,n, let r denote the approximation to i the value of the integral over the sub-interval [a ,b ] in the i i partition of [a,b] and e be the corresponding absolute error i b i n / -- estimate. Then, | g(x)w(x)dx~=r and RESULT= > r unless D01ANF / i -- i a i=1 i terminates while testing for divergence of the integral (see Piessens et al [3] Section 3.4.3). In this case, RESULT (and ABSERR) are taken to be the values returned from the extrapolation process. The value of n is returned in IW(1), and the values a , b , e and r are stored consecutively in the i i i i array W, that is: a =W(i), i b =W(n+i), i e =W(2n+i) and i r =W(3n+i). i 9. Example To compute 1 / |ln(x)sin(10(pi)x)dx. / 0 The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXd01apf}{NAG On-line Documentation: d01apf} \beginscroll \begin{verbatim} D01APF(3NAG) Foundation Library (12/10/92) D01APF(3NAG) D01 -- Quadrature D01APF D01APF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose D01APF is an adaptive integrator which calculates an approximation to the integral of a function g(x)w(x) over a finite interval [a,b]: b / I= |g(x)w(x)dx / a where the weight function w has end-point singularities of algebraico-logarithmic type. 2. Specification SUBROUTINE D01APF (G, A, B, ALFA, BETA, KEY, EPSABS, 1 EPSREL, RESULT, ABSERR, W, LW, IW, LIW, 2 IFAIL) INTEGER KEY, LW, IW(LIW), LIW, IFAIL DOUBLE PRECISION G, A, B, ALFA, BETA, EPSABS, EPSREL, 1 RESULT, ABSERR, W(LW) EXTERNAL G 3. Description D01APF is based upon the QUADPACK routine QAWSE (Piessens et al [3]) and integrates a function of the form g(x)w(x), where the weight function w(x) may have algebraico-logarithmic singularities at the end-points a and/or b. The strategy is a modification of that in D01AKF. We start by bisecting the original interval and applying modified Clenshaw-Curtis integration of orders 12 and 24 to both halves. Clenshaw-Curtis integration is then used on all sub-intervals which have a or b as one of their end-points (Piessens et al [2]). On the other sub-intervals Gauss-Kronrod (7-15 point) integration is carried out. A 'global' acceptance criterion (as defined by Malcolm and Simpson [1]) is used. The local error estimation control is described by Piessens et al [3]. 4. References [1] Malcolm M A and Simpson R B (1976) Local Versus Global Strategies for Adaptive Quadrature. ACM Trans. Math. Softw. 1 129--146. [2] Piessens R, Mertens I and Branders M (1974) Integration of Functions having End-point Singularities. Angewandte Informatik. 16 65--68. [3] Piessens R, De Doncker-Kapenga E, Uberhuber C and Kahaner D (1983) QUADPACK, A Subroutine Package for Automatic Integration. Springer-Verlag. 5. Parameters 1: G -- DOUBLE PRECISION FUNCTION, supplied by the user. External Procedure G must return the value of the function g at a given point X. Its specification is: DOUBLE PRECISION FUNCTION G (X) DOUBLE PRECISION X 1: X -- DOUBLE PRECISION Input On entry: the point at which the function g must be evaluated. G must be declared as EXTERNAL in the (sub)program from which D01APF is called. Parameters denoted as Input must not be changed by this procedure. 2: A -- DOUBLE PRECISION Input On entry: the lower limit of integration, a. 3: B -- DOUBLE PRECISION Input On entry: the upper limit of integration, b. Constraint: B > A. 4: ALFA -- DOUBLE PRECISION Input On entry: the parameter (alpha) in the weight function. Constraint: ALFA>-1. 5: BETA -- DOUBLE PRECISION Input On entry: the parameter (beta) in the weight function. Constraint: BETA>-1. 6: KEY -- INTEGER Input On entry: indicates which weight function is to be used: (alpha) (beta) if KEY = 1, w(x)=(x-a) (b-x) (alpha) (beta) if KEY = 2, w(x)=(x-a) (b-x) ln(x-a) (alpha) (beta) if KEY = 3, w(x)=(x-a) (b-x) ln(b-x) (alpha) (beta) if KEY = 4, w(x)=(x-a) (b-x) ln(x-a)ln(b-x) Constraint: KEY = 1, 2, 3 or 4 7: EPSABS -- DOUBLE PRECISION Input On entry: the absolute accuracy required. If EPSABS is negative, the absolute value is used. See Section 7. 8: EPSREL -- DOUBLE PRECISION Input On entry: the relative accuracy required. If EPSREL is negative, the absolute value is used. See Section 7. 9: RESULT -- DOUBLE PRECISION Output On exit: the approximation to the integral I. 10: ABSERR -- DOUBLE PRECISION Output On exit: an estimate of the modulus of the absolute error, which should be an upper bound for |I-RESULT|. 11: W(LW) -- DOUBLE PRECISION array Output On exit: details of the computation, as described in Section 8. 12: LW -- INTEGER Input On entry: the dimension of the array W as declared in the (sub)program from which D01APF is called. The value of LW (together with that of LIW below) imposes a bound on the number of sub-intervals into which the interval of integration may be divided by the routine. The number of sub-intervals cannot exceed LW/4. The more difficult the integrand, the larger LW should be. Suggested value: LW = 800 to 2000 is adequate for most problems. Constraint: LW >= 8. 13: IW(LIW) -- INTEGER array Output On exit: IW(1) contains the actual number of sub-intervals used. The rest of the array is used as workspace. 14: LIW -- INTEGER Input On entry: the dimension of the array IW as declared in the (sub)program from which D01APF is called. The number of sub-intervals into which the interval of integration may be divided cannot exceed LIW. Suggested value: LIW = LW/4. Constraint: LIW >= 2. 15: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. Users who are unfamiliar with this parameter should refer to the Essential Introduction for details. On exit: IFAIL = 0 unless the routine detects an error or gives a warning (see Section 6). For this routine, because the values of output parameters may be useful even if IFAIL /=0 on exit, users are recommended to set IFAIL to -1 before entry. It is then essential to test the value of IFAIL on exit. 6. Error Indicators and Warnings Errors or warnings specified by the routine: If on entry IFAIL = 0 or -1, explanatory error messages are output on the current error message unit (as defined by X04AAF). IFAIL= 1 The maximum number of subdivisions allowed with the given workspace has been reached without the accuracy requirements being achieved. Look at the integrand in order to determine the integration difficulties. If the position of a discontinuity or a singularity of algebraico-logarithmic type within the interval can be determined, the interval must be split up at this point and the integrator called on the subranges. If necessary, another integrator, which is designed for handling the difficulty involved, must be used. Alternatively consider relaxing the accuracy requirements specified by EPSABS and EPSREL, or increasing the amount of workspace. IFAIL= 2 Round-off error prevents the requested tolerance from being achieved. Consider requesting less accuracy. IFAIL= 3 Extremely bad local integrand behaviour causes a very strong subdivision around one (or more) points of the interval. The same advice applies as in the case of IFAIL = 1. IFAIL= 4 On entry B <= A, or ALFA <= -1, or BETA <= -1, or KEY < 1, or KEY > 4. IFAIL= 5 On entry LW < 8, or LIW < 2. 7. Accuracy The routine cannot guarantee, but in practice usually achieves, the following accuracy: |I-RESULT|<=tol, where tol=max{|EPSABS|,|EPSREL|*|I|}, and EPSABS and EPSREL are user-specified absolute and relative error tolerances. Moreover it returns the quantity ABSERR which, in normal circumstances, satisfies: |I-RESULT|<=ABSERR<=tol. 8. Further Comments The time taken by the routine depends on the integrand and on the accuracy required. If IFAIL /= 0 on exit, then the user may wish to examine the contents of the array W, which contains the end-points of the sub-intervals used by D01APF along with the integral contributions and error estimates over these sub-intervals. Specifically, for i=1,2,...,n, let r denote the approximation to i the value of the integral over the sub-interval [a ,b ] in the i i partition of [a,b] and e be the corresponding absolute error i b i n / -- estimate. Then, | f(x)w(x)dx~=r and RESULT= > r . The value of / i -- i a i=1 i n is returned in IW(1), and the values a , b , e and r are i i i i stored consecutively in the array W, that is: a =W(i), i b =W(n+i), i e =W(2n+i), i r =W(3n+i). i 9. Example To compute: 1 / |ln(x)cos(10(pi)x)dx and / 0 1 / sin(10x) | --------dx. / 0 \/x(1-x) The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXd01aqf}{NAG On-line Documentation: d01aqf} \beginscroll \begin{verbatim} D01AQF(3NAG) Foundation Library (12/10/92) D01AQF(3NAG) D01 -- Quadrature D01AQF D01AQF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose D01AQF calculates an approximation to the Hilbert transform of a function g(x) over [a,b]: b / g(x) I= | ----dx / x-c a for user-specified values of a, b and c. 2. Specification SUBROUTINE D01AQF (G, A, B, C, EPSABS, EPSREL, RESULT, 1 ABSERR, W, LW, IW, LIW, IFAIL) INTEGER LW, IW(LIW), LIW, IFAIL DOUBLE PRECISION G, A, B, C, EPSABS, EPSREL, RESULT, 1 ABSERR, W(LW) EXTERNAL G 3. Description D01AQF is based upon the QUADPACK routine QAWC (Piessens et al [3]) and integrates a function of the form g(x)w(x), where the weight function 1 w(x)= --- x-c is that of the Hilbert transform. (If a= 4. 11: IW(LIW) -- INTEGER array Output On exit: IW(1) contains the actual number of sub-intervals used. The rest of the array is used as workspace. 12: LIW -- INTEGER Input On entry: the dimension of the array IW as declared in the (sub)program from which D01AQF is called. The number of sub-intervals into which the interval of integration may be divided cannot exceed LIW. Suggested value: LIW = LW/4. Constraint: LIW >= 1. 13: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. Users who are unfamiliar with this parameter should refer to the Essential Introduction for details. On exit: IFAIL = 0 unless the routine detects an error or gives a warning (see Section 6). For this routine, because the values of output parameters may be useful even if IFAIL /=0 on exit, users are recommended to set IFAIL to -1 before entry. It is then essential to test the value of IFAIL on exit. 6. Error Indicators and Warnings Errors or warnings specified by the routine: If on entry IFAIL = 0 or -1, explanatory error messages are output on the current error message unit (as defined by X04AAF). IFAIL= 1 The maximum number of subdivisions allowed with the given workspace has been reached without the accuracy requirements being achieved. Look at the integrand in order to determine the integration difficulties. Another integrator which is designed for handling the type of difficulty involved, must be used. Alternatively consider relaxing the accuracy requirements specified by EPSABS and EPSREL, or increasing the workspace. IFAIL= 2 Round-off error prevents the requested tolerance from being achieved. Consider requesting less accuracy. IFAIL= 3 Extremely bad local behaviour of g(x) causes a very strong subdivision around one (or more) points of the interval. The same advice applies as in the case of IFAIL = 1. IFAIL= 4 On entry C = A or C = B. IFAIL= 5 On entry LW < 4, or LIW < 1. 7. Accuracy The routine cannot guarantee, but in practice usually achieves, the following accuracy: |I-RESULT|<=tol, where tol=max{|EPSABS|,|EPSREL|*|I|}, and EPSABS and EPSREL are user-specified absolute and relative error tolerances. Moreover it returns the quantity ABSERR which, in normal circumstances satisfies: |I-RESULT|<=ABSERR<=tol. 8. Further Comments The time taken by the routine depends on the integrand and on the accuracy required. If IFAIL /= 0 on exit, then the user may wish to examine the contents of the array W, which contains the end-points of the sub-intervals used by D01AQF along with the integral contributions and error estimates over these sub-intervals. Specifically, for i=1,2,...,n, let r denote the approximation to i the value of the integral over the sub-interval [a ,b ] in the i i partition of [a,b] and e be the corresponding absolute error i b i n / -- estimate. Then, | g(x)w(x)dx~=r and RESULT= > r . The value of / i -- i a i=1 i n is returned in IW(1), and the values a , b , e and r are i i i i stored consecutively in the array W, that is: a =W(i), i b =W(n+i), i e =W(2n+i) and i r =W(3n+i). i 9. Example To compute the Cauchy principal value of 1 / dx | ----------------. / 2 2 1 -1 (x +0.01 )(x- -) 2 The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXd01asf}{NAG On-line Documentation: d01asf} \beginscroll \begin{verbatim} D01ASF(3NAG) Foundation Library (12/10/92) D01ASF(3NAG) D01 -- Quadrature D01ASF D01ASF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose D01ASF calculates an approximation to the sine or the cosine transform of a function g over [a,infty): infty / I= | g(x)||,sin((omega)x)dx or I= / a infty / | g(x)||,cos((omega)x)dx / a (for a user-specified value of (omega)). 2. Specification SUBROUTINE D01ASF (G, A, OMEGA, KEY, EPSABS, RESULT, 1 ABSERR, LIMLST, LST, ERLST, RSLST, 2 IERLST, W, LW, IW, LIW, IFAIL) INTEGER KEY, LIMLST, LST, IERLST(LIMLST), LW, IW 1 (LIW), LIW, IFAIL DOUBLE PRECISION G, A, OMEGA, EPSABS, RESULT, ABSERR, ERLST 1 (LIMLST), RSLST(LIMLST), W(LW) EXTERNAL G 3. Description D01ASF is based upon the QUADPACK routine QAWFE (Piessens et al [2]). It is an adaptive routine, designed to integrate a function of the form g(x)w(x) over a semi-infinite interval, where w(x) is either sin((omega)x) or cos((omega)x). Over successive intervals C =[a+(k-1)c,a+kc], k=1,2,...,LST k integration is performed by the same algorithm as is used by D01ANF. The intervals C are of constant length k c={2[|(omega)|]+1}(pi)/|(omega)|, (omega)/=0 where [|(omega)|] represents the largest integer less than or equal to |(omega)|. Since c equals an odd number of half periods, the integral contributions over succeeding intervals will alternate in sign when the function g is positive and monotonically decreasing over [a,infty). The algorithm, described by [2], incorporates a global acceptance criterion (as defined by Malcolm and Simpson [1]) together with the (epsilon)-algorithm (Wynn [3]) to perform extrapolation. The local error estimation is described by Piessens et al [2]. If (omega)=0 and KEY = 1, the routine uses the same algorithm as D01AMF (with EPSREL = 0.0). In contrast to the other routines in Chapter D01, D01ASF works only with a user-specified absolute error tolerance (EPSABS). Over the interval C it attempts to satisfy the absolute accuracy k requirement EPSA =U *EPSABS k k k-1 where U =(1-p)p , for k=1,2,... and p=0.9. k However, when difficulties occur during the integration over the kth sub-interval C such that the error flag IERLST(k) is non- k zero, the accuracy requirement over subsequent intervals is relaxed. See Piessens et al [2] for more details. 4. References [1] Malcolm M A and Simpson R B (1976) Local Versus Global Strategies for Adaptive Quadrature. ACM Trans. Math. Softw. 1 129--146. [2] Piessens R, De Doncker-Kapenga E, Uberhuber C and Kahaner D (1983) QUADPACK, A Subroutine Package for Automatic Integration. Springer-Verlag. [3] Wynn P (1956) On a Device for Computing the e (S ) m n Transformation. Math. Tables Aids Comput. 10 91--96. 5. Parameters 1: G -- DOUBLE PRECISION FUNCTION, supplied by the user. External Procedure G must return the value of the function g at a given point. Its specification is: DOUBLE PRECISION FUNCTION G (X) DOUBLE PRECISION X 1: X -- DOUBLE PRECISION Input On entry: the point at which the function g must be evaluated. G must be declared as EXTERNAL in the (sub)program from which D01ASF is called. Parameters denoted as Input must not be changed by this procedure. 2: A -- DOUBLE PRECISION Input On entry: the lower limit of integration, a. 3: OMEGA -- DOUBLE PRECISION Input On entry: the parameter (omega) in the weight function of the transform. 4: KEY -- INTEGER Input On entry: indicates which integral is to be computed: if KEY = 1, w(x)=cos((omega)x); if KEY = 2, w(x)=sin((omega)x). Constraint: KEY = 1 or 2. 5: EPSABS -- DOUBLE PRECISION Input On entry: the absolute accuracy requirement. If EPSABS is negative, the absolute value is used. See Section 7. 6: RESULT -- DOUBLE PRECISION Output On exit: the approximation to the integral I. 7: ABSERR -- DOUBLE PRECISION Output On exit: an estimate of the modulus of the absolute error, which should be an upper bound for |I-RESULT|. 8: LIMLST -- INTEGER Input On entry: an upper bound on the number of intervals C k needed for the integration. Suggested value: LIMLST = 50 is adequate for most problems. Constraint: LIMLST >= 3. 9: LST -- INTEGER Output On exit: the number of intervals C actually used for the k integration. 10: ERLST(LIMLST) -- DOUBLE PRECISION array Output On exit: ERLST(k) contains the error estimate corresponding to the integral contribution over the interval C , for k k=1,2,...,LST. 11: RSLST(LIMLST) -- DOUBLE PRECISION array Output On exit: RSLST(k) contains the integral contribution over the interval C for k=1,2,...,LST. k 12: IERLST(LIMLST) -- INTEGER array Output On exit: IERLST(k) contains the error flag corresponding to RSLST(k), for k=1,2,...,LST. See Section 6. 13: W(LW) -- DOUBLE PRECISION array Workspace 14: LW -- INTEGER Input On entry: the dimension of the array W as declared in the (sub)program from which D01ASF is called. The value of LW (together with that of LIW below) imposes a bound on the number of sub-intervals into which each interval C may be divided by the routine. The number of k sub-intervals cannot exceed LW/4. The more difficult the integrand, the larger LW should be. Suggested value: a value in the range 800 to 2000 is adequate for most problems. Constraint: LW >= 4. 15: IW(LIW) -- INTEGER array Output On exit: IW(1) contains the maximum number of sub-intervals actually used for integrating over any of the intervals C . k The rest of the array is used as workspace. 16: LIW -- INTEGER Input On entry: the dimension of the array IW as declared in the (sub)program from which D01ASF is called. The number of sub-intervals into which each interval C may k be divided cannot exceed LIW/2. Suggested value: LIW = LW/2. Constraint: LIW >= 2. 17: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. Users who are unfamiliar with this parameter should refer to the Essential Introduction for details. On exit: IFAIL = 0 unless the routine detects an error or gives a warning (see Section 6). For this routine, because the values of output parameters may be useful even if IFAIL /=0 on exit, users are recommended to set IFAIL to -1 before entry. It is then essential to test the value of IFAIL on exit. 6. Error Indicators and Warnings Errors or warnings specified by the routine: If on entry IFAIL = 0 or -1, explanatory error messages are output on the current error message unit (as defined by X04AAF). IFAIL= 1 The maximum number of subdivisions allowed with the given workspace has been reached without the accuracy requirements being achieved. Look at the integrand in order to determine the integration difficulties. If the position of a local difficulty within the interval can be determined (e.g. a singularity of the integrand or its derivative, a peak, a discontinuity, etc) you will probably gain from splitting up the interval at this point and calling D01ASF on the infinite subrange and an appropriate integrator on the finite subrange. Alternatively, consider relaxing the accuracy requirements specified by EPSABS or increasing the amount of workspace. IFAIL= 2 Round-off error prevents the requested tolerance from being achieved. The error may be underestimated. Consider requesting less accuracy. IFAIL= 3 Extremely bad local integrand behaviour causes a very strong subdivision around one (or more) points of the interval. The same advice applies as in the case of IFAIL = 1. IFAIL= 4 The requested tolerance cannot be achieved, because the extrapolation does not increase the accuracy satisfactorily; the returned result is the best which can be obtained. The same advice applies as in the case of IFAIL = 1. IFAIL= 5 The integral is probably divergent, or slowly convergent. Please note that divergence can occur with any non-zero value of IFAIL. IFAIL= 6 On entry KEY < 1, or KEY > 2, or LIMLST < 3. IFAIL= 7 Bad integration behaviour occurs within one or more of the intervals C . Location and type of the difficulty involved k can be determined from the vector IERLST (see below). IFAIL= 8 Maximum number of intervals C (= LIMLST) allowed has been k achieved. Increase the value of LIMLST to allow more cycles. IFAIL= 9 The extrapolation table constructed for convergence acceleration of the series formed by the integral contribution over the intervals C , does not converge to the k required accuracy. IFAIL= 10 On entry LW < 4, or LIW < 2. In the cases IFAIL = 7, 8 or 9, additional information about the cause of the error can be obtained from the array IERLST, as follows: IERLST(k)=1 The maximum number of subdivisions = min(LW/4,LIW/2) has been achieved on the kth interval. IERLST(k)=2 Occurrence of round-off error is detected and prevents the tolerance imposed on the kth interval from being achieved. IERLST(k)=3 Extremely bad integrand behaviour occurs at some points of the kth interval. IERLST(k)=4 The integration procedure over the kth interval does not converge (to within the required accuracy) due to round-off in the extrapolation procedure invoked on this interval. It is assumed that the result on this interval is the best which can be obtained. IERLST(k)=5 The integral over the kth interval is probably divergent or slowly convergent. It must be noted that divergence can occur with any other value of IERLST(k). 7. Accuracy The routine cannot guarantee, but in practice usually achieves, the following accuracy: |I-RESULT|<=|EPSABS|, where EPSABS is the user-specified absolute error tolerance. Moreover, it returns the quantity ABSERR, which, in normal circumstances, satisfies |I-RESULT|<=ABSERR<=|EPSABS|. 8. Further Comments None. 9. Example To compute infty / 1 | ---cos((pi)x/2)dx. / 0 \/x The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXd01bbf}{NAG On-line Documentation: d01bbf} \beginscroll \begin{verbatim} D01BBF(3NAG) D01BBF D01BBF(3NAG) D01 -- Quadrature D01BBF D01BBF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. Note for users via the AXIOM system: the interface to this routine has been enhanced for use with AXIOM and is slightly different to that offered in the standard version of the Foundation Library. 1. Purpose D01BBF returns the weights and abscissae appropriate to a Gaussian quadrature formula with a specified number of abscissae. The formulae provided are Gauss-Legendre, Gauss-Rational, Gauss- Laguerre and Gauss-Hermite. 2. Specification SUBROUTINE D01BBF (A, B, ITYPE, N, WEIGHT, ABSCIS, GTYPE, 1 IFAIL) INTEGER ITYPE, N, GTYPE, IFAIL DOUBLE PRECISION A, B, WEIGHT(N), ABSCIS(N) 3. Description This routine returns the weights and abscissae for use in the Gaussian quadrature of a function f(x). The quadrature takes the form n -- S= > w f(x ) -- i i i=1 where w are the weights and x are the abscissae (see Davis and i i Rabinowitz [1], Froberg [2], Ralston [3] or Stroud and Secrest [4]). Weights and abscissae are available for Gauss-Legendre, Gauss- Rational, Gauss-Laguerre and Gauss-Hermite quadrature, and for a selection of values of n (see Section 5). (a) Gauss-Legendre Quadrature: b / S~= |f(x)dx / a where a and b are finite and it will be exact for any function of the form 2n-1 -- i f(x)= > c x -- i i=0 (b) Gauss-Rational quadrature: infty a / / S~= | f(x)dx (a+b>0) or S~= | f(x)dx (a+b<0) / / a -infty and will be exact for any function of the form 2n-1 -- i > c (x+b) 2n+1 c -- 2n+1-i -- i i=0 f(x)= > ------= ------------------ -- i 2n+1 i=2 (x+b) (x+b) (c) Gauss-Laguerre quadrature, adjusted weights option: infty a / / S~= | f(x)dx (b>0<) or S~= | f(x)dx (b<0) / / a -infty and will be exact for any function of the form 2n-1 -bx -- i f(x)=e > c x -- i i=0 (d) Gauss-Hermite quadrature, adjusted weights option: +infty / S~= | f(x)dx / -infty and will be exact for any function of the form 2 2n-1 -b(x-a) -- i f(x)=e > c x (b>0) -- i i=0 (e) Gauss-Laguerre quadrature, normal weights option: infty a / -bx / -bx S~= | e f(x)dx (b>0) or S~= | e f(x)dx / / a -infty (b<0) and will be exact for any function of the form 2n-1 -- i f(x)= > c x -- i i=0 (f) Gauss-Hermite quadrature, normal weights option: +infty 2 / -b(x-a) S~= | e f(x)dx / -infty and will be exact for any function of the form: 2n-1 -- i f(x)= > c x -- i i=0 Note: that the Gauss-Legendre abscissae, with a=-1, b=+1, are the zeros of the Legendre polynomials; the Gauss-Laguerre abscissae, with a=0, b=1, are the zeros of the Laguerre polynomials; and the Gauss-Hermite abscissae, with a=0, b=1, are the zeros of the Hermite polynomials. 4. References [1] Davis P J and Rabinowitz P (1967) Numerical Integration. Blaisdell Publishing Company. 33--52. [2] Froberg C E (1965) Introduction to Numerical Analysis. Addison-Wesley. 181--187. [3] Ralston A (1965) A First Course in Numerical Analysis. McGraw-Hill. 87--90. [4] Stroud A H and Secrest D (1966) Gaussian Quadrature Formulas. Prentice-Hall. 5. Parameters 1: A -- DOUBLE PRECISION Input 2: B -- DOUBLE PRECISION Input On entry: the quantities a and b as described in the appropriate subsection of Section 3. 3: ITYPE -- INTEGER Input On entry: indicates the type of weights for Gauss-Laguerre or Gauss-Hermite quadrature (see Section 3): if ITYPE = 1, adjusted weights will be returned; if ITYPE = 0, normal weights will be returned. Constraint: ITYPE = 0 or 1. For Gauss-Legendre or Gauss-Rational quadrature, this parameter is not used. 4: N -- INTEGER Input On entry: the number of weights and abscissae to be returned, n. Constraint: N = 1,2,3,4,5,6,8,10,12,14,16,20,24,32,48 or 64. 5: WEIGHT(N) -- DOUBLE PRECISION array Output On exit: the N weights. For Gauss-Laguerre and Gauss- Hermite quadrature, these will be the adjusted weights if ITYPE = 1, and the normal weights if ITYPE = 0. 6: ABSCIS(N) -- DOUBLE PRECISION array Output On exit: the N abscissae. 7: GTYPE -- INTEGER Input On entry: The value of GTYPE indicates which quadrature formula are to be used: GTYPE = 0 for Gauss-Legendre weights and abscissae; GTYPE = 1 for Gauss-Rational weights and abscissae; GTYPE = 2 for Gauss-Laguerre weights and abscissae; GTYPE = 3 for Gauss-Hermite weights and abscissae. 8: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. On exit: IFAIL = 0 unless the routine detects an error (see Section 6). 6. Error Indicators and Warnings Errors detected by the routine: IFAIL= 1 The N-point rule is not among those stored. If the soft fail option is used, the weights and abscissae returned will be those for the largest valid value of N less than the requested value, and the excess elements of WEIGHT and ABSCIS (i.e., up to the requested N) will be filled with zeros. IFAIL= 2 The value of A and/or B is invalid. Gauss-Rational: A + B = 0 Gauss-Laguerre: B = 0 Gauss-Hermite: B <= 0 If the soft fail option is used the weights and abscissae are returned as zero. IFAIL= 3 Laguerre and Hermite normal weights only: underflow is occurring in evaluating one or more of the normal weights. If the soft fail option is used, the underflowing weights are returned as zero. A smaller value of N must be used; or adjusted weights should be used (ITYPE = 1). In the latter case, take care that underflow does not occur when evaluating the integrand appropriate for adjusted weights. IFAIL=4 GTYPE < 0 or GTYPE > 3 7. Accuracy The weights and abscissae are stored for standard values of A and B to full machine accuracy. 8. Further Comments Timing is negligible. 9. Example This example program returns the abscissae and (adjusted) weights for the six-point Gauss-Laguerre formula. The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXd01fcf}{NAG On-line Documentation: d01fcf} \beginscroll \begin{verbatim} D01FCF(3NAG) Foundation Library (12/10/92) D01FCF(3NAG) D01 -- Quadrature D01FCF D01FCF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose D01FCF attempts to evaluate a multi-dimensional integral (up to 15 dimensions), with constant and finite limits, to a specified relative accuracy, using an adaptive subdivision strategy. 2. Specification SUBROUTINE D01FCF (NDIM, A, B, MINPTS, MAXPTS, FUNCTN, 1 EPS, ACC, LENWRK, WRKSTR, FINVAL, IFAIL) INTEGER NDIM, MINPTS, MAXPTS, LENWRK, IFAIL DOUBLE PRECISION A(NDIM), B(NDIM), FUNCTN, EPS, ACC, WRKSTR 1 (LENWRK), FINVAL EXTERNAL FUNCTN 3. Description The routine returns an estimate of a multi-dimensional integral over a hyper-rectangle (i.e., with constant limits), and also an estimate of the relative error. The user sets the relative accuracy required, supplies the integrand as a function subprogram (FUNCTN), and also sets the minimum and maximum acceptable number of calls to FUNCTN (in MINPTS and MAXPTS). The routine operates by repeated subdivision of the hyper- rectangular region into smaller hyper-rectangles. In each subregion, the integral is estimated using a seventh-degree rule, and an error estimate is obtained by comparison with a fifth- degree rule which uses a subset of the same points. The fourth differences of the integrand along each co-ordinate axis are evaluated, and the subregion is marked for possible future subdivision in half along that co-ordinate axis which has the largest absolute fourth difference. If the estimated errors, totalled over the subregions, exceed the requested relative error (or if fewer than MINPTS calls to FUNCTN have been made), further subdivision is necessary, and is performed on the subregion with the largest estimated error, that subregion being halved along the appropriate co-ordinate axis. The routine will fail if the requested relative error level has not been attained by the time MAXPTS calls to FUNCTN have been made; or, if the amount LENWRK of working storage is insufficient. A formula for the recommended value of LENWRK is given in Section 5. If a smaller value is used, and is exhausted in the course of execution, the routine switches to a less efficient mode of operation; only if this mode also breaks down is insufficient storage reported. D01FCF is based on the HALF subroutine developed by van Dooren and de Ridder [1]. It uses a different basic rule, described by Genz and Malik [2]. 4. References [1] Van Dooren P and De Ridder L (1976) An Adaptive Algorithm for Numerical Integration over an N-dimensional Cube. J. Comput. Appl. Math. 2 207--217. [2] Genz A C and Malik A A (1980) An Adaptive Algorithm for Numerical Integration over an N-dimensional Rectangular Region. J. Comput. Appl. Math. 6 295--302. 5. Parameters 1: NDIM -- INTEGER Input On entry: the number of dimensions of the integral, n. Constraint: 2 <= NDIM <= 15. 2: A(NDIM) -- DOUBLE PRECISION array Input On entry: the lower limits of integration, a , for i i=1,2,...,n. 3: B(NDIM) -- DOUBLE PRECISION array Input On entry: the upper limits of integration, b , for i i=1,2,...,n. 4: MINPTS -- INTEGER Input/Output On entry: MINPTS must be set to the minimum number of integrand evaluations to be allowed. On exit: MINPTS contains the actual number of integrand evaluations used by D01FCF. 5: MAXPTS -- INTEGER Input On entry: the maximum number of integrand evaluations to be allowed. Constraints: MAXPTS >= MINPTS MAXPTS >= (alpha), NDIM 2 where (alpha)=2 +2*NDIM +2*NDIM+1. 6: FUNCTN -- DOUBLE PRECISION FUNCTION, supplied by the user. External Procedure FUNCTN must return the value of the integrand f at a given point. Its specification is: DOUBLE PRECISION FUNCTION FUNCTN (NDIM,Z) INTEGER NDIM DOUBLE PRECISION Z(NDIM) 1: NDIM -- INTEGER Input On entry: the number of dimensions of the integral, n. 2: Z(NDIM) -- DOUBLE PRECISION array Input On entry: the co-ordinates of the point at which the integrand must be evaluated. FUNCTN must be declared as EXTERNAL in the (sub)program from which D01FCF is called. Parameters denoted as Input must not be changed by this procedure. 7: EPS -- DOUBLE PRECISION Input On entry: the relative error acceptable to the user. When the solution is zero or very small relative accuracy may not be achievable but the user may still set EPS to a reasonable value and check for the error exit IFAIL = 2. Constraint: EPS > 0.0. 8: ACC -- DOUBLE PRECISION Output On exit: the estimated relative error in FINVAL. 9: LENWRK -- INTEGER Input On entry: the dimension of the array WRKSTR as declared in the (sub)program from which D01FCF is called. Suggested value: for maximum efficiency, LENWRK >= (NDIM+2)*(1+MAXPTS/(alpha)) (see parameter MAXPTS for (alpha)). If LENWRK is less than this, the routine will usually run less efficiently and may fail. Constraint: LENWRK>=2*NDIM+4. 10: WRKSTR(LENWRK) -- DOUBLE PRECISION array Workspace 11: FINVAL -- DOUBLE PRECISION Output On exit: the best estimate obtained for the integral. 12: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. Users who are unfamiliar with this parameter should refer to the Essential Introduction for details. On exit: IFAIL = 0 unless the routine detects an error or gives a warning (see Section 6). For this routine, because the values of output parameters may be useful even if IFAIL /=0 on exit, users are recommended to set IFAIL to -1 before entry. It is then essential to test the value of IFAIL on exit. To suppress the output of an error message when soft failure occurs, set IFAIL to 1. 6. Error Indicators and Warnings Errors or warnings specified by the routine: IFAIL= 1 On entry NDIM < 2, or NDIM > 15, or MAXPTS is too small, or LENWRK<2*NDIM+4, or EPS <= 0.0. IFAIL= 2 MAXPTS was too small to obtain the required relative accuracy EPS. On soft failure, FINVAL and ACC contain estimates of the integral and the relative error, but ACC will be greater than EPS. IFAIL= 3 LENWRK was too small. On soft failure, FINVAL and ACC contain estimates of the integral and the relative error, but ACC will be greater than EPS. 7. Accuracy A relative error estimate is output through the parameter ACC. 8. Further Comments Execution time will usually be dominated by the time taken to evaluate the integrand FUNCTN, and hence the maximum time that could be taken will be proportional to MAXPTS. 9. Example This example program estimates the integral 2 1 1 1 1 4z z exp(2z z ) / / / / 1 3 1 3 | | | | ---------------dz dz dz dz =0.575364. / / / / 2 4 3 2 1 0 0 0 0 (1+z +z ) 2 4 The accuracy requested is one part in 10,000. The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXd01gaf}{NAG On-line Documentation: d01gaf} \beginscroll \begin{verbatim} D01GAF(3NAG) Foundation Library (12/10/92) D01GAF(3NAG) D01 -- Quadrature D01GAF D01GAF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose D01GAF integrates a function which is specified numerically at four or more points, over the whole of its specified range, using third-order finite-difference formulae with error estimates, according to a method due to Gill and Miller. 2. Specification SUBROUTINE D01GAF (X, Y, N, ANS, ER, IFAIL) INTEGER N, IFAIL DOUBLE PRECISION X(N), Y(N), ANS, ER 3. Description This routine evaluates the definite integral x n / I= | y(x)dx, / x 1 where the function y is specified at the n-points x ,x ,...,x , 1 2 n which should be all distinct, and in either ascending or descending order. The integral between successive points is calculated by a four-point finite-difference formula centred on the interval concerned, except in the case of the first and last intervals, where four-point forward and backward difference formulae respectively are employed. If n is less than 4, the routine fails. An approximation to the truncation error is integrated and added to the result. It is also returned separately to give an estimate of the uncertainty in the result. The method is due to Gill and Miller. 4. References [1] Gill P E and Miller G F (1972) An Algorithm for the Integration of Unequally Spaced Data. Comput. J. 15 80--83. 5. Parameters 1: X(N) -- DOUBLE PRECISION array Input On entry: the values of the independent variable, i.e., the x ,x ,...,x . Constraint: either X(1) < X(2) <... < X(N) or 1 2 n X(1) > X(2) >... > X(N). 2: Y(N) -- DOUBLE PRECISION array Input On entry: the values of the dependent variable y at the i points x , for i=1,2,...,n. i 3: N -- INTEGER Input On entry: the number of points, n. Constraint: N >= 4. 4: ANS -- DOUBLE PRECISION Output On exit: the estimate of the integral. 5: ER -- DOUBLE PRECISION Output On exit: an estimate of the uncertainty in ANS. 6: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. On exit: IFAIL = 0 unless the routine detects an error (see Section 6). 6. Error Indicators and Warnings Errors detected by the routine: IFAIL= 1 Indicates that fewer than four-points have been supplied to the routine. IFAIL= 2 Values of X are neither strictly increasing nor strictly decreasing. IFAIL= 3 Two points have the same X-value. No error is reported arising from the relative magnitudes of ANS and ER on return, due to the difficulty when the true answer is zero. 7. Accuracy No accuracy level is specified by the user before calling the routine but on return ABS(ER) is an approximation to, but not necessarily a bound for, |I-ANS|. If on exit IFAIL > 0, both ANS and ER are returned as zero. 8. Further Comments The time taken by the routine depends on the number of points supplied, n. In their paper, Gill and Miller [1] do not add the quantity ER to ANS before return. However, extensive tests have shown that a dramatic reduction in the error often results from such addition. In other cases, it does not make an improvement, but these tend to be cases of low accuracy in which the modified answer is not significantly inferior to the unmodified one. The user has the option of recovering the Gill-Miller answer by subtracting ER from ANS on return from the routine. 9. Example The example program evaluates the integral 1 / 4 | ----dx=(pi) / 2 0 1+x reading in the function values at 21 unequally-spaced points. The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXd01gbf}{NAG On-line Documentation: d01gbf} \beginscroll \begin{verbatim} D01GBF(3NAG) Foundation Library (12/10/92) D01GBF(3NAG) D01 -- Quadrature D01GBF D01GBF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose D01GBF returns an approximation to the integral of a function over a hyper-rectangular region, using a Monte Carlo method. An approximate relative error estimate is also returned. This routine is suitable for low accuracy work. 2. Specification SUBROUTINE D01GBF (NDIM, A, B, MINCLS, MAXCLS, FUNCTN, 1 EPS, ACC, LENWRK, WRKSTR, FINEST, IFAIL) INTEGER NDIM, MINCLS, MAXCLS, LENWRK, IFAIL DOUBLE PRECISION A(NDIM), B(NDIM), FUNCTN, EPS, ACC, WRKSTR 1 (LENWRK), FINEST EXTERNAL FUNCTN 3. Description D01GBF uses an adaptive Monte Carlo method based on the algorithm described by Lautrup [1]. It is implemented for integrals of the form: b b b 1 2 n / / / | | ... | f(x ,x ,...,x )dx ...dx dx . / / / 1 2 n n 2 1 a a a 1 2 n Upon entry, unless LENWRK has been set to the minimum value 10*NDIM, the routine subdivides the integration region into a number of equal volume subregions. Inside each subregion the integral and the variance are estimated by means of pseudo-random sampling. All contributions are added together to produce an estimate for the whole integral and total variance. The variance along each co-ordinate axis is determined and the routine uses this information to increase the density and change the widths of the sub-intervals along each axis, so as to reduce the total variance. The total number of subregions is then increased by a factor of two and the program recycles for another iteration. The program stops when a desired accuracy has been reached or too many integral evaluations are needed for the next cycle. 4. References [1] Lautrup B (1971) An Adaptive Multi-dimensional Integration Procedure. Proc. 2nd Coll. on Advanced Methods in Theoretical Physics, Marseille. 5. Parameters 1: NDIM -- INTEGER Input On entry: the number of dimensions of the integral, n. Constraint: NDIM >= 1. 2: A(NDIM) -- DOUBLE PRECISION array Input On entry: the lower limits of integration, a , for i i=1,2,...,n. 3: B(NDIM) -- DOUBLE PRECISION array Input On entry: the upper limits of integration, b , for i i=1,2,...,n. 4: MINCLS -- INTEGER Input/Output On entry: MINCLS must be set: either to the minimum number of integrand evaluations to be allowed, in which case MINCLS >= 0; or to a negative value. In this case the routine assumes that a previous call had been made with the same parameters NDIM, A and B and with either the same integrand (in which case D01GBF continues calculation) or a similar integrand (in which case D01GBF begins the calculation with the subdivision used in the last iteration of the previous call) . See also WRKSTR. On exit: MINCLS contains the number of integrand evaluations actually used by D01GBF. 5: MAXCLS -- INTEGER Input On entry: the maximum number of integrand evaluations to be allowed. In the continuation case this is the number of new integrand evaluations to be allowed. These counts do not include zero integrand values. Constraints: MAXCLS > MINCLS, MAXCLS >= 4*(NDIM+1). 6: FUNCTN -- DOUBLE PRECISION FUNCTION, supplied by the user. External Procedure FUNCTN must return the value of the integrand f at a given point. Its specification is: DOUBLE PRECISION FUNCTION FUNCTN (NDIM, X) INTEGER NDIM DOUBLE PRECISION X(NDIM) 1: NDIM -- INTEGER Input On entry: the number of dimensions of the integral, n. 2: X(NDIM) -- DOUBLE PRECISION array Input On entry: the co-ordinates of the point at which the integrand must be evaluated. FUNCTN must be declared as EXTERNAL in the (sub)program from which D01GBF is called. Parameters denoted as Input must not be changed by this procedure. 7: EPS -- DOUBLE PRECISION Input On entry: the relative accuracy required. Constraint: EPS >= 0.0. 8: ACC -- DOUBLE PRECISION Output On exit: the estimated relative accuracy of FINEST. 9: LENWRK -- INTEGER Input On entry: the dimension of the array WRKSTR as declared in the (sub)program from which D01GBF is called. For maximum efficiency, LENWRK should be about 1/NDIM 3*NDIM*(MAXCLS/4) +7*NDIM. If LENWRK is given the value 10*NDIM then the subroutine uses only one iteration of a crude Monte Carlo method with MAXCLS sample points. Constraint: LENWRK >= 10*NDIM. 10: WRKSTR(LENWRK) -- DOUBLE PRECISION array Input/Output On entry: if MINCLS<0, WRKSTR must be unchanged from the previous call of D01GBF - except that for a new integrand WRKSTR(LENWRK) must be set to 0.0. See also MINCLS. On exit: WRKSTR contains information about the current sub- interval structure which could be used in later calls of D01GBF. In particular, WRKSTR(j) gives the number of sub- intervals used along the jth co-ordinate axis. 11: FINEST -- DOUBLE PRECISION Output On exit: the best estimate obtained for the integral. 12: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. Users who are unfamiliar with this parameter should refer to the Essential Introduction for details. On exit: IFAIL = 0 unless the routine detects an error or gives a warning (see Section 6). For this routine, because the values of output parameters may be useful even if IFAIL /=0 on exit, users are recommended to set IFAIL to -1 before entry. It is then essential to test the value of IFAIL on exit. To suppress the output of an error message when soft failure occurs, set IFAIL to 1. 6. Error Indicators and Warnings Errors or warnings specified by the routine: IFAIL= 1 On entry NDIM < 1, or MINCLS >= MAXCLS, or LENWRK < 10*NDIM, or MAXCLS < 4*(NDIM+1), or EPS < 0.0. IFAIL= 2 MAXCLS was too small for D01GBF to obtain the required relative accuracy EPS. In this case D01GBF returns a value of FINEST with estimated relative error ACC, but ACC will be greater than EPS. This error exit may be taken before MAXCLS non-zero integrand evaluations have actually occurred, if the routine calculates that the current estimates could not be improved before MAXCLS was exceeded. 7. Accuracy A relative error estimate is output through the parameter ACC. The confidence factor is set so that the actual error should be less than ACC 90% of the time. If a user desires a higher confidence level then a smaller value of EPS should be used. 8. Further Comments The running time for D01GBF will usually be dominated by the time used to evaluate the integrand FUNCTN, so the maximum time that could be used is approximately proportional to MAXCLS. For some integrands, particularly those that are poorly behaved in a small part of the integration region, D01GBF may terminate with a value of ACC which is significantly smaller than the actual relative error. This should be suspected if the returned value of MINCLS is small relative to the expected difficulty of the integral. Where this occurs, D01GBF should be called again, but with a higher entry value of MINCLS (e.g. twice the returned value) and the results compared with those from the previous call. The exact values of FINEST and ACC on return will depend (within statistical limits) on the sequence of random numbers generated within D01GBF by calls to G05CAF. Separate runs will produce identical answers unless the part of the program executed prior to calling D01GBF also calls (directly or indirectly) routines from Chapter G05, and the series of such calls differs between runs. If desired, the user may ensure the identity or difference between runs of the results returned by D01GBF, by calling G05CBF or G05CCF respectively, immediately before calling D01GBF. 9. Example This example program calculates the integral 1 1 1 1 4x x exp(2x x ) / / / / 1 3 1 3 | | | | ---------------dx dx dx dx =0.575364. / / / / 2 1 2 3 4 0 0 0 0 (1+x +x ) 2 4 The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXd02}{NAG On-line Documentation: d02} \beginscroll \begin{verbatim} D02(3NAG) Foundation Library (12/10/92) D02(3NAG) D02 -- Ordinary Differential Equations Introduction -- D02 Chapter D02 Ordinary Differential Equations 1. Scope of the Chapter This chapter is concerned with the numerical solution of ordinary differential equations. There are two main types of problem, those in which all boundary conditions are specified at one point (initial-value problems), and those in which the boundary conditions are distributed between two or more points (boundary- value problems and eigenvalue problems). Routines are available for initial-value problems, two-point boundary-value problems and Sturm-Liouville eigenvalue problems. 2. Background to the Problems For most of the routines in this chapter a system of ordinary differential equations must be written in the form y' = f (x,y ,y ,...,y ), 1 1 1 2 n y' = f (x,y ,y ,...,y ), 2 2 1 2 n .. ... y' = f (x,y ,y ,...y ), n n 1 2 n that is the system must be given in first-order form. The n dependent variables (also, the solution) y ,y ,...,y are 1 2 n functions of the independent variable x, and the differential equations give expressions for the first derivatives y' =dy /dx i i in terms of x and y ,y ,...,y . For a system of n first-order 1 2 n equations, n associated boundary conditions are usually required to define the solution. A more general system may contain derivatives of higher order, but such systems can almost always be reduced to the first-order form by introducing new variables. For example, suppose we have the third-order equation 2 z'''+zz''+k(l-z' )=0. We write y =z, y =z', y =z'', and the third order equation may 1 2 3 then be written as the system of first-order equations y' = y 1 2 y' = y 2 3 2 y' = -y y -k(l-y ). 3 1 3 2 For this system n = 3 and we require 3 boundary conditions in order to define the solution. These conditions must specify values of the dependent variables at certain points. For example, we have an initial-value problem if the conditions are: y = 0 at x=0 1 y = 0 at x=0 2 y = 0.1 at x=0. 3 These conditions would enable us to integrate the equations numerically from the point x=0 to some specified end-point. We have a boundary-value problem if the conditions are: y = 0 at x=0 1 y = 0 at x=0 2 y = 1 at x=10. 2 These conditions would be sufficient to define a solution in the range 0<=x<=10, but the problem could not be solved by direct integration (see Section 2.2 ). More general boundary conditions are permitted in the boundary-value case. 2.1. Initial-value Problems To solve first-order systems, initial values of the dependent variables y , for i=1,2,...,n must be supplied at a given point, i a. Also a point, b, at which the values of the dependent variables are required, must be specified. The numerical solution is then obtained by a step-by-step calculation which approximates values of the variables y , for i=1,2,...,n at finite intervals i over the required range [a,b]. The routines in this chapter adjust the step length automatically to meet specified accuracy tolerances. Although the accuracy tests used are reliable over each step individually, in general an accuracy requirement cannot be guaranteed over a long range. For many problems there may be no serious accumulation of error, but for unstable systems small perturbations of the solution will often lead to rapid divergence of the calculated values from the true values. A simple check for stability is to carry out trial calculations with different tolerances; if the results differ appreciably the system is probably unstable. Over a short range, the difficulty may possibly be overcome by taking sufficiently small tolerances, but over a long range it may be better to try to reformulate the problem. A special class of initial-value problems are those for which the solutions contain rapidly decaying transient terms. Such problems are called stiff; an alternative way of describing them is to say that certain eigenvalues of the Jacobian matrix (ddf /ddy ) have i j large negative real parts when compared to others. These problems require special methods for efficient numerical solution; the methods designed for non-stiff problems when applied to stiff problems tend to be very slow, because they need small step lengths to avoid numerical instability. A full discussion is given in Hall and Watt [6] and a discussion of the methods for stiff problems is given in Berzins, Brankin and Gladwell [1]. 2.2. Boundary-value Problems A full discussion of the design of the methods and codes for boundary-value problems is given in Gladwell [4]. In general, a system of nonlinear differential equations with boundary conditions given at two or more points cannot be guaranteed to have a solution. The solution has to be determined iteratively (if it exists). Finite-difference equations are set up on a mesh of points and estimated values for the solution at the grid points are chosen. Using these estimated values as starting values a Newton iteration is used to solve the finite-difference equations. The accuracy of the solution is then improved by deferred corrections or the addition of points to the mesh or a combination of both. Good initial estimates of the solution may be required in some cases but results may be difficult to compute when the solution varies very rapidly over short ranges. A discussion is given in Chapters 9 and 11 of Gladwell and Sayers [5] and Chapter 4 of Childs et al [2]. 2.3. Eigenvalue Problems Sturm-Liouville problems of the form (p(x)y')'+q(x,(lambda))y=0 with appropriate boundary conditions given at two points, can be solved by a Scaled Pruefer method. In this method the differential equation is transformed to another which can be solved for a specified eigenvalue by a shooting method. A discussion is given in Chapter 11 of Gladwell and Sayers [5] and a complete description is given in Pryce [7]. 2.6. References [1] Berzins M, Brankin R W and Gladwell I (1987) Design of the Stiff Integrators in the NAG Library. Technical Report. TR14/87 NAG. [2] (1979) Codes for Boundary-value Problems in Ordinary Differential Equations. Lecture Notes in Computer Science. (ed Childs B, Scott M, Daniel J W, Denman E and Nelson P) 76 Springer-Verlag. [3] Gladwell I (1979) Initial Value Routines in the NAG Library. ACM Trans Math Softw. 5 386--400. [4] Gladwell I (1987) The NAG Library Boundary Value Codes. Numerical Analysis Report. 134 Manchester University. [5] (1980) Computational Techniques for Ordinary Differential Equations. (Gladwell I and Sayers D K) Academic Press. [6] Hall G and Watt J M (eds) (1976) Modern Numerical Methods for Ordinary Differential Equations. Clarendon Press. [7] Pryce J D (1986) Error Estimation for Phase-function Shooting Methods for Sturm-Liouville Problems. J. Num. Anal. 6 103--123. 3. Recommendations on Choice and Use of Routines There are no routines which deal directly with COMPLEX equations. These may however be transformed to larger systems of real equations of the required form. Split each equation into its real and imaginary parts and solve for the real and imaginary parts of each component of the solution. Whilst this process doubles the size of the system and may not always be appropriate it does make available for use the full range of routines provided presently. 3.1. Initial-value Problems For simple first-order problems with low accuracy requirements, that is problems on a short range of integration, with derivative functions f which are inexpensive to calculate and where only a i few correct figures are required, the best routines to use are likely to be the Runge-Kutta-Merson (RK) routines, D02BBF and D02BHF. For larger problems, over long ranges or with high accuracy requirements the variable-order, variable-step Adams routine D02CJF should usually be preferred. For stiff equations, that is those with rapidly decaying transient solutions, the Backward Differentiation Formula (BDF) variable-order, variable- step routine D02EJF should be used. There are four routines for initial-value problems, two of which use the Runge-Kutta-Merson method: D02BBF integrates a system of first order ordinary differential equations over a range with intermediate output and a choice of error control D02BHF integrates a system of first order ordinary differential equations with a choice of error control until a position is determined where a function of the solution is zero. one uses an Adams method: D02CJF combines the functionality of D02BBF and D02BHF and one uses a BDF method: D02EJF combines the functionality of D02BBF and D02BHF. 3.2. Boundary-value Problems D02GAF may be used for simple boundary-value problems with assigned boundary values. The user may find that convergence is difficult to achieve using D02GAF since only specifying the unknown boundary values and the position of the finite-difference mesh is permitted. In such cases the user may use D02RAF which permits specification of an initial estimate for the solution at all mesh points and allows the calculation to be influenced in other ways too. D02RAF is designed to solve a general nonlinear two-point boundary value problem with nonlinear boundary conditions. A routine, D02GBF, is also supplied specifically for the general linear two-point boundary-value problem written in a standard The user is advised to use interpolation routines from the E01 Chapter to obtain solution values at points not on the final mesh. 3.3. Eigenvalue Problems There is one general purpose routine for eigenvalue problems, D02KEF. It may be used to solve regular or singular second-order Sturm-Liouville problems on a finite or infinite range. Discontinous coefficient functions can be treated and eigenfunctions can be computed. D02 -- Ordinary Differential Equations Contents -- D02 Chapter D02 Ordinary Differential Equations D02BBF ODEs, IVP, Runge-Kutta-Merson method, over a range, intermediate output D02BHF ODEs, IVP, Runge-Kutta-Merson method, until function of solution is zero D02CJF ODEs, IVP, Adams method, until function of solution is zero, intermediate output D02EJF ODEs, stiff IVP, BDF method, until function of solution is zero, intermediate output D02GAF ODEs, boundary value problem, finite difference technique with deferred correction, simple nonlinear problem D02GBF ODEs, boundary value problem, finite difference technique with deferred correction, general linear problem D02KEF 2nd order Sturm-Liouville problem, regular/singular system, finite/infinite range, eigenvalue and eigenfunction, user-specified break-points D02RAF ODEs, general nonlinear boundary value problem, finite difference technique with deferred correction, continuation facility \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXd02bbf}{NAG On-line Documentation: d02bbf} \beginscroll \begin{verbatim} D02BBF(3NAG) D02BBF D02BBF(3NAG) D02 -- Ordinary Differential Equations D02BBF D02BBF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. Note for users via the AXIOM system: the interface to this routine has been enhanced for use with AXIOM and is slightly different to that offered in the standard version of the Foundation Library. 1. Purpose D02BBF integrates a system of first-order ordinary differential equations over an interval with suitable initial conditions, using a Runge-Kutta-Merson method, and returns the solution at points specified by the user. 2. Specification SUBROUTINE D02BBF (X, XEND, M, N, Y, TOL, IRELAB, RESULT, 1 FCN, OUTPUT, W, IFAIL) INTEGER M, N, IRELAB, IFAIL DOUBLE PRECISION X, XEND, Y(N), TOL, W(N,7), RESULT(M,N) EXTERNAL FCN, OUTPUT 3. Description The routine integrates a system of ordinary differential equations y' =f (x,y ,y ,...,y ) i=1,2,...,n i i 1 2 n from x = X to x = XEND using a Merson form of the Runge-Kutta method. The system is defined by a subroutine FCN supplied by the user, which evaluates f in terms of x and y ,y ,...,y , and the i 1 2 n values of y ,y ,...,y must be given at x = X. 1 2 n The solution is returned via the user-supplied routine OUTPUT at a set of points specified by the user. This solution is obtained by quintic Hermite interpolation on solution values produced by the Runge-Kutta method. The accuracy of the integration and, indirectly, the interpolation is controlled by the parameters TOL and IRELAB. For a description of Runge-Kutta methods and their practical implementation see Hall and Watt [1]. 4. References [1] Hall G and Watt J M (eds) (1976) Modern Numerical Methods for Ordinary Differential Equations. Clarendon Press. 5. Parameters 1: X -- DOUBLE PRECISION Input/Output On entry: X must be set to the initial value of the independent variable x. On exit: XEND, unless an error has occurred, when it contains the value of x at the error. 2: XEND -- DOUBLE PRECISION Input On entry: the final value of the independent variable. If XEND < X on entry, integration will proceed in a negative direction. 3: M -- INTEGER Input On entry: the first dimension of the array RESULT. This will usually be equal to the number of points at which the solution is required. Constraint: M > 0. 4: N -- INTEGER Input On entry: the number of differential equations. Constraint: N > 0. 5: Y(N) -- DOUBLE PRECISION array Input/Output On entry: the initial values of the solution y ,y ,...,y . 1 2 n On exit: the computed values of the solution at the final value of X. 6: TOL -- DOUBLE PRECISION Input/Output On entry: TOL must be set to a positive tolerance for controlling the error in the integration. D02BBF has been designed so that, for most problems, a reduction in TOL leads to an approximately proportional reduction in the error in the solution at XEND. The relation between changes in TOL and the error at intermediate output points is less clear, but for TOL small enough the error at intermediate output points should also be approximately proportional to TOL. However, the actual relation between TOL and the accuracy achieved cannot be guaranteed. The user is strongly recommended to call D02BBF with more than one value for TOL and to compare the results obtained to estimate their accuracy. In the absence of any prior knowledge, the user might compare the results obtained by -p -p-1 calling D02BBF with TOL=10.0 and TOL=10.0 if p correct decimal digits in the solution are required. Constraint: TOL > 0.0. On exit: normally unchanged. However if the range X to XEND is so short that a small change in TOL is unlikely to make any change in the computed solution then, on return, TOL has its sign changed. This should be treated as a warning that the computed solution is likely to be more accurate than would be produced by using the same value of TOL on a longer range of integration. 7: IRELAB -- INTEGER Input On entry: IRELAB determines the type of error control. At each step in the numerical solution an estimate of the local error, EST, is made. For the current step to be accepted the following condition must be satisfied: IRELAB = 0 EST=10.0<=TOL*max{1.0,|y |,|y |,...,|y |}; 1 2 n IRELAB = 1 EST <= TOL; IRELAB = 2 EST<=TOL*max{(epsilon),|y |,|y |,...,|y |}, where 1 2 n (epsilon) is machine precision. If the appropriate condition is not satisfied, the step size is reduced and the solution is recomputed on the current step. If the user wishes to measure the error in the computed solution in terms of the number of correct decimal places, then IRELAB should be given the value 1 on entry, whereas if the error requirement is in terms of the number of correct significant digits, then IRELAB should be given the value 2. Where there is no preference in the choice of error test IRELAB = 0 will result in a mixed error test. Constraint: 0 <= IRELAB <= 2. 8: RESULT(M,N) -- DOUBLE PRECISION array Output On exit: the computed values of the solution at the points given by OUTPUT. 9: FCN -- SUBROUTINE, supplied by the user. External Procedure FCN must evaluate the functions f (i.e., the derivatives i y' ) for given values of its arguments x,y ,...,y . i 1 n Its specification is: SUBROUTINE FCN (X, Y, F) DOUBLE PRECISION X, Y(n), F(n) where n is the actual value of N in the call of D02BBF. 1: X -- DOUBLE PRECISION Input On entry: the value of the argument x. 2: Y(*) -- DOUBLE PRECISION array Input On entry: the value of the argument y , for i i=1,2,...,n. 3: F(*) -- DOUBLE PRECISION array Output On exit: the value of f , for i=1,2,...,n. i FCN must be declared as EXTERNAL in the (sub)program from which D02BBF is called. Parameters denoted as Input must not be changed by this procedure. 10: OUTPUT -- SUBROUTINE, supplied by the user. External Procedure OUTPUT allows the user to have access to intermediate values of the computed solution at successive points specified by the user. These solution values may be returned to the user via the array RESULT if desired (this is a non-standard feature added for use with the AXIOM system). OUTPUT is initially called by D02BBF with XSOL = X (the initial value of x). The user must reset XSOL to the next point where OUTPUT is to be called, and so on at each call to OUTPUT. If, after a call to OUTPUT, the reset point XSOL is beyond XEND, D02BBF will integrate to XEND with no further calls to OUTPUT; if a call to OUTPUT is required at the point XSOL = XEND, then XSOL must be given precisely the value XEND. Its specification is: SUBROUTINE OUTPUT(XSOL,Y,COUNT,M,N,RESULT) DOUBLE PRECISION Y(N),RESULT(M,N),XSOL INTEGER M,N,COUNT 1: XSOL -- DOUBLE PRECISION Input/Output On entry: the current value of the independent variable x. On exit: the next value of x at which OUTPUT is to be called. 2: Y(N) -- DOUBLE PRECISION array Input On entry: the computed solution at the point XSOL. 3: COUNT -- INTEGER Input/Output On entry: Zero if OUTPUT has not been called before, or the previous value of COUNT. On exit: A new value of COUNT: this can be used to keep track of the number of times OUTPUT has been called. 4: M -- INTEGER Input On entry: The first dimension of RESULT. 5: N -- INTEGER Input On entry: The dimension of Y. 6: RESULT(M,N) -- DOUBLE PRECISION array Input/Output On entry: the previous contents of RESULT. On exit: RESULT may be used to return the values of the intermediate solutions to the user. OUTPUT must be declared as EXTERNAL in the (sub)program from which D02BBF is called. Parameters denoted as Input must not be changed by this procedure. 11: W(N,7) -- DOUBLE PRECISION array Workspace 12: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. On exit: IFAIL = 0 unless the routine detects an error (see Section 6). 6. Error Indicators and Warnings Errors detected by the routine: IFAIL= 1 On entry TOL <= 0, or N <= 0, or IRELAB /= 0, 1 or 2. IFAIL= 2 With the given value of TOL, no further progress can be made across the integration range from the current point x = X, or the dependence of the error on TOL would be lost if further progress across the integration range were attempted (see Section 8 for a discussion of this error exit). The components Y(1),Y(2),...,Y(n) contain the computed values of the solution at the current point x = X. IFAIL= 3 TOL is too small for the routine to take an initial step (see Section 8). X and Y(1),Y(2),...,Y(n) retain their initial values. IFAIL= 4 X = XEND and XSOL /= X after the initial call to OUTPUT. IFAIL= 5 A value of XSOL returned by OUTPUT lies behind the previous value of XSOL in the direction of integration. IFAIL= 6 A serious error has occurred in an internal call to D02PAF(*). Check all subroutine calls and array dimensions. Seek expert help. IFAIL= 7 A serious error has occurred in an internal call to D02XAF(*). Check all subroutine calls and array dimensions. Seek expert help. 7. Accuracy The accuracy depends on TOL, on the mathematical properties of the differential system, on the length of the range of integration and on the method. It can be controlled by varying TOL but the approximate proportionality of the error to TOL holds only for a restricted range of values of TOL. For TOL too large, the underlying theory may break down and the result of varying TOL may be unpredictable. For TOL too small, rounding errors may affect the solution significantly and an error exit with IFAIL = 2 or IFAIL = 3 is possible. At the intermediate output points the same remarks apply. For large values of TOL it is possible that the errors at some intermediate output points may be much larger than at XEND. In any case, it must not be expected that the error will have the same size at all output points. At any point, it is a combination of the errors arising from the integration of the differential equation and the interpolation. The effect of combining these errors will vary, though in most cases the integration error will dominate. The user who requires a more reliable estimate of the accuracy achieved than can be obtained by varying TOL, is recommended to call D02BDF(*) where both the solution and a global error estimate are computed. 8. Further Comments The time taken by the routine depends on the complexity and mathematical properties of the system of differential equations defined by FCN, on the range, the tolerance and the number of calls to OUTPUT. There is also an overhead of the form a+b*n where a and b are machine-dependent computing times. If the routine fails with IFAIL = 3, then it can be called again with a larger value of TOL (if this has not already been tried). If the accuracy requested is really needed and cannot be obtained with this routine, the system may be very stiff (see below) or so badly scaled that it cannot be solved to the required accuracy. If the routine fails with IFAIL = 2, it is probable that it has been called with a value of TOL which is so small that the solution cannot be obtained on the range X to XEND. This can happen for well-behaved systems and very small values of TOL. The user should, however, consider whether there is a more fundamental difficulty. For example: (a) in the region of a singularity (infinite value) of the solution, the routine will usually stop with IFAIL = 2, unless overflow occurs first. If overflow occurs using D02BBF, D02PAF(*) can be used instead to trap the increasing solution before overflow occurs. In any case, numerical integration cannot be continued through a singularity, and analytic treatment should be considered; (b) for 'stiff' equations, where the solution contains rapidly decaying components, the routine will use very small steps in x (internally to D02BBF) to preserve stability. This will usually exhibit itself by making the computing time excessively long, or occasionally by an exit with IFAIL = 2. Merson's method is not efficient in such cases, and the user should try using D02EBF(*) (Backward Differentiation Formula). To determine whether a problem is stiff, D02BDF(*) may be used. For well-behaved systems with no difficulties such as stiffness or singularities, the Merson method should work well for low accuracy calculations (three or four figures). For higher accuracy calculations or where FCN is costly to evaluate, Merson's method may not be appropriate and a computationally less expensive method may be D02CBF(*) which uses an Adams method. Users with problems for which D02BBF is not sufficiently general should consider using D02PAF(*) with D02XAF(*). D02PAF(*) is a more general Merson routine with many facilities including more general error control options and several criteria for interrupting the calculations. D02XAF(*) interpolates on values produced by D02PAF(*). 9. Example To integrate the following equations (for a projectile) y'=tan(phi) -0.032tan(phi) 0.02v v'= --------------- -------- v cos(phi) -0.032 (phi)'= ------ 2 v over an interval X = 0.0 to XEND = 8.0, starting with values y=0.0, v=0.5 and (phi)=(pi)/5 and printing the solution at steps of 1.0. We write y=Y(1), v=Y(2) and (phi)=Y(3), and we set TOL=1.0E-4 and TOL=1.0E-5 in turn so that we may compare the solutions. The value of (pi) is obtained by using X01AAF(*). Note the careful construction of routine OUT to ensure that the value of XEND is printed. The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXd02bhf}{NAG On-line Documentation: d02bhf} \beginscroll \begin{verbatim} D02BHF(3NAG) Foundation Library (12/10/92) D02BHF(3NAG) D02 -- Ordinary Differential Equations D02BHF D02BHF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose D02BHF integrates a system of first-order ordinary differential equations over an interval with suitable initial conditions, using a Runge-Kutta-Merson method, until a user-specified function of the solution is zero. 2. Specification SUBROUTINE D02BHF (X, XEND, N, Y, TOL, IRELAB, HMAX, FCN, 1 G, W, IFAIL) INTEGER N, IRELAB, IFAIL DOUBLE PRECISION X, XEND, Y(N), TOL, HMAX, G, W(N,7) EXTERNAL FCN, G 3. Description The routine advances the solution of a system of ordinary differential equations y' =f (x,y ,y ,...,y ), i=1,2,...,n, i i 1 2 n from x = X towards x = XEND using a Merson form of the Runge- Kutta method. The system is defined by a subroutine FCN supplied by the user, which evaluates f in terms of x and y ,y ,...,y i 1 2 n (see Section 5), and the values of y ,y ,...,y must be given at 1 2 n x = X. As the integration proceeds, a check is made on the function g(x,y) specified by the user, to determine an interval where it changes sign. The position of this sign change is then determined accurately by interpolating for the solution and its derivative. It is assumed that g(x,y) is a continuous function of the variables, so that a solution of g(x,y) = 0 can be determined by searching for a change in sign in g(x,y). The accuracy of the integration and, indirectly, of the determination of the position where g(x,y) = 0, is controlled by the parameter TOL. For a description of Runge-Kutta methods and their practical implementation see Hall and Watt [1]. 4. References [1] Hall G and Watt J M (eds) (1976) Modern Numerical Methods for Ordinary Differential Equations. Clarendon Press. 5. Parameters 1: X -- DOUBLE PRECISION Input/Output On entry: X must be set to the initial value of the independent variable x. On exit: the point where g(x,y) = 0. 0 unless an error has occurred, when it contains the value of x at the error. In particular, if g(x,y)/=0.0 anywhere on the range X to XEND, it will contain XEND on exit. 2: XEND -- DOUBLE PRECISION Input On entry: the final value of the independent variable x. If XEND < X on entry, integration proceeds in a negative direction. 3: N -- INTEGER Input On entry: the number of differential equations, n. Constraint: N > 0. 4: Y(N) -- DOUBLE PRECISION array Input/Output On entry: the initial values of the solution y ,y ,...,y . 1 2 n On exit: the computed values of the solution at the final point x = X. 5: TOL -- DOUBLE PRECISION Input/Output On entry: TOL must be set to a positive tolerance for controlling the error in the integration and in the determination of the position where g(x,y) = 0.0. D02BHF has been designed so that, for most problems, a reduction in TOL leads to an approximately proportional reduction in the error in the solution obtained in the integration. The relation between changes in TOL and the error in the determination of the position where g(x,y) = 0. 0 is less clear, but for TOL small enough the error should be approximately proportional to TOL. However, the actual relation between TOL and the accuracy cannot be guaranteed. The user is strongly recommended to call D02BHF with more than one value for TOL and to compare the results obtained to estimate their accuracy. In the absence of any prior knowledge the user might compare results obtained by calling -p -p-1 D02BHF with TOL=10.0 and TOL=10.0 if p correct decimal digits in the solution are required. Constraint: TOL > 0.0. On exit: normally unchanged. However if the range from x = X to the position where g(x,y) = 0.0 (or to the final value of x if an error occurs) is so short that a small change in TOL is unlikely to make any change in the computed solution, then TOL is returned with its sign changed. To check results returned with TOL < 0.0, D02BHF should be called again with a positive value of TOL whose magnitude is considerably smaller than that of the previous call. 6: IRELAB -- INTEGER Input On entry: IRELAB determines the type of error control. At each step in the numerical solution an estimate of the local error, EST, is made. For the current step to be accepted the following condition must be satisfied: IRELAB = 0 EST<=TOL*max{1.0,|y |,|y |,...,|y |}; 1 2 n IRELAB = 1 EST <= TOL; IRELAB = 2 EST<=TOL*max{(epsilon),|y |,|y |,...,|y |}, 1 2 n where (epsilon) is machine precision. If the appropriate condition is not satisfied, the step size is reduced and the solution recomputed on the current step. If the user wishes to measure the error in the computed solution in terms of the number of correct decimal places, then IRELAB should be given the value 1 on entry, whereas if the error requirement is in terms of the number of correct significant digits, then IRELAB should be given the value 2. Where there is no preference in the choice of error test, IRELAB = 0 will result in a mixed error test. It should be borne in mind that the computed solution will be used in evaluating g(x,y). Constraint: 0 <= IRELAB <= 2. 7: HMAX -- DOUBLE PRECISION Input On entry: if HMAX = 0.0, no special action is taken. If HMAX /= 0.0, a check is made for a change in sign of g(x,y) at steps not greater than |HMAX|. This facility should be used if there is any chance of 'missing' the change in sign by checking too infrequently. For example, if two changes of sign of g(x,y) are expected within a distance h, say, of each other, then a suitable value for HMAX might be HMAX = h/2. If only one change of sign in g(x,y) is expected on the range X to XEND, then the choice HMAX = 0.0 is most appropriate. 8: FCN -- SUBROUTINE, supplied by the user. External Procedure FCN must evaluate the functions f (i.e., the derivatives i y' ) for given values of its arguments x,y ,...,y . i 1 n Its specification is: SUBROUTINE FCN (X, Y, F) DOUBLE PRECISION X, Y(n), F(n) 1: X -- DOUBLE PRECISION Input On entry: the value of the argument x. 2: Y(*) -- DOUBLE PRECISION array Input On entry: the value of the argument y , for i i=1,2,...,n. 3: F(*) -- DOUBLE PRECISION array Output On exit: the value of f , for i=1,2,...,n. i FCN must be declared as EXTERNAL in the (sub)program from which D02BHF is called. Parameters denoted as Input must not be changed by this procedure. 9: G -- DOUBLE PRECISION FUNCTION, supplied by the user. External Procedure G must evaluate the function g(x,y) at a specified point. Its specification is: DOUBLE PRECISION FUNCTION G (X, Y) DOUBLE PRECISION X, Y(n) where n is the actual value of N in the call of D02BHF. 1: X -- DOUBLE PRECISION Input On entry: the value of the independent variable x. 2: Y(*) -- DOUBLE PRECISION array Input On entry: the value of y , for i=1,2,...,n. i G must be declared as EXTERNAL in the (sub)program from which D02BHF is called. Parameters denoted as Input must not be changed by this procedure. 10: W(N,7) -- DOUBLE PRECISION array Workspace 11: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. On exit: IFAIL = 0 unless the routine detects an error (see Section 6). 6. Error Indicators and Warnings Errors detected by the routine: IFAIL= 1 On entry TOL <= 0.0, or N <= 0, or IRELAB /= 0, 1 or 2. IFAIL= 2 With the given value of TOL, no further progress can be made across the integration range from the current point x = X, or dependence of the error on TOL would be lost if further progress across the integration range were attempted (see Section 8 for a discussion of this error exit). The components Y(1),Y(2),...,Y(n) contain the computed values of the solution at the current point x = X. No point at which g(x,y) changes sign has been located up to the point x = X. IFAIL= 3 TOL is too small for the routine to take an initial step (see Section 8). X and Y(1),Y(2),...,Y(n) retain their initial values. IFAIL= 4 At no point in the range X to XEND did the function g(x,y) change sign. It is assumed that g(x,y) = 0.0 has no solution. IFAIL= 5 A serious error has occurred in an internal call to C05AZF(*). Check all subroutine calls and array dimensions. Seek expert help. IFAIL= 6 A serious error has occurred in an internal call to D02PAF(*). Check all subroutine calls and array dimensions. Seek expert help. IFAIL= 7 A serious error has occurred in an internal call to D02XAF(*). Check all subroutine calls and array dimensions. Seek expert help. 7. Accuracy The accuracy depends on TOL, on the mathematical properties of the differential system, on the position where g(x,y) = 0.0 and on the method. It can be controlled by varying TOL but the approximate proportionality of the error to TOL holds only for a restricted range of values of TOL. For TOL too large, the underlying theory may break down and the result of varying TOL may be unpredictable. For TOL too small, rounding error may affect the solution significantly and an error exit with IFAIL = 2 or IFAIL = 3 is possible. The accuracy may also be restricted by the properties of g(x,y). The user should try to code G without introducing any unnecessary cancellation errors. 8. Further Comments The time taken by the routine depends on the complexity and mathematical properties of the system of differential equations defined by FCN, the complexity of G, on the range, the position of the solution and the tolerance. There is also an overhead of the form a+b*n where a and b are machine-dependent computing times. For some problems it is possible that D02BHF will return IFAIL = 4 because of inaccuracy of the computed values Y, leading to inaccuracy in the computed values of g(x,y) used in the search for the solution of g(x,y) = 0.0. This difficulty can be overcome by reducing TOL sufficiently, and if necessary, by choosing HMAX sufficiently small. If possible, the user should choose XEND well beyond the expected point where g(x,y) = 0.0; for example make |XEND-X| about 50 larger than the expected range. As a simple check, if, with XEND fixed, a change in TOL does not lead to a significant change in Y at XEND, then inaccuracy is not a likely source of error. If the routine fails with IFAIL = 3, then it could be called again with a larger value of TOL if this has not already been tried. If the accuracy requested is really needed and cannot be obtained with this routine, the system may be very stiff (see below) or so badly scaled that it cannot be solved to the required accuracy. If the routine fails with IFAIL = 2, it is likely that it has been called with a value of TOL which is so small that a solution cannot be obtained on the range X to XEND. This can happen for well-behaved systems and very small values of TOL. The user should, however, consider whether there is a more fundamental difficulty. For example: (a) in the region of a singularity (infinite value) of the solution, the routine will usually stop with IFAIL = 2, unless overflow occurs first. If overflow occurs using D02BHF, D02PAF(*) can be used instead to trap the increasing solution, before overflow occurs. In any case, numerical integration cannot be continued through a singularity, and analytical treatment should be considered; (b) for 'stiff' equations, where the solution contains rapidly decaying components, the routine will compute in very small steps in x (internally to D02BHF) to preserve stability. This will usually exhibit itself by making the computing time excessively long, or occasionally by an exit with IFAIL = 2. Merson's method is not efficient in such cases, and the user should try D02EHF(*) which uses a Backward Differentiation Formula method. To determine whether a problem is stiff, D02BDF(*) may be used. For well-behaved systems with no difficulties such as stiffness or singularities, the Merson method should work well for low accuracy calculations (three or four figures). For high accuracy calculations or where FCN is costly to evaluate, Merson's method may not be appropriate and a computationally less expensive method may be D02CHF(*) which uses an Adams method. For problems for which D02BHF is not sufficiently general, the user should consider D02PAF(*). D02PAF(*) is a more general Merson routine with many facilities including more general error control options and several criteria for interrupting the calculations. D02PAF(*) can be combined with the rootfinder C05AZF(*) and the interpolation routine D02XAF(*) to solve equations involving y ,y ,...,y and their derivatives. 1 2 n D02BHF can also be used to solve an equation involving x, y ,y ,...,y and the derivatives of y ,y ,...,y . For example in 1 2 n 1 2 n Section 9, D02BHF is used to find a value of X > 0.0 where Y(1) = 0.0. It could instead be used to find a turning-point of y by 1 replacing the function g(x,y) in the program by: DOUBLE PRECISION FUNCTION G(X,Y) DOUBLE PRECISION X,Y(3),F(3) CALL FCN(X,Y,F) G = F(1) RETURN END This routine is only intended to locate the first zero of g(x,y). If later zeros are required, users are strongly advised to construct their own more general root finding routines as discussed above. 9. Example To find the value X > 0.0 at which y=0.0, where y, v, (phi) are defined by y'=tan(phi) -0.032tan(phi) 0.02v v'= --------------- -------- v cos(phi) 0.032 (phi)'=- ----- 2 v and where at X = 0.0 we are given y=0.5, v=0.5 and (phi)=(pi)/5. We write y=Y(1), v=Y(2) and (phi)=Y(3) and we set TOL=1.0E-4 and TOL=1.0E-5 in turn so that we can compare the solutions. We expect the solution X ~= 7.3 and so we set XEND = 10.0 to avoid determining the solution of y=0.0 too near the end of the range of integration. The value of (pi) is obtained by using X01AAF(*). The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXd02cjf}{NAG On-line Documentation: d02cjf} \beginscroll \begin{verbatim} D02CJF(3NAG) D02CJF D02CJF(3NAG) D02 -- Ordinary Differential Equations D02CJF D02CJF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. Note for users via the AXIOM system: the interface to this routine has been enhanced for use with AXIOM and is slightly different to that offered in the standard version of the Foundation Library. 1. Purpose D02CJF integrates a system of first-order ordinary differential equations over a range with suitable initial conditions, using a variable-order, variable-step Adams method until a user-specified function, if supplied, of the solution is zero, and returns the solution at points specified by the user, if desired. 2. Specification SUBROUTINE D02CJF (X, XEND, M, N, Y, FCN, TOL, RELABS, 1 RESULT, OUTPUT, G, W, IFAIL) INTEGER M, N, IFAIL DOUBLE PRECISION X, XEND, Y(N), TOL, G, W(28+21*N), RESULT(M,N) CHARACTER*1 RELABS EXTERNAL FCN, OUTPUT, G 3. Description The routine advances the solution of a system of ordinary differential equations y' =f (x,y ,y ,...,y ), i=1,2,...,n, i i 1 2 n from x = X to x = XEND using a variable-order, variable-step Adams method. The system is defined by a subroutine FCN supplied by the user, which evaluates f in terms of x and y ,y ,...,y . i 1 2 n The initial values of y ,y ,...,y must be given at x = X. 1 2 n The solution is returned via the user-supplied routine OUTPUT at points specified by the user, if desired: this solution is 1 obtained by C interpolation on solution values produced by the method. As the integration proceeds a check can be made on the user-specified function g(x,y) to determine an interval where it changes sign. The position of this sign change is then determined 1 accurately by C interpolation to the solution. It is assumed that g(x,y) is a continuous function of the variables, so that a solution of g(x,y)=0.0 can be determined by searching for a change in sign in g(x,y). The accuracy of the integration, the interpolation and, indirectly, of the determination of the position where g(x,y)=0.0, is controlled by the parameters TOL and RELABS. For a description of Adams methods and their practical implementation see Hall and Watt [1]. 4. References [1] Hall G and Watt J M (eds) (1976) Modern Numerical Methods for Ordinary Differential Equations. Clarendon Press. 5. Parameters 1: X -- DOUBLE PRECISION Input/Output On entry: the initial value of the independent variable x. Constraint: X /= XEND. On exit: if g is supplied by the user, it contains the point where g(x,y)=0.0, unless g(x,y)/=0.0 anywhere on the range X to XEND, in which case, X will contain XEND. If g is not supplied by the user it contains XEND, unless an error has occurred, when it contains the value of x at the error. 2: XEND -- DOUBLE PRECISION Input On entry: the final value of the independent variable. If XEND < X, integration proceeds in the negative direction. Constraint: XEND /= X. 3: M -- INTEGER Input On entry: the first dimension of the array RESULT. This will usually be equal to the number of points at which the solution is required. Constraint: M > 0. 4: N -- INTEGER Input On entry: the number of differential equations. Constraint: N >= 1. 5: Y(N) -- DOUBLE PRECISION array Input/Output On entry: the initial values of the solution y ,y ,...,y 1 2 n at x = X. On exit: the computed values of the solution at the final point x = X. 6: FCN -- SUBROUTINE, supplied by the user. External Procedure FCN must evaluate the functions f (i.e., the derivatives i y' ) for given values of their arguments x,y ,y ,...,y . i 1 2 n Its specification is: SUBROUTINE FCN (X, Y, F) DOUBLE PRECISION X, Y(n), F(n) where n is the actual value of N in the call of D02CJF. 1: X -- DOUBLE PRECISION Input On entry: the value of the independent variable x. 2: Y(*) -- DOUBLE PRECISION array Input On entry: the value of the variable y , for i i=1,2,...,n. 3: F(*) -- DOUBLE PRECISION array Output On exit: the value of f , for i=1,2,...,n. i FCN must be declared as EXTERNAL in the (sub)program from which D02CJF is called. Parameters denoted as Input must not be changed by this procedure. 7: TOL -- DOUBLE PRECISION Input On entry: a positive tolerance for controlling the error in the integration. Hence TOL affects the determination of the position where g(x,y)=0.0, if g is supplied. D02CJF has been designed so that, for most problems, a reduction in TOL leads to an approximately proportional reduction in the error in the solution. However, the actual relation between TOL and the accuracy achieved cannot be guaranteed. The user is strongly recommended to call D02CJF with more than one value for TOL and to compare the results obtained to estimate their accuracy. In the absence of any prior knowledge, the user might compare the results obtained -p -p-1 by calling D02CJF with TOL=10.0 and TOL=10.0 where p correct decimal digits are required in the solution. Constraint: TOL > 0.0. 8: RELABS -- CHARACTER*1 Input On entry: the type of error control. At each step in the numerical solution an estimate of the local error, EST, is made. For the current step to be accepted the following condition must be satisfied: ______________________________ / n / -- 2 EST= / > (e /((tau) *|y |+(tau) )) <=1.0 / -- i r i a \/ i=1 where (tau) and (tau) are defined by r a RELABS (tau) (tau) r a 'M' TOL TOL 'A' 0.0 TOL 'R' TOL (epsilon) 'D' TOL TOL where (epsilon) is a small machine-dependent number and e i is an estimate of the local error at y , computed i internally. If the appropriate condition is not satisfied, the step size is reduced and the solution is recomputed on the current step. If the user wishes to measure the error in the computed solution in terms of the number of correct decimal places, then RELABS should be set to 'A' on entry, whereas if the error requirement is in terms of the number of correct significant digits, then RELABS should be set to 'R'. If the user prefers a mixed error test, then RELABS should be set to 'M', otherwise if the user has no preference, RELABS should be set to the default 'D'. Note that in this case 'D' is taken to be 'M'. Constraint: RELABS = 'M', 'A', 'R', 'D'. 9: RESULT(M,N) -- DOUBLE PRECISION array Output On exit: the computed values of the solution at the points given by OUTPUT. 10: OUTPUT -- SUBROUTINE, supplied by the user. External Procedure OUTPUT allows the user to have access to intermediate values of the computed solution at successive points specified by the user. These solution values may be returned to the user via the array RESULT if desired (this is a non-standard feature added for use with the AXIOM system). OUTPUT is initially called by D02CJF with XSOL = X (the initial value of x). The user must reset XSOL to the next point where OUTPUT is to be called, and so on at each call to OUTPUT. If, after a call to OUTPUT, the reset point XSOL is beyond XEND, D02CJF will integrate to XEND with no further calls to OUTPUT; if a call to OUTPUT is required at the point XSOL = XEND, then XSOL must be given precisely the value XEND. Its specification is: SUBROUTINE OUTPUT(XSOL,Y,COUNT,M,N,RESULT) DOUBLE PRECISION Y(N),RESULT(M,N),XSOL INTEGER M,N,COUNT 1: XSOL -- DOUBLE PRECISION Input/Output On entry: the current value of the independent variable x. On exit: the next value of x at which OUTPUT is to be called. 2: Y(N) -- DOUBLE PRECISION array Input On entry: the computed solution at the point XSOL. 3: COUNT -- INTEGER Input/Output On entry: Zero if OUTPUT has not been called before, or the previous value of COUNT. On exit: A new value of COUNT: this can be used to keep track of the number of times OUTPUT has been called. 4: M -- INTEGER Input On entry: The first dimension of RESULT. 5: N -- INTEGER Input On entry: The dimension of Y. 6: RESULT(M,N) -- DOUBLE PRECISION array Input/Output On entry: the previous contents of RESULT. On exit: RESULT may be used to return the values of the intermediate solutions to the user. OUTPUT must be declared as EXTERNAL in the (sub)program from which D02CJF is called. Parameters denoted as Input must not be changed by this procedure. 11: G -- DOUBLE PRECISION FUNCTION, supplied by the user. External Procedure G must evaluate the function g(x,y) for specified values x,y . It specifies the function g for which the first position x where g(x,y)=0 is to be found. If the user does not require the root finding option, the actual argument G mustbe the dummy routine D02CJW. (D02CJW is included in the NAG Foundation Library and so need not be supplied by the user). Its specification is: DOUBLE PRECISION FUNCTION G (X, Y) DOUBLE PRECISION X, Y(n) where n is the actual value of N in the call of D02CJF. 1: X -- DOUBLE PRECISION Input On entry: the value of the independent variable x. 2: Y(*) -- DOUBLE PRECISION array Input On entry: the value of the variable y , for i i=1,2,...,n. G must be declared as EXTERNAL in the (sub)program from which D02CJF is called. Parameters denoted as Input must not be changed by this procedure. 12: W(28+21*N) -- DOUBLE PRECISION array Workspace 13: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. On exit: IFAIL = 0 unless the routine detects an error (see Section 6). 6. Error Indicators and Warnings Errors detected by the routine: If on entry IFAIL = 0 or -1, explanatory error messages are output on the current error message unit (as defined by X04AAF). IFAIL= 1 On entry TOL <= 0.0, or N <= 0, or RELABS /= 'M', 'A', 'R' or 'D'. or X = XEND. IFAIL= 2 With the given value of TOL, no further progress can be made across the integration range from the current point x = X. (See Section 8 for a discussion of this error exit.) The components Y(1),Y(2),...,Y(N) contain the computed values of the solution at the current point x = X. If the user has supplied g, then no point at which g(x,y) changes sign has been located up to the point x = X. IFAIL= 3 TOL is too small for D02CJF to take an initial step. X and Y (1),Y(2),...,Y(N) retain their initial values. IFAIL= 4 XSOL has not been reset or XSOL lies behind X in the direction of integration, after the initial call to OUTPUT, if the OUTPUT option was selected. IFAIL= 5 A value of XSOL returned by OUTPUT has not been reset or lies behind the last value of XSOL in the direction of integration, if the OUTPUT option was selected. IFAIL= 6 At no point in the range X to XEND did the function g(x,y) change sign, if g was supplied. It is assumed that g(x,y)=0 has no solution. IFAIL= 7 A serious error has occurred in an internal call. Check all subroutine calls and array sizes. Seek expert help. 7. Accuracy The accuracy of the computation of the solution vector Y may be controlled by varying the local error tolerance TOL. In general, a decrease in local error tolerance should lead to an increase in accuracy. Users are advised to choose RELABS = 'M' unless they have a good reason for a different choice. If the problem is a root-finding one, then the accuracy of the root determined will depend on the properties of g(x,y). The user should try to code G without introducing any unnecessary cancellation errors. 8. Further Comments If more than one root is required then D02QFF(*) should be used. If the routine fails with IFAIL = 3, then it can be called again with a larger value of TOL if this has not already been tried. If the accuracy requested is really needed and cannot be obtained with this routine, the system may be very stiff (see below) or so badly scaled that it cannot be solved to the required accuracy. If the routine fails with IFAIL = 2, it is probable that it has been called with a value of TOL which is so small that a solution cannot be obtained on the range X to XEND. This can happen for well-behaved systems and very small values of TOL. The user should, however, consider whether there is a more fundamental difficulty. For example: (a) in the region of a singularity (infinite value) of the solution, the routine will usually stop with IFAIL = 2, unless overflow occurs first. Numerical integration cannot be continued through a singularity, and analytic treatment should be considered; (b) for 'stiff' equations where the solution contains rapidly decaying components, the routine will use very small steps in x (internally to D02CJF) to preserve stability. This will exhibit itself by making the computing time excessively long, or occasionally by an exit with IFAIL = 2. Adams methods are not efficient in such cases, and the user should try D02EJF. 9. Example We illustrate the solution of four different problems. In each case the differential system (for a projectile) is y'=tan(phi) -0.032tan(phi) 0.02v v'= --------------- -------- v cos(phi) -0.032 (phi)'= ------ 2 v over an interval X = 0.0 to XEND = 10.0 starting with values y=0.5, v=0.5 and (phi)=(pi)/5. We solve each of the following problems with local error tolerances 1.0E-4 and 1.0E-5. (i) To integrate to x=10.0 producing output at intervals of 2.0 until a point is encountered where y=0.0. (ii) As (i) but with no intermediate output. (iii) As (i) but with no termination on a root-finding condition. (iv) As (i) but with no intermediate output and no root-finding termination condition. The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXd02ejf}{NAG On-line Documentation: d02ejf} \beginscroll \begin{verbatim} D02EJF(3NAG) D02EJF D02EJF(3NAG) D02 -- Ordinary Differential Equations D02EJF D02EJF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. Note for users via the AXIOM system: the interface to this routine has been enhanced for use with AXIOM and is slightly different to that offered in the standard version of the Foundation Library. 1. Purpose D02EJF integrates a stiff system of first-order ordinary differential equations over an interval with suitable initial conditions, using a variable-order, variable-step method implementing the Backward Differentiation Formulae (BDF), until a user-specified function, if supplied, of the solution is zero, and returns the solution at points specified by the user, if desired. 2. Specification SUBROUTINE D02EJF (X, XEND, M, N, Y, FCN, PEDERV, TOL, 1 RELABS, OUTPUT, G, W, IW, RESULT, IFAIL) INTEGER M, N, IW, IFAIL DOUBLE PRECISION X, XEND, Y(N), TOL, G, W(IW), RESULT(M,N) CHARACTER*1 RELABS EXTERNAL FCN, PEDERV, OUTPUT, G 3. Description The routine advances the solution of a system of ordinary differential equations y' =f (x,y ,y ,...,y ), i=1,2,...,n, i i 1 2 n from x = X to x = XEND using a variable-order, variable-step method implementing the BDF. The system is defined by a subroutine FCN supplied by the user, which evaluates f in terms i of x and y ,y ,...,y (see Section 5). The initial values of 1 2 n y ,y ,...,y must be given at x = X. 1 2 n The solution is returned via the user-supplied routine OUTPUT at points specified by the user, if desired: this solution is 1 obtained by C interpolation on solution values produced by the method. As the integration proceeds a check can be made on the user-specified function g(x,y) to determine an interval where it changes sign. The position of this sign change is then determined 1 accurately by C interpolation to the solution. It is assumed that g(x,y) is a continuous function of the variables, so that a solution of g(x,y) = 0.0 can be determined by searching for a change in sign in g(x,y). The accuracy of the integration, the interpolation and, indirectly, of the determination of the position where g(x,y) = 0.0, is controlled by the parameters TOL and RELABS. The Jacobian of the system y'=f(x,y) may be supplied in routine PEDERV, if it is available. For a description of BDF and their practical implementation see Hall and Watt [1]. 4. References [1] Hall G and Watt J M (eds) (1976) Modern Numerical Methods for Ordinary Differential Equations. Clarendon Press. 5. Parameters 1: X -- DOUBLE PRECISION Input/Output On entry: the initial value of the independent variable x. Constraint: X /= XEND On exit: if G is supplied by the user, X contains the point where g(x,y) = 0.0, unless g(x,y) /= 0. 0 anywhere on the range X to XEND, in which case, X will contain XEND. If G is not supplied X contains XEND, unless an error has occured, when it contains the value of x at the error. 2: XEND -- DOUBLE PRECISION Input On entry: the final value of the independent variable. If XEND < X, integration proceeds in the negative direction. Constraint: XEND /= X. 3: M -- INTEGER Input On entry: the first dimension of the array RESULT. This will usually be equal to the number of points at which the solution is required. Constraint: M > 0. 4: N -- INTEGER Input On entry: the number of differential equations, n. Constraint: N >= 1. 5: Y(N) -- DOUBLE PRECISION array Input/Output On entry: the initial values of the solution y ,y ,...,y 1 2 n at x = X. On exit: the computed values of the solution at the final point x = X. 6: FCN -- SUBROUTINE, supplied by the user. External Procedure FCN must evaluate the functions f (i.e., the derivatives i y' ) for given values of their arguments x,y ,y ,...,y . i 1 2 n Its specification is: SUBROUTINE FCN (X, Y, F) DOUBLE PRECISION X, Y(n), F(n) where n is the actual value of N in the call of D02EJF. 1: X -- DOUBLE PRECISION Input On entry: the value of the independent variable x. 2: Y(*) -- DOUBLE PRECISION array Input On entry: the value of the variable y , for i i=1,2,...,n. 3: F(*) -- DOUBLE PRECISION array Output On exit: the value of f , for i=1,2,...,n. i FCN must be declared as EXTERNAL in the (sub)program from which D02EJF is called. Parameters denoted as Input must not be changed by this procedure. 7: PEDERV -- SUBROUTINE, supplied by the user. External Procedure PEDERV must evaluate the Jacobian of the system (that is, ddf i the partial derivatives ----) for given values of the ddy j variables x,y ,y ,...,y . 1 2 n Its specification is: SUBROUTINE PEDERV (X, Y, PW) DOUBLE PRECISION X, Y(n), PW(n,n) where n is the actual value of N in the call of D02EJF. 1: X -- DOUBLE PRECISION Input On entry: the value of the independent variable x. 2: Y(*) -- DOUBLE PRECISION array Input On entry: the value of the variable y , for i i=1,2,...,n. 3: PW(n,*) -- DOUBLE PRECISION array Output ddf i On exit: the value of ----, for i,j=1,2,...,n. ddy j If the user does not wish to supply the Jacobian, the actual argument PEDERV must be the dummy routine D02EJY . (D02EJY is included in the NAG Foundation Library and so need not be supplied by the user. The name may be implementation dependent: see the User's Note for your implementation for details). PEDERV must be declared as EXTERNAL in the (sub)program from which D02EJF is called. Parameters denoted as Input must not be changed by this procedure. 8: TOL -- DOUBLE PRECISION Input/Output On entry: TOL must be set to a positive tolerance for controlling the error in the integration. Hence TOL affects the determination of the position where g(x,y) = 0.0, if G is supplied. D02EJF has been designed so that, for most problems, a reduction in TOL leads to an approximately proportional reduction in the error in the solution. However, the actual relation between TOL and the accuracy achieved cannot be guaranteed. The user is strongly recommended to call D02EJF with more than one value for TOL and to compare the results obtained to estimate their accuracy. In the absence of any prior knowledge, the user might compare the results obtained -p -p-1 by calling D02EJF with TOL=10 and TOL=10 if p correct decimal digits are required in the solution. Constraint: TOL > 0.0. On exit: normally unchanged. However if the range X to XEND is so short that a small change in TOL is unlikely to make any change in the computed solution, then, on return, TOL has its sign changed. 9: RELABS -- CHARACTER*1 Input On entry: the type of error control. At each step in the numerical solution an estimate of the local error, EST, is made. For the current step to be accepted the following condition must be satisfied: ________________________________ / n / 1 -- 2 EST= / - > (e /((tau) *|y |+(tau) )) <=1.0 / n -- i r i a \/ i=1 where (tau) and (tau) are defined by r a RELABS (tau) (tau) r a 'M' TOL TOL 'A' 0.0 TOL 'R' TOL (epsilon) 'D' TOL (epsilon) where (epsilon) is a small machine-dependent number and e i is an estimate of the local error at y , computed i internally. If the appropriate condition is not satisfied, the step size is reduced and the solution is recomputed on the current step. If the user wishes to measure the error in the computed solution in terms of the number of correct decimal places, then RELABS should be set to 'A' on entry, whereas if the error requirement is in terms of the number of correct significant digits, then RELABS should be set to 'R'. If the user prefers a mixed error test, then RELABS should be set to 'M', otherwise if the user has no preference, RELABS should be set to the default 'D'. Note that in this case 'D' is taken to be 'R'. Constraint: RELABS = 'A', 'M', 'R' or 'D'. 10: OUTPUT -- SUBROUTINE, supplied by the user. External Procedure OUTPUT allows the user to have access to intermediate values of the computed solution at successive points specified by the user. These solution values may be returned to the user via the array RESULT if desired (this is a non-standard feature added for use with the AXIOM system). OUTPUT is initially called by D02EJF with XSOL = X (the initial value of x). The user must reset XSOL to the next point where OUTPUT is to be called, and so on at each call to OUTPUT. If, after a call to OUTPUT, the reset point XSOL is beyond XEND, D02EJF will integrate to XEND with no further calls to OUTPUT; if a call to OUTPUT is required at the point XSOL = XEND, then XSOL must be given precisely the value XEND. Its specification is: SUBROUTINE OUTPUT(XSOL,Y,COUNT,M,N,RESULT) DOUBLE PRECISION Y(N),RESULT(M,N),XSOL INTEGER M,N,COUNT 1: XSOL -- DOUBLE PRECISION Input/Output On entry: the current value of the independent variable x. On exit: the next value of x at which OUTPUT is to be called. 2: Y(N) -- DOUBLE PRECISION array Input On entry: the computed solution at the point XSOL. 3: COUNT -- INTEGER Input/Output On entry: Zero if OUTPUT has not been called before, or the previous value of COUNT. On exit: A new value of COUNT: this can be used to keep track of the number of times OUTPUT has been called. 4: M -- INTEGER Input On entry: The first dimension of RESULT. 5: N -- INTEGER Input On entry: The dimension of Y. 6: RESULT(M,N) -- DOUBLE PRECISION array Input/Output On entry: the previous contents of RESULT. On exit: RESULT may be used to return the values of the intermediate solutions to the user. OUTPUT must be declared as EXTERNAL in the (sub)program from which D02EJF is called. Parameters denoted as Input must not be changed by this procedure. 11: G -- DOUBLE PRECISION FUNCTION, supplied by the user. External Procedure G must evaluate the function g(x,y) for specified values x,y . It specifies the function g for which the first position x where g(x,y) = 0 is to be found. Its specification is: DOUBLE PRECISION FUNCTION G (X, Y) DOUBLE PRECISION X, Y(n) where n is the actual value of N in the call of D02EJF. 1: X -- DOUBLE PRECISION Input On entry: the value of the independent variable x. 2: Y(*) -- DOUBLE PRECISION array Input On entry: the value of the variable y , for i i=1,2,...,n. If the user does not require the root finding option, the actual argument G must be the dummy routine D02EJW. (D02EJW is included in the NAG Foundation Library and so need not be supplied by the user). G must be declared as EXTERNAL in the (sub)program from which D02EJF is called. Parameters denoted as Input must not be changed by this procedure. 12: W(IW) -- DOUBLE PRECISION array Workspace 13: IW -- INTEGER Input On entry: the dimension of the array W as declared in the (sub)program from which D02EJF is called. Constraint: IW>=(12+N)*N+50. 14: RESULT(M,N) -- DOUBLE PRECISION array Output On exit: the computed values of the solution at the points given by OUTPUT. 15: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. On exit: IFAIL = 0 unless the routine detects an error (see Section 6). 6. Error Indicators and Warnings Errors detected by the routine: If on entry IFAIL = 0 or -1, explanatory error messages are output on the current error message unit (as defined by X04AAF). IFAIL= 1 On entry TOL <= 0.0, or X = XEND, or N <= 0, or RELABS /= 'M', 'A', 'R', 'D'. or IW<(12+N)*N+50. IFAIL= 2 With the given value of TOL, no further progress can be made across the integration range from the current point x = X. (See Section 5 for a discussion of this error test.) The components Y(1),Y(2),...,Y(n) contain the computed values of the solution at the current point x = X. If the user has supplied G, then no point at which g(x,y) changes sign has been located up to the point x = X. IFAIL= 3 TOL is too small for D02EJF to take an initial step. X and Y (1),Y(2),...,Y(n) retain their initial values. IFAIL= 4 XSOL lies behind X in the direction of integration, after the initial call to OUTPUT, if the OUTPUT option was selected. IFAIL= 5 A value of XSOL returned by OUTPUT lies behind the last value of XSOL in the direction of integration, if the OUTPUT option was selected. IFAIL= 6 At no point in the range X to XEND did the function g(x,y) change sign, if G was supplied. It is assumed that g(x,y) = 0 has no solution. IFAIL= 7 A serious error has occurred in an internal call to C05AZF(*). Check all subroutine calls and array dimensions. Seek expert help. IFAIL= 8 A serious error has occurred in an internal call to D02XKF(*). Check all subroutine calls and array dimensions. Seek expert help. IFAIL= 9 A serious error has occurred in an internal call to D02NMF(*). Check all subroutine calls and array dimensions. Seek expert help. 7. Accuracy The accuracy of the computation of the solution vector Y may be controlled by varying the local error tolerance TOL. In general, a decrease in local error tolerance should lead to an increase in accuracy. Users are advised to choose RELABS = 'R' unless they have a good reason for a different choice. It is particularly appropriate if the solution decays. If the problem is a root-finding one, then the accuracy of the ddg ddg root determined will depend strongly on --- and ----, for ddx ddy i i=1,2,...,n. Large values for these quantities may imply large errors in the root. 8. Further Comments If more than one root is required, then to determine the second and later roots D02EJF may be called again starting a short distance past the previously determined roots. Alternatively the user may construct his own root finding code using D02QDF(*) (or the routines of the subchapter D02M-D02N), D02XKF(*) and C05AZF(*). If it is easy to code, the user should supply the routine PEDERV. However, it is important to be aware that if PEDERV is coded incorrectly, a very inefficient integration may result and possibly even a failure to complete the integration (IFAIL = 2). 9. Example We illustrate the solution of five different problems. In each case the differential system is the well-known stiff Robertson problem. 4 a'= -0.04a-10 bc 4 7 2 b'= 0.04a-10 bc -3*10 b 7 2 c'= 3*10 b with initial conditions a=1.0, b=c=0.0 at x=0.0. We solve each of the following problems with local error tolerances 1.0E-3 and 1.0 E-4. (i) To integrate to x=10.0 producing output at intervals of 2.0 until a point is encountered where a=0.9. The Jacobian is calculated numerically. (ii) As (i) but with the Jacobian calculated analytically. (iii) As (i) but with no intermediate output. (iv) As (i) but with no termination on a root-finding condition. (v) Integrating the equations as in (i) but with no intermediate output and no root-finding termination condition. The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXd02gaf}{NAG On-line Documentation: d02gaf} \beginscroll \begin{verbatim} D02GAF(3NAG) Foundation Library (12/10/92) D02GAF(3NAG) D02 -- Ordinary Differential Equations D02GAF D02GAF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose D02GAF solves the two-point boundary-value problem with assigned boundary values for a system of ordinary differential equations, using a deferred correction technique and a Newton iteration. 2. Specification SUBROUTINE D02GAF (U, V, N, A, B, TOL, FCN, MNP, X, Y, NP, 1 W, LW, IW, LIW, IFAIL) INTEGER N, MNP, NP, LW, IW(LIW), LIW, IFAIL DOUBLE PRECISION U(N,2), V(N,2), A, B, TOL, X(MNP), Y 1 (N,MNP), W(LW) EXTERNAL FCN 3. Description D02GAF solves a two-point boundary-value problem for a system of n differential equations in the interval [a,b]. The system is written in the form y' =f (x,y ,y ,...,y ) , i=1,2,...,n (1) i i 1 2 n and the derivatives are evaluated by a subroutine FCN supplied by the user. Initially, n boundary values of the variables y must i be specified (assigned), some at a and some at b. The user also supplies estimates of the remaining n boundary values and all the boundary values are used in constructing an initial approximation to the solution. This approximate solution is corrected by a finite-difference technique with deferred correction allied with a Newton iteration to solve the finite-difference equations. The technique used is described fully in Pereyra [1]. The Newton ddf i iteration requires a Jacobian matrix ---- and this is calculated ddy j by numerical differentiation using an algorithm described in Curtis et al [2]. The user supplies an absolute error tolerance and may also supply an initial mesh for the construction of the finite-difference equations (alternatively a default mesh is used). The algorithm constructs a solution on a mesh defined by adding points to the initial mesh. This solution is chosen so that the error is everywhere less than the user's tolerance and so that the error is approximately equidistributed on the final mesh. The solution is returned on this final mesh. If the solution is required at a few specific points then these should be included in the initial mesh. If on the other hand the solution is required at several specific points then the user should use the interpolation routines provided in Chapter E01 if these points do not themselves form a convenient mesh. 4. References [1] Pereyra V (1979) PASVA3: An Adaptive Finite-Difference Fortran Program for First Order Nonlinear, Ordinary Boundary Problems. Codes for Boundary Value Problems in Ordinary Differential Equations. Lecture Notes in Computer Science. (ed B Childs, M Scott, J W Daniel, E Denman and P Nelson) 76 Springer-Verlag. [2] Curtis A R, Powell M J D and Reid J K (1974) On the Estimation of Sparse Jacobian Matrices. J. Inst. Maths Applics. 13 117--119. 5. Parameters 1: U(N,2) -- DOUBLE PRECISION array Input On entry: U(i,1) must be set to the known (assigned) or estimated values of y at a and U(i,2) must be set to the i known or estimated values of y at b, for i=1,2,...,n. i 2: V(N,2) -- DOUBLE PRECISION array Input On entry: V(i,j) must be set to 0.0 if U(i,j) is a known (assigned) value and to 1.0 if U(i,j) is an estimated value, i=1,2,...,n; j=1,2. Constraint: precisely N of the V(i,j) must be set to 0.0, i.e., precisely N of the U(i,j) must be known values, and these must not be all at a or all at b. 3: N -- INTEGER Input On entry: the number of equations. Constraint: N >= 2. 4: A -- DOUBLE PRECISION Input On entry: the left-hand boundary point, a. 5: B -- DOUBLE PRECISION Input On entry: the right-hand boundary point, b. Constraint: B > A. 6: TOL -- DOUBLE PRECISION Input On entry: a positive absolute error tolerance. If a=x 0.0. 7: FCN -- SUBROUTINE, supplied by the user. External Procedure FCN must evaluate the functions f (i.e., the derivatives i y' ) at the general point x. i Its specification is: SUBROUTINE FCN (X, Y, F) DOUBLE PRECISION X, Y(n), F(n) where n is the actual value of N in the call of D02GAF. 1: X -- DOUBLE PRECISION Input On entry: the value of the argument x. 2: Y(*) -- DOUBLE PRECISION array Input On entry: the value of the argument y , for i i=1,2,...,n. 3: F(*) -- DOUBLE PRECISION array Output On exit: the values of f , for i=1,2,...,n. i FCN must be declared as EXTERNAL in the (sub)program from which D02GAF is called. Parameters denoted as Input must not be changed by this procedure. 8: MNP -- INTEGER Input On entry: the maximum permitted number of mesh points. Constraint: MNP >= 32. 9: X(MNP) -- DOUBLE PRECISION array Input/Output On entry: if NP >= 4 (see NP below), the first NP elements must define an initial mesh. Otherwise the elements of X need not be set. Constraint: A=X(1)=4 (3) On exit: X(1),X(2),...,X(NP) define the final mesh (with the returned value of NP) satisfying the relation (3). 10: Y(N,MNP) -- DOUBLE PRECISION array Output On exit: the approximate solution z (x ) satisfying (2), on j i the final mesh, that is Y(j,i)=z (x ) , i=1,2,...,NP;j=1,2,...,n, j i where NP is the number of points in the final mesh. The remaining columns of Y are not used. 11: NP -- INTEGER Input/Output On entry: determines whether a default or user-supplied mesh is used. If NP = 0, a default value of 4 for NP and a corresponding equispaced mesh X(1),X(2),...,X(NP) are used. If NP >= 4, then the user must define an initial mesh using the array X as described. Constraint: NP = 0 or 4 <= NP <= MNP. On exit: the number of points in the final (returned) mesh. 12: W(LW) -- DOUBLE PRECISION array Workspace 13: LW -- INTEGER Input On entry: the length of the array W as declared in the 2 2 calling (sub)program. Constraint: LW>=MNP*(3N +6N+2)+4N +4N 14: IW(LIW) -- INTEGER array Workspace 15: LIW -- INTEGER Input On entry: the length of the array IW as declared in the 2 calling (sub)program. Constraint: LIW>=MNP*(2N+1)+N +4N+2. 16: IFAIL -- INTEGER Input/Output For this routine, the normal use of IFAIL is extended to control the printing of error and warning messages as well as specifying hard or soft failure (see the Essential Introduction). Before entry, IFAIL must be set to a value with the decimal expansion cba, where each of the decimal digits c, b and a must have a value of 0 or 1. a=0 specifies hard failure, otherwise soft failure; b=0 suppresses error messages, otherwise error messages will be printed (see Section 6); c=0 suppresses warning messages, otherwise warning messages will be printed (see Section 6). The recommended value for inexperienced users is 110 (i.e., hard failure with all messages printed). Unless the routine detects an error (see Section 6), IFAIL contains 0 on exit. 6. Error Indicators and Warnings Errors detected by the routine: If on entry IFAIL = 0 or -1, explanatory error messages are output on the current error message unit (as defined by X04AAF). IFAIL= 1 One or more of the parameters N, TOL, NP, MNP, LW or LIW has been incorrectly set, or B <= A, or the condition (3) on X is not satisfied, or the number of known boundary values (specified by V) is not N. IFAIL= 2 The Newton iteration has failed to converge. This could be due to there being too few points in the initial mesh or to the initial approximate solution being too inaccurate. If this latter reason is suspected the user should use subroutine D02RAF instead. If the warning 'Jacobian matrix is singular' is printed this could be due to specifying zero estimated boundary values and these should be varied. This warning could also be printed in the unlikely event of the Jacobian matrix being calculated inaccurately. If the user cannot make changes to prevent the warning then subroutine D02RAF should be used. IFAIL= 3 The Newton iteration has reached round-off level. It could be, however, that the answer returned is satisfactory. This error might occur if too much accuracy is requested. IFAIL= 4 A finer mesh is required for the accuracy requested; that is MNP is not large enough. IFAIL= 5 A serious error has occurred in a call to D02GAF. Check all array subscripts and subroutine parameter lists in calls to D02GAF. Seek expert help. 7. Accuracy The solution returned by the routine will be accurate to the user's tolerance as defined by the relation (2) except in extreme circumstances. If too many points are specified in the initial mesh, the solution may be more accurate than requested and the error may not be approximately equidistributed. 8. Further Comments The time taken by the routine depends on the difficulty of the problem, the number of mesh points used (and the number of different meshes used), the number of Newton iterations and the number of deferred corrections. The user is strongly recommended to set IFAIL to obtain self- explanatory error messages, and also monitoring information about the course of the computation. The user may select the channel numbers on which this output is to appear by calls of X04AAF (for error messages) or X04ABF (for monitoring information) - see Section 9 for an example. Otherwise the default channel numbers will be used, as specified in the implementation document. A common cause of convergence problems in the Newton iteration is the user specifying too few points in the initial mesh. Although the routine adds points to the mesh to improve accuracy it is unable to do so until the solution on the initial mesh has been calculated in the Newton iteration. If the user specifies zero known and estimated boundary values, the routine constructs a zero initial approximation and in many cases the Jacobian is singular when evaluated for this approximation, leading to the breakdown of the Newton iteration. The user may be unable to provide a sufficiently good choice of initial mesh and estimated boundary values, and hence the Newton iteration may never converge. In this case the continuation facility provided in D02RAF is recommended. In the case where the user wishes to solve a sequence of similar problems, the final mesh from solving one case is strongly recommended as the initial mesh for the next. 9. Example We solve the differential equation 2 y'''=-yy''-(beta)(1-y' ) with boundary conditions y(0)=y'(0)=0, y'(10)=1 for (beta)=0.0 and (beta)=0.2 to an accuracy specified by TOL=1.0 E-3. We solve first the simpler problem with (beta)=0.0 using an equispaced mesh of 26 points and then we solve the problem with (beta)=0.2 using the final mesh from the first problem. Note the call to X04ABF prior to the call to D02GAF. The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXd02gbf}{NAG On-line Documentation: d02gbf} \beginscroll \begin{verbatim} D02GBF(3NAG) Foundation Library (12/10/92) D02GBF(3NAG) D02 -- Ordinary Differential Equations D02GBF D02GBF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose D02GBF solves a general linear two-point boundary value problem for a system of ordinary differential equations using a deferred correction technique. 2. Specification SUBROUTINE D02GBF (A, B, N, TOL, FCNF, FCNG, C, D, GAM, 1 MNP, X, Y, NP, W, LW, IW, LIW, IFAIL) INTEGER N, MNP, NP, LW, IW(LIW), LIW, IFAIL DOUBLE PRECISION A, B, TOL, C(N,N), D(N,N), GAM(N), X(MNP), 1 Y(N,MNP), W(LW) EXTERNAL FCNF, FCNG 3. Description D02GBF solves the linear two-point boundary value problem for a system of n ordinary differential equations in the interval [a,b]. The system is written in the form y'=F(x)y+g(x) (1) and the boundary conditions are written in the form Cy(a)+Dy(b)=(gamma) (2) Here F(x), C and D are n by n matrices, and g(x) and (gamma) are n-component vectors. The approximate solution to (1) and (2) is found using a finite-difference method with deferred correction. The algorithm is a specialisation of that used in subroutine D02RAF which solves a nonlinear version of (1) and (2). The nonlinear version of the algorithm is described fully in Pereyra [1]. The user supplies an absolute error tolerance and may also supply an initial mesh for the construction of the finite-difference equations (alternatively a default mesh is used). The algorithm constructs a solution on a mesh defined by adding points to the initial mesh. This solution is chosen so that the error is everywhere less than the user's tolerance and so that the error is approximately equidistributed on the final mesh. The solution is returned on this final mesh. If the solution is required at a few specific points then these should be included in the initial mesh. If, on the other hand, the solution is required at several specific points, then the user should use the interpolation routines provided in Chapter E01 if these points do not themselves form a convenient mesh. 4. References [1] Pereyra V (1979) PASVA3: An Adaptive Finite-Difference Fortran Program for First Order Nonlinear, Ordinary Boundary Problems. Codes for Boundary Value Problems in Ordinary Differential Equations. Lecture Notes in Computer Science. (ed B Childs, M Scott, J W Daniel, E Denman and P Nelson) 76 Springer-Verlag. 5. Parameters 1: A -- DOUBLE PRECISION Input On entry: the left-hand boundary point, a. 2: B -- DOUBLE PRECISION Input On entry: the right-hand boundary point, b. Constraint: B > A. 3: N -- INTEGER Input On entry: the number of equations; that is n is the order of system (1). Constraint: N >= 2. 4: TOL -- DOUBLE PRECISION Input On entry: a positive absolute error tolerance. If a=x 0.0. 5: FCNF -- SUBROUTINE, supplied by the user. External Procedure FCNF must evaluate the matrix F(x) in (1) at a general point x. Its specification is: SUBROUTINE FCN (X, F) DOUBLE PRECISION X, F(n,n) where n is the actual value of N in the call of D02GBF. 1: X -- DOUBLE PRECISION Input On entry: the value of the independent variable x. 2: F(n,n) -- DOUBLE PRECISION array Output On exit: the (i,j)th element of the matrix F(x), for i,j=1,2,...,n. (See Section 9 for an example.) FCN must be declared as EXTERNAL in the (sub)program from which D02GBF is called. Parameters denoted as Input must not be changed by this procedure. 6: FCNG -- SUBROUTINE, supplied by the user. External Procedure FCNG must evaluate the vector g(x) in (1) at a general point x. Its specification is: SUBROUTINE FCNG (X, G) DOUBLE PRECISION X, G(n) where n is the actual value of N in the call of D02GBF. 1: X -- DOUBLE PRECISION Input On entry: the value of the independent variable x. 2: G(*) -- DOUBLE PRECISION array Output On exit: the ith element of the vector g(x), for i=1,2,...,n. (See Section Section 9 for an example.) FCNG must be declared as EXTERNAL in the (sub)program from which D02GBF is called. Parameters denoted as Input must not be changed by this procedure. 7: C(N,N) -- DOUBLE PRECISION array Input/Output 8: D(N,N) -- DOUBLE PRECISION array Input/Output 9: GAM(N) -- DOUBLE PRECISION array Input/Output On entry: the arrays C and D must be set to the matrices C and D in (2). GAM must be set to the vector (gamma) in (2). On exit: the rows of C and D and the components of GAM are re-ordered so that the boundary conditions are in the order: (i) conditions on y(a) only; (ii) condition involving y(a) and y(b); and (iii) conditions on y(b) only. The routine will be slightly more efficient if the arrays C, D and GAM are ordered in this way before entry, and in this event they will be unchanged on exit. Note that the problems (1) and (2) must be of boundary value type, that is neither C nor D may be identically zero. Note also that the rank of the matrix [C,D] must be n for the problem to be properly posed. Any violation of these conditions will lead to an error exit. 10: MNP -- INTEGER Input On entry: the maximum permitted number of mesh points. Constraint: MNP >= 32. 11: X(MNP) -- DOUBLE PRECISION array Input/Output On entry: if NP >= 4 (see NP below), the first NP elements must define an initial mesh. Otherwise the elements of x need not be set. Constraint: A=X(1)=4. (4) On exit: X(1),X(2),...,X(NP) define the final mesh (with the returned value of NP) satisfying the relation (4). 12: Y(N,MNP) -- DOUBLE PRECISION array Output On exit: the approximate solution z(x) satisfying (3), on the final mesh, that is Y(j,i)=z (x ) , i=1,2,...,NP;j=1,2,...,n j i where NP is the number of points in the final mesh. The remaining columns of Y are not used. 13: NP -- INTEGER Input/Output On entry: determines whether a default mesh or user-supplied mesh is used. If NP = 0, a default value of 4 for NP and a corresponding equispaced mesh X(1),X(2),...,X(NP) are used. If NP >= 4, then the user must define an initial mesh X as in (4) above. On exit: the number of points in the final (returned) mesh. 14: W(LW) -- DOUBLE PRECISION array Workspace 15: LW -- INTEGER Input On entry: the length of the array W, Constraint: 2 2 LW>=MNP*(3N +5N+2)+3N +5N. 16: IW(LIW) -- INTEGER array Workspace 17: LIW -- INTEGER Input On entry: the length of the array IW. Constraint: LIW>=MNP*(2N+1)+N. 18: IFAIL -- INTEGER Input/Output For this routine, the normal use of IFAIL is extended to control the printing of error and warning messages as well as specifying hard or soft failure (see the Essential Introduction). Before entry, IFAIL must be set to a value with the decimal expansion cba, where each of the decimal digits c, b and a must have a value of 0 or 1. a=0 specifies hard failure, otherwise soft failure; b=0 suppresses error messages, otherwise error messages will be printed (see Section 6); c=0 suppresses warning messages, otherwise warning messages will be printed (see Section 6). The recommended value for inexperienced users is 110 (i.e., hard failure with all messages printed). Unless the routine detects an error (see Section 6), IFAIL contains 0 on exit. 6. Error Indicators and Warnings Errors detected by the routine: For each error, an explanatory error message is output on the current error message unit (as defined by X04AAF), unless suppressed by the value of IFAIL on entry. IFAIL= 1 One or more of the parameters N, TOL, NP, MNP, LW or LIW is incorrectly set, B <= A or the condition (4) on X is not satisfied. IFAIL= 2 There are three possible reasons for this error exit to be taken: (i) one of the matrices C or D is identically zero (that is the problem is of initial value and not boundary value type). In this case, IW(1) = 0 on exit; (ii) a row of C and the corresponding row of D are identically zero (that is the boundary conditions are rank deficient). In this case, on exit IW(1) contains the index of the first such row encountered; and (iii more than n of the columns of the n by 2n matrix [C,D ) ] are identically zero (that is the boundary conditions are rank deficient). In this case, on exit IW(1) contains minus the number of non-identically zero columns. IFAIL= 3 The routine has failed to find a solution to the specified accuracy. There are a variety of possible reasons including: (i) the boundary conditions are rank deficient, which may be indicated by the message that the Jacobian is singular. However this is an unlikely explanation for the error exit as all rank deficient boundary conditions should lead instead to error exits with either IFAIL = 2 or IFAIL = 5; see also (iv) below; (ii) not enough mesh points are permitted in order to attain the required accuracy. This is indicated by NP = MNP on return from a call to D02GBF. This difficulty may be aggravated by a poor initial choice of mesh points; (iii) the accuracy requested cannot be attained on the computer being used; and (iv) an unlikely combination of values of F(x) has led to a singular Jacobian. The error should not persist if more mesh points are allowed. IFAIL= 4 A serious error has occurred in a call to D02GBF. Check all array subscripts and subroutine parameter lists in calls to D02GBF. Seek expert help. IFAIL= 5 There are two possible reasons for this error exit which occurs when checking the rank of the boundary conditions by reduction to a row echelon form: (i) at least one row of the n by 2n matrix [C,D] is a linear combination of the other rows and hence the boundary conditions are rank deficient. The index of the first such row encountered is given by IW(1) on exit; and (ii) as (i) but the rank deficiency implied by this error exit has only been determined up to a numerical tolerance. Minus the index of the first such row encountered is given by IW(1) on exit. In case (ii) above there is some doubt as to the rank deficiency of the boundary conditions. However even if the boundary conditions are not rank deficient they are not posed in a suitable form for use with this routine. For example, if ((gamma) ) (1 0 ) (1 0) ( 1) C=(1 (epsilon)) , D=(1 0) , (gamma)=((gamma) ) ( 2) and (epsilon) is small enough, this error exit is likely to be taken. A better form for the boundary conditions in this case would be ((gamma) ) ( 1 ) (1 0) (1 0) ( -1 ) C=(0 1) , D=(0 0) , (gamma)=((epsilon) ((gamma) -(gamma) )) ( 2 1 ) 7. Accuracy The solution returned by the routine will be accurate to the user's tolerance as defined by the relation (3) except in extreme circumstances. If too many points are specified in the initial mesh, the solution may be more accurate than requested and the error may not be approximately equidistributed. 8. Further Comments The time taken by the routine depends on the difficulty of the problem, the number of mesh points (and meshes) used and the number of deferred corrections. The user is strongly recommended to set IFAIL to obtain self- explanatory error messages, and also monitoring information about the course of the computation. The user may select the channel numbers on which this output is to appear by calls of X04AAF (for error messages) or X04ABF (for monitoring information) - see Section 9 for an example. Otherwise the default channel numbers will be used, as specified in the implementation document. In the case where the user wishes to solve a sequence of similar problems, the use of the final mesh from one case is strongly recommended as the initial mesh for the next. 9. Example We solve the problem (written as a first order system) (epsilon)y''+y'=0 with boundary conditions y(0)=0, y(1)=1 -1 -2 for the cases (epsilon)=10 and (epsilon)=10 using the default initial mesh in the first case, and the final mesh of the first case as initial mesh for the second (more difficult) case. We give the solution and the error at each mesh point to illustrate the accuracy of the method given the accuracy request TOL=1.0E-3. Note the call to X04ABF prior to the call to D02GBF. The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXd02kef}{NAG On-line Documentation: d02kef} \beginscroll \begin{verbatim} D02KEF(3NAG) D02KEF D02KEF(3NAG) D02 -- Ordinary Differential Equations D02KEF D02KEF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose D02KEF finds a specified eigenvalue of a regular singular second- order Sturm-Liouville system on a finite or infinite range, using a Pruefer transformation and a shooting method. It also reports values of the eigenfunction and its derivatives. Provision is made for discontinuities in the coefficient functions or their derivatives. 2. Specification SUBROUTINE D02KEF (XPOINT, M, MATCH, COEFFN, BDYVAL, K, 1 TOL, ELAM, DELAM, HMAX, MAXIT, MAXFUN, 2 MONIT, REPORT, IFAIL) INTEGER M, MATCH, K, MAXIT, MAXFUN, IFAIL DOUBLE PRECISION XPOINT(M), TOL, ELAM, DELAM, HMAX(2,M) EXTERNAL COEFFN, BDYVAL, MONIT, REPORT 3. Description D02KEF has essentially the same purpose as D02KDF(*) with minor modifications to enable values of the eigenfunction to be obtained after convergence to the eigenvalue has been achieved. ~~~~~~~~ It first finds a specified eigenvalue (lambda) of a Sturm- Liouville system defined by a self-adjoint differential equation of the second-order (p(x)y')'+q(x;(lambda))y=0, a= |f(x)g(x)w(x)dx, / a with respect to which the eigenfunctions are mutually orthogonal. This normalisation of y is only approximate, but experience shows that S generally differs from unity by only one or two per cent. During this final integration the REPORT routine supplied by the user is called at each integration mesh point x. Sufficient information is returned to permit the user to compute y(x) and y'(x) for printing or plotting. For reasons described in Section 8.2, D02KEF passes across to REPORT, not y and y', but the Pru efer variables (beta), (phi) and (rho) on which the numerical method is based. Their relationship to y and y' is given by the equations ______ ( (rho)) ( (phi)) p(x)y'=\/(beta)exp( -----)cos( -----); ( 2 ) ( 2 ) 1 ( (rho)) ( (phi)) y= --------exp( -----)sin( -----). ______ ( 2 ) ( 2 ) \/(beta) For the theoretical basis of the numerical method to be valid, the following conditions should hold on the coefficient functions: (a) p(x) must be non-zero and of one sign throughout the interval (a,b); and, ddq (b) ---------- must be of one sign throughout (a,b) for all dd(lambda) relevant values of (lambda), and must not be identically zero as x varies, for any (lambda). Points of discontinuity in the functions p and q or their derivatives are allowed, and should be included as 'break-points' in the array XPOINT. A good account of the theory of Sturm-Liouville systems, with some description of Pruefer transformations, is given in Birkhoff and Rota [4], Chapter X. An introduction for the user of Pruefer transformations for the numerical solution of eigenvalue problems arising from physics and chemistry is Bailey [2]. The scaled Pruefer method is fairly recent, and is described in a short note by Pryce [6] and in some detail in the technical report [5]. 4. References [1] Abramowitz M and Stegun I A (1968) Handbook of Mathematical Functions. Dover Publications. [2] Bailey P B (1966) Sturm-Liouville Eigenvalues via a Phase Function. SIAM J. Appl. Math . 14 242--249. [3] Banks D O and Kurowski I (1968) Computation of Eigenvalues of Singular Sturm-Liouville Systems. Math. Computing. 22 304--310. [4] Birkhoff G and Rota G C (1962) Ordinary Differential Equations. Ginn & Co., Boston and New York. [5] Pryce J D (1981) Two codes for Sturm-Liouville problems. Technical Report CS-81-01. Dept of Computer Science, Bristol University . [6] Pryce J D and Hargrave B A (1977) The Scale Pruefer Method for one-parameter and multi-parameter eigenvalue problems in ODEs. Inst. Math. Appl., Numerical Analysis Newsletter. 1(3) 5. Parameters 1: XPOINT(M) -- DOUBLE PRECISION array Input On entry: the points where the boundary conditions computed by BDYVAL are to be imposed, and also any break-points, i.e., XPOINT(1) to XPOINT(m) must contain values x ,...,x 1 m such that x <=x x . Similarly x =x (=b) 2 1 m-1 m if the right-hand end-point is regular, and x = 4. 3: MATCH -- INTEGER Input/Output On entry: MATCH must be set to the index of the 'break-point ' to be used as the matching point (see Section 8.3). If MATCH is set to a value outside the range [2,m-1] then a default value is taken, corresponding to the break-point nearest the centre of the interval [XPOINT(2),XPOINT(m-1)]. On exit: the index of the break-point actually used as the matching point. 4: COEFFN -- SUBROUTINE, supplied by the user. External Procedure COEFFN must compute the values of the coefficient functions p(x) and q(x;(lambda)) for given values of x and (lambda). Section 3 states conditions which p and q must satisfy. Its specification is: SUBROUTINE COEFFN (P, Q, DQDL, X, ELAM, JINT) DOUBLE PRECISION P, Q, DQDL, X, ELAM INTEGER JINT 1: P -- DOUBLE PRECISION Output On exit: the value of p(x) for the current value of x. 2: Q -- DOUBLE PRECISION Output On exit: the value of q(x;(lambda)) for the current value of x and the current trial value of (lambda). 3: DQDL -- DOUBLE PRECISION Output ddq On exit: the value of ----------(x;(lambda)) for the dd(lambda) current value of x and the current trial value of (lambda). However DQDL is only used in error estimation and an approximation (say to within 20%) will suffice. 4: X -- DOUBLE PRECISION Input On entry: the current value of x. 5: ELAM -- DOUBLE PRECISION Input On entry: the current trial value of the eigenvalue parameter (lambda). 6: JINT -- INTEGER Input On entry: the index j of the sub-interval i (see j specification of XPOINT) in which x lies. See Section 8.4 and Section 9 for examples. COEFFN must be declared as EXTERNAL in the (sub)program from which D02KEF is called. Parameters denoted as Input must not be changed by this procedure. 5: BDYVAL -- SUBROUTINE, supplied by the user. External Procedure BDYVAL must define the boundary conditions. For each end- point, BDYVAL must return (in YL or YR) values of y(x) and p(x)y'(x) which are consistent with the boundary conditions at the end-points; only the ratio of the values matters. Here x is a given point (XL or XR) equal to, or close to, the end-point. For a regular end-point (a, say), x=a; and a boundary condition of the form c y(a)+c y'(a)=0 1 2 can be handled by returning constant values in YL, e.g. YL(1)=c and YL(2)=-c p(a). 2 1 For a singular end-point however, YL(1) and YL(2) will in general be functions of XL and ELAM, and YR(1) and YR(2) functions of XR and ELAM, usually derived analytically from a power-series or asymptotic expansion. Examples are given in Section 8.5 Section 9. Its specification is: SUBROUTINE BDYVAL (XL, XR, ELAM, YL, YR)) DOUBLE PRECISION XL, XR, ELAM, YL(3), YR(3) 1: XL -- DOUBLE PRECISION Input On entry: if a is a regular end-point of the system (so that a=x =x ), then XL contains a. If a is a singular 1 2 point (so that a<=x = 0. 7: TOL -- DOUBLE PRECISION Input On entry: the tolerance parameter which determines the accuracy of the computed eigenvalue. The error estimate held in DELAM on exit satisfies the mixed absolute/relative error test DELAM<=TOL*max(1.0,|ELAM|) (*) where ELAM is the final estimate of the eigenvalue. DELAM is usually somewhat smaller than the right-hand side of (*) but not several orders of magnitude smaller. Constraint: TOL > 0.0. 8: ELAM -- DOUBLE PRECISION Input/Output ~~~~~~~~ On entry: an initial estimate of the eigenvalue (lambda). On exit: the final computed estimate, whether or not an error occurred. 9: DELAM -- DOUBLE PRECISION Input/Output On entry: an indication of the scale of the problem in the (lambda)-direction. DELAM holds the initial 'search step' (positive or negative). Its value is not critical but the first two trial evaluations are made at ELAM and ELAM + DELAM, so the routine will work most efficiently if the eigenvalue lies between these values. A reasonable choice (if a closer bound is not known) is half the distance between adjacent eigenvalues in the neighbourhood of the one sought. In practice, there will often be a problem, similar to the one in hand but with known eigenvalues, which will help one to choose initial values for ELAM and DELAM. If DELAM = 0.0 on entry, it is given the default value of 0.25*max(1.0,|ELAM|). On exit: with IFAIL = 0, DELAM holds an estimate of the absolute error in the computed ~~~~~~~~ eigenvalue, that is |(lambda)-ELAM|~=DELAM. (In Section 8.2 we discuss the assumptions under which this is true.) The true error is rarely more than twice, or less than a tenth, of the estimated error. With IFAIL /= 0, DELAM may hold an estimate of the error, or its initial value, depending on the value of IFAIL. See Section 6 for further details. 10: HMAX(2,M) -- DOUBLE PRECISION array Input/Output On entry: HMAX(1,j) a maximum step size to be used by the differential equation code in the jth sub-interval i (as j described in the specification of parameter XPOINT), for j=1,2,...,m-3. If it is zero the routine generates a maximum step size internally. It is recommended that HMAX(1,j) be set to zero unless the coefficient functions p and q have features (such as a narrow peak) within the jth sub-interval that could be ' missed' if a long step were taken. In such a case HMAX(1,j) should be set to about half the distance over which the feature should be observed. Too small a value will increase the computing time for the routine. See Section 8 for further suggestions. The rest of the array is used as workspace. On exit: HMAX( 1,m-1) and HMAX(1,m) contain the sensitivity coefficients (sigma) ,(sigma) , described in Section 8.6. Other entries l r contain diagnostic output in case of an error (see Section 6 ). 11: MAXIT -- INTEGER Input/Output On entry: a bound on n , the number of root-finding r iterations allowed, that is the number of trial values of (lambda) that are used. If MAXIT <= 0, no such bound is assumed. (See also under MAXFUN.) Suggested value: MAXIT = 0. On exit: MAXIT will have been decreased by the number of iterations actually performed, whether or not it was positive on entry. 12: MAXFUN -- INTEGER Input On entry: a bound on n , the number of calls to COEFFN made f in any one root-finding iteration. If MAXFUN <= 0, no such bound is assumed. Suggested value: MAXFUN = 0. MAXFUN and MAXIT may be used to limit the computational cost of a call to D02KEF, which is roughly proportional to n *n . r f 13: MONIT -- SUBROUTINE, supplied by the user. External Procedure MONIT is called by D02KEF at the end of each root-finding iteration and allows the user to monitor the course of the computation by printing out the parameters (see Section 8 for an example). If no monitoring is required, the dummy subroutine D02KAY may be used. (D02KAY is included in the NAG Foundation Library). Its specification is: SUBROUTINE MONIT (MAXIT, IFLAG, ELAM, FINFO) INTEGER MAXIT, IFLAG DOUBLE PRECISION ELAM, FINFO(15) 1: MAXIT -- INTEGER Input On entry: the current value of the parameter MAXIT of D02KEF; this is decreased by one at each iteration. 2: IFLAG -- INTEGER Input On entry: IFLAG describes what phase the computation is in, as follows: IFLAG < 0 an error occurred in the computation of the ' miss-distance' at this iteration; an error exit from D02KEF with IFAIL =-IFLAG will follow. IFLAG = 1 the routine is trying to bracket the eigenvalue ~~~~~~~~ (lambda). IFLAG = 2 the routine is converging to the eigenvalue ~~~~~~~~ (lambda) (having already bracketed it). 3: ELAM -- DOUBLE PRECISION Input On entry: the current trial value of (lambda). 4: FINFO(15) -- DOUBLE PRECISION array Input On entry: information about the behaviour of the shooting method, and diagnostic information in the case of errors. It should not normally be printed in full if no error has occurred (that is, if IFLAG > 0), though the first few components may be of interest to the user. In case of an error (IFLAG < 0) all the components of FINFO should be printed. The contents of FINFO are as follows: FINFO(1): the current value of the 'miss-distance' or ' residual' function f((lambda)) on which the shooting method is based. FINFO(1) is set to zero if IFLAG < 0. FINFO(2): an estimate of the quantity dd(lambda) defined as follows. Consider the perturbation in the miss-distance f((lambda)) that would result if the local error, in the solution of the differential equation, were always positive and equal to its maximum permitted value. Then dd(lambda) is the perturbation in (lambda) that would have the same effect on f((lambda)) . Thus, at the zero of f((lambda)),|dd(lambda)| is an approximate bound on the perturbation of the zero (that is the eigenvalue) caused by errors in numerical solution. If dd(lambda) is very large then it is possible that there has been a programming error in COEFFN such that q is independent of (lambda). If this is the case, an error exit with IFAIL = 5 should follow. FINFO(2) is set to zero if IFLAG < 0. FINFO(3): the number of internal iterations, using the same value of (lambda) and tighter accuracy tolerances, needed to bring the accuracy (that is the value of dd(lambda)) to an acceptable value. Its value should normally be 1.0, and should almost never exceed 2.0. FINFO(4): the number of calls to COEFFN at this iteration. FINFO(5): the number of successful steps taken by the internal differential equation solver at this iteration. A step is successful if it is used to advance the integration (cf. COUT(8) in specification of D02PAF(*)). FINFO(6): the number of unsuccessful steps used by the internal integrator at this iteration (cf. COUT(9) in specification of D02PAF(*)). FINFO(7): the number of successful steps at the maximum step size taken by the internal integrator at this iteration (cf. COUT(3) in specification of D02PAF(*)). FINFO(8): is not used. FINFO(9) to FINFO(15): set to zero, unless IFLAG < 0 in which case they hold the following values describing the point of failure: FINFO(9): contains the index of the sub-interval where failure occurred, in the range 1 to m-3. In case of an error in BDYVAL, it is set to 0 or m-2 depending on whether the left or right boundary condition caused the error. FINFO(10): the value f the independent variable x, the point at which error occurred. In case of an error in BDYVAL, it is set to the value of XL or XR as appropriate (see the specification of BDYVAL). FINFO(11), FINFO(12), FINFO(13): the current values of the Pruefer dependent variables (beta), (phi) and (rho) respectively. These are set to zero in case of an error in BDYVAL. FINFO(14): the local-error tolerance being used by the internal integrator at the point of failure. This is set to zero in the case of an error in BDYVAL. FINFO(15): the last integration mesh point. This is set to zero in the case of an error in BDYVAL. MONIT must be declared as EXTERNAL in the (sub)program from which D02KEF is called. Parameters denoted as Input must not be changed by this procedure. 14: REPORT -- SUBROUTINE, supplied by the user. External Procedure This routine provides the means by which the user may compute the eigenfunction y(x) and its derivative at each integration mesh point x. (See Section 8 for an example). Its specification is: SUBROUTINE REPORT (X, V, JINT) INTEGER JINT DOUBLE PRECISION X, V(3) 1: X -- DOUBLE PRECISION Input On entry: the current value of the independent variable x. See Section 8.3 for the order in which values of x are supplied. 2: V(3) -- DOUBLE PRECISION array Input On entry: V(1), V(2), V(3) hold the current values of the Pruefer variables (beta), (phi), (rho) respectively. 3: JINT -- INTEGER Input On entry: JINT indicates the sub-interval between break-points in which X lies exactly as for the routine COEFFN, except that at the extreme left end-point (when x = XPOINT(2)) JINT is set to 0 and at the extreme right end-point (when x=x =XPOINT(m-1)) JINT is set to r m-2. REPORT must be declared as EXTERNAL in the (sub)program from which D02KEF is called. Parameters denoted as Input must not be changed by this procedure. 15: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. On exit: IFAIL = 0 unless the routine detects an error (see Section 6). 6. Error Indicators and Warnings Errors detected by the routine: IFAIL= 1 A parameter error. All parameters (except IFAIL) are left unchanged. The reason for the error is shown by the value of HMAX(2,1) as follows: HMAX(2,1) = 1: M < 4; HMAX(2,1) = 2: K < 0; HMAX(2,1) = 3: TOL <= 0.0; HMAX(2,1) = 4: XPOINT(1) to XPOINT(m) are not in ascending order. HMAX(2,2) gives the position i in XPOINT where this was detected. IFAIL= 2 At some call to BDYVAL, invalid values were returned, that is, either YL(1) = YL(2) = 0.0, or YR(1) = YR(2) = 0.0 (a programming error in BDYVAL). See the last call of MONIT for details. This error exit will also occur if p(x) is zero at the point where the boundary condition is imposed. Probably BDYVAL was called with XL equal to a singular end-point a or with XR equal to a singular end-point b. IFAIL= 3 At some point between XL and XR the value of p(x) computed by COEFFN became zero or changed sign. See the last call of MONIT for details. IFAIL= 4 MAXIT > 0 on entry, and after MAXIT iterations the eigenvalue had not been found to the required accuracy. IFAIL= 5 The 'bracketing' phase (with parameter IFLAG of MONIT equal to 1) failed to bracket the eigenvalue within ten iterations. This is caused by an error in formulating the problem (for example, q is independent of (lambda)), or by very poor initial estimates of ELAM, DELAM. On exit ELAM and ELAM + DELAM give the end-points of the interval within which no eigenvalue was located by the routine. IFAIL= 6 MAXFUN > 0 on entry, and the last iteration was terminated because more than MAXFUN calls to COEFFN were used. See the last call of MONIT for details. IFAIL= 7 To obtain the desired accuracy the local error tolerance was set so small at the start of some sub-interval that the differential equation solver could not choose an initial step size large enough to make significant progress. See the last call of MONIT for diagnostics. IFAIL= 8 At some point inside a sub-interval the step size in the differenital equation solver was reduced to a value too small to make significant progress (for the same reasons as with IFAIL = 7). This could be due to pathological behaviour of p(x) and q(x;(lambda)) or to an unreasonable accuracy requirement or to the current value of (lambda) making the equations 'stiff'. See the last call of MONIT for details. IFAIL= 9 TOL is too small for the problem being solved and the machine precision is being used. The final value of ELAM should be a very good approximation to the eigenvalue. IFAIL= 10 C05AZF(*), called by D02KEF, has terminated with the error exit corresponding to a pole of the residual function f((lambda)). This error exit should not occur, but if it does, try solving the problem again with a smaller value for TOL. IFAIL= 11 A serious error has occurred in an internal call to D02KDY. Check all subroutine calls and array dimensions. Seek expert help. IFAIL= 12 A serious error has occurred in an internal call to C05AZF(*). Check all subroutine calls and array dimensions. Seek expert help. HMAX(2,1) holds the failure exit number from the routine where the failure occurred. In the case of a failure in C05AZF(*), HMAX(2,2) holds the value of parameter IND of C05AZF(*). Note: error exits with IFAIL = 2, 3, 6, 7, 8, 11 are caused by being unable to set up or solve the differential equation at some iteration, and will be immediately preceded by a call of MONIT giving diagnostic information. For other errors, diagnostic information is contained in HMAX(2,j), for j=1,2,...,m, where appropriate. 7. Accuracy See the discussion in Section 8.2. 8. Further Comments 8.1. Timing The time taken by the routine depends on the complexity of the coefficient functions, whether they or their derivatives are rapidly changing, the tolerance demanded, and how many iterations are needed to obtain convergence. The amount of work per iteration is roughly doubled when TOL is divided by 16. To make the most economical use of the routine, one should try to obtain good initial values for ELAM and DELAM, and, where appropriate, good asymptotic formulae. The boundary matching points should not be set unnecessarily close to singular points. The extra time needed to compute the eigenfunction is principally the cost of one additional integration once the eigenvalue has been found. 8.2. General Description of the Algorithm A shooting method, for differential equation problems containing unknown parameters, relies on the construction of a 'miss- distance function', which for given trial values of the parameters measures how far the conditions of the problem are from being met. The problem is then reduced to one of finding the values of the parameters for which the miss-distance function is zero, that is to a root-finding process. Shooting methods differ mainly in how the miss-distance is defined. This routine defines a miss-distance f((lambda)) based on the rotation around the origin of the point P(x)=(p(x)y'(x),y(x)) in the Phase Plane as the solution proceeds from a to b. The boundary-conditions define the ray (i.e., two-sided line through the origin) on which p(x) should start, and the ray on which it should finish. The eigenvalue index k defines the total number of half-turns it should make. Numerical solution is actually done by matching point x=c. Then f((lambda)) is taken as the angle between the rays to the two resulting points P (c) and P (c). A a b relative scaling of the py' and y axes, based on the behaviour of the coefficient functions p and q, is used to improve the numerical behaviour. Please see figure in printed Reference Manual The resulting function f((lambda)) is monotonic over - ddq infty<(lambda)0 and decreasing dd(lambda) ddq if ----------<0, with a unique zero at the desired eigenvalue dd(lambda) ~~~~~~~~ (lambda). The routine measures f((lambda)) in units of a half- turn. This means that as (lambda) increases, f((lambda)) varies by about 1 as each eigenvalue is passed. (This feature implies that the values of f((lambda)) at successive iterations - especially in the early stages of the iterative process - can be used with suitable extrapolation or interpolation to help the choice of initial estimates for eigenvalues near to the one currently being found.) The routine actually computes a value for f((lambda)) with errors, arising from the local errors of the differential equation code and from the asymptotic formulae provided by the user if singular points are involved. However, the error estimate output in DELAM is usually fairly realistic, in that the actual ~~~~~~~~ error |(lambda)-ELAM| is within an order of magnitude of DELAM. We pass the values of (beta), (phi), (rho) across through REPORT rather than converting them to values of y, y' inside D02KEF, for the following reasons. First, there may be cases where auxiliary quantities can be more accurately computed from the Pruefer variables than from y and y'. Second, in singular problems on an infinite interval y and y' may underflow towards the end of the range, whereas the Pruefer variables remain well-behaved. Third, with high-order eigenvalues (and therefore highly oscillatory eigenfunctions) the eigenfunction may have a complete oscillation (or more than one oscillation) between two mesh points, so that values of y and y' at mesh points give a very poor representation of the curve. The probable behaviour of the Pruefer variables in this case is that (beta) and (rho) vary slowly whilst (phi) increases quickly: for all three Pruefer variables linear interpolation between the values at adjacent mesh points is probably sufficiently accurate to yield acceptable intermediate values of (beta), (phi), (rho) (and hence of y,y') for graphical purposes. Similar considerations apply to the exponentially decaying 'tails Here (phi) has approximately constant value whilst (rho) increases rapidly in the direction of integration, though the step length is generally fairly small over such a range. If the solution is output through REPORT at x-values which are too widely spaced, the step length can be controlled by choosing HMAX suitably, or, preferably, by reducing TOL. Both these choices will lead to more accurate eigenvalues and eigenfunctions but at some computational cost. 8.3. The Position of the Shooting Matching Point c This point is always one of the values x in array XPOINT. It may i be specified using the parameter MATCH. The default value is chosen to be the value of that x ,2<=i<=m-1, that lies closest to i the middle of the interval [x ,x ]. If there is a tie, the 2 m-1 rightmost candidate is chosen. In particular if there are no break-points then c=x (=x ) - that is the shooting is from left m-1 3 to right in this case. A break-point may be inserted purely to move c to an interior point of the interval, even though the form of the equations does not require it. This often speeds up convergence especially with singular problems. Note that the shooting method used by the code integrates first from the left-hand end x , then from the right-hand end x , to l r meet at the matching point c in the middle. This will of course be reflected in printed or graphical output. The diagram shows a possible sequence of nine mesh points (tau) through (tau) in 1 9 the order in which they appear, assuming there are just two sub- intervals (so m=5). Figure 1 Please see figure in printed Reference Manual Since the shooting method usually fails to match up the two 'legs p(x)y' or both, at the matching point c. The code in fact 'shares large jump does not imply an inaccurate eigenvalue, but implies either (a) a badly chosen matching point: if q(x;(lambda)) has a ' humped' shape, c should be chosen near the maximum value of q, especially if q is negative at the ends of the interval. (b) An inherently ill-conditioned problem, typically one where another eigenvalue is pathologically close to the one being sought. In this case it is extremely difficult to obtain an accurate eigenfunction. In Section 9 below, we find the 11th eigenvalue and corresponding eigenfunction of the equation 2 y''+((lambda)-x-2/x )y=0 on 04), over some interval 'approximating to (-infty,infty)', say [-20, 20]. Here we need break-points at +- 4, forming three sub- intervals i =[-20,-4], i =[-4,4],i =[4,20]. The code could be 1 2 3 SUBROUTINE COEFFN(P,Q,DQDL,X,ELAM,JINT) IF (JINT.EQ.2) THEN Q = ELAM + X*X - 10.0E0 ELSE Q = ELAM + 6.0E0/ABS(X) ENDIF P = 1.0E0 DQDL = 1.0E0 RETURN END The array XPOINT would contain the values x , -20.0, -4.0, +4.0, 1 +20.0, x and m would be 6. The choice of appropriate values for 6 x and x depends on the form of the asymptotic formula computed 1 6 by BDYVAL and the technique is discussed in the next subsection. Example C 2 -100x y''+(lambda)(1-2e )y=0, over -10<=x<=10 Here q(x;(lambda)) is nearly constant over the range except for a sharp inverted spike over approximately -0.1<=x<=0.1. There is a danger that the routine will build up to a large step size and ' step over' the spike without noticing it. By using break-points - say at +- 0.5 - one can restrict the step size near the spike without impairing the efficiency elsewhere. The code for COEFFN could be SUBROUTINE COEFFN(P,Q,DQDL,X,ELAM,JINT) P = 1.0E0 DQDL = 1.0E0 - 2.0E0*EXP(-100.0E0*X*X) Q = ELAM*DQDL RETURN END XPOINT might contain -0.0, -10.0, -0.5, 0.5, 10.0, 10.0 (assuming +- 10 are regular points) and m would be 6. HMAX(1,j), j=1,2,3 might contain 0.0, 0.1 and 0.0. 8.5. Examples of Boundary Conditions at Singular Points Quoting from Bailey [2] page 243: 'Usually... the differential equation has two essentially different types of solution near a singular point, and the boundary condition there merely serves to distinguish one kind from the other. This is the case in all the standard examples of mathematical physics.' In most cases the behaviour of the ratio p(x)y'/y near the point is quite different for the two types of solution. Essentially what the user provides through his BDYVAL routine is an approximation to this ratio, valid as x tends to the singular point (SP). The user must decide (a) how accurate to make this approximation or asymptotic formula, for example how many terms of a series to use, and (b) where to place the boundary matching point (BMP) at which the numerical solution of the differential equation takes over from the asymptotic formula. Taking the BMP closer to the SP will generally improve the accuracy of the asymptotic formula, but will make the computation more expensive as the Pruefer differential equations generally become progressively more ill- behaved as the SP is approached. The user is strongly recommended to experiment with placing the BMPs. In many singular problems quite crude asymptotic formulae will do. To help the user avoid needlessly accurate formulae, D02KEF outputs two 'sensitivity coefficients' (sigma) ,(sigma) which estimate how much the l r errors at the BMP's affect the computed eigenvalue. They are described in detail below, see Section 8.6. Example of coding BDYVAL: The example below illustrates typical situations: ( 2 ) y''+((lambda)-x- --)y=0, for 0=1. k+k' (b) The same index k can cause convergence to any of several eigenvalues depending on the initial values of ELAM and DELAM. It is up to the user to take suitable precautions - for instance by varying the position of the BMP's in the light of his knowledge of the asymptotic behaviour of the eigenfunction at different eigenvalues. 9. Example To find the 11th eigenvalue and eigenfunction of the example of Section 8.5, using the simple asymptotic formulae for the boundary conditions. Comparison of the results from this example program with the corresponding results from D02KDF(*) example program shows that similar output is produced from the routine MONIT, followed by the eigenfunction values from REPORT, and then a further line of information from MONIT (corresponding to the integration to find the eigenfunction). Final information is printed within the example program exactly as with D02KDF(*). 3 _ Note the discrepancy at the matching point c(=\/4 , the maximum of q(x;(lambda)), in this case) between the solutions obtained by integrations from left and right end-points. The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXd02raf}{NAG On-line Documentation: d02raf} \beginscroll \begin{verbatim} D02RAF(3NAG) Foundation Library (12/10/92) D02RAF(3NAG) D02 -- Ordinary Differential Equations D02RAF D02RAF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose D02RAF solves the two-point boundary-value problem with general boundary conditions for a system of ordinary differential equations, using a deferred correction technique and Newton iteration. 2. Specification SUBROUTINE D02RAF (N, MNP, NP, NUMBEG, NUMMIX, TOL, INIT, 1 X, Y, IY, ABT, FCN, G, IJAC, JACOBF, 2 JACOBG, DELEPS, JACEPS, JACGEP, WORK, 3 LWORK, IWORK, LIWORK, IFAIL) INTEGER N, MNP, NP, NUMBEG, NUMMIX, INIT, IY, 1 IJAC, LWORK, IWORK(LIWORK), LIWORK, IFAIL DOUBLE PRECISION TOL, X(MNP), Y(IY,MNP), ABT(N), DELEPS, 1 WORK(LWORK) EXTERNAL FCN, G, JACOBF, JACOBG, JACEPS, JACGEP 3. Description D02RAF solves a two-point boundary-value problem for a system of n ordinary differential equations in the interval (a,b) with b>a. The system is written in the form y '=f (x,y ,y ,...,y ) , i=1,2,...,n (1) i i 1 2 n and the derivatives f are evaluated by a subroutine FCN supplied i by the user. With the differential equations (1) must be given a system of n (nonlinear) boundary conditions g (y(a),y(b))=0 , i=1,2,...,n i where T y(x)=[y (x),y (x),...,y (x)] . (2) 1 2 n The functions g are evaluated by a subroutine G supplied by the i user. The solution is computed using a finite-difference technique with deferred correction allied to a Newton iteration to solve the finite-difference equations. The technique used is described fully in Pereyra [1]. The user must supply an absolute error tolerance and may also supply an initial mesh for the finite-difference equations and an initial approximate solution (alternatively a default mesh and approximation are used). The approximate solution is corrected using Newton iteration and deferred correction. Then, additional points are added to the mesh and the solution is recomputed with the aim of making the error everywhere less than the user's tolerance and of approximately equidistributing the error on the final mesh. The solution is returned on this final mesh. If the solution is required at a few specific points then these should be included in the initial mesh. If, on the other hand, the solution is required at several specific points then the user should use the interpolation routines provided in Chapter E01 if these points do not themselves form a convenient mesh. The Newton iteration requires Jacobian matrices ( ddf ) ( ddg ) ( ddg ) ( i) ( i ) ( i ) ( ----), ( -------) and ( -------). ( ddy ) ( ddy (a)) ( ddy (b)) ( j) ( j ) ( j ) These may be supplied by the user through subroutines JACOBF for ( ddf ) ( i) ( ----) and JACOBG for the others. Alternatively the Jacobians ( ddy ) ( j) may be calculated by numerical differentiation using the algorithm described in Curtis et al [2]. For problems of the type (1) and (2) for which it is difficult to determine an initial approximation from which the Newton iteration will converge, a continuation facility is provided. The user must set up a family of problems y'=f(x,y,(epsilon)), g(y(a),y(b),(epsilon))=0, (3) T where f=[f ,f ,...,f ] etc, and where (epsilon) is a 1 2 n continuation parameter. The choice (epsilon)=0 must give a problem (3) which is easy to solve and (epsilon)=1 must define the problem whose solution is actually required. The routine solves a sequence of problems with (epsilon) values 0=(epsilon) <(epsilon) <...<(epsilon) =1. (4) 1 2 p The number p and the values (epsilon) are chosen by the routine i so that each problem can be solved using the solution of its ddf predecessor as a starting approximation. Jacobians ----------- dd(epsilon) ddg and ----------- are required and they may be supplied by the dd(epsilon) user via routines JACEPS and JACGEP respectively or may be computed by numerical differentiation. 4. References [1] Pereyra V (1979) PASVA3: An Adaptive Finite-Difference Fortran Program for First Order Nonlinear, Ordinary Boundary Problems. Codes for Boundary Value Problems in Ordinary Differential Equations. Lecture Notes in Computer Science. (ed B Childs, M Scott, J W Daniel, E Denman and P Nelson) 76 Springer-Verlag. [2] Curtis A R, Powell M J D and Reid J K (1974) On the Estimation of Sparse Jacobian Matrices. J. Inst. Maths Applics. 13 117--119. 5. Parameters 1: N -- INTEGER Input On entry: the number of differential equations, n. Constraint: N > 0. 2: MNP -- INTEGER Input On entry: MNP must be set to the maximum permitted number of points in the finite-difference mesh. If LWORK or LIWORK (see below) is too small then internally MNP will be replaced by the maximum permitted by these values. (A warning message will be output if on entry IFAIL is set to obtain monitoring information.) Constraint: MNP >= 32. 3: NP -- INTEGER Input/Output On entry: NP must be set to the number of points to be used in the initial mesh. Constraint: 4 <= NP <= MNP. On exit: the number of points in the final mesh. 4: NUMBEG -- INTEGER Input On entry: the number of left-hand boundary conditions (that is the number involving y(a) only). Constraint: 0 <= NUMBEG < N. 5: NUMMIX -- INTEGER Input On entry: the number of coupled boundary conditions (that is the number involving both y(a) and y(b)). Constraint: 0 <= NUMMIX <= N - NUMBEG. 6: TOL -- DOUBLE PRECISION Input On entry: a positive absolute error tolerance. If a=x 0.0. 7: INIT -- INTEGER Input On entry: indicates whether the user wishes to supply an initial mesh and approximate solution (INIT /= 0) or whether default values are to be used, (INIT = 0). 8: X(MNP) -- DOUBLE PRECISION array Input/Output On entry: the user must set X(1) = a and X(NP) = b. If INIT = 0 on entry a default equispaced mesh will be used, otherwise the user must specify a mesh by setting X(i)=x , i for i = 2,3,...NP-1. Constraints: X(1) < X(NP), if INIT = 0, X(1) < X(2) <... < X(NP), if INIT /= 0. On exit: X(1),X(2),...,X(NP) define the final mesh (with the returned value of NP) and X(1) = a and X(NP) = b. 9: Y(IY,MNP) -- DOUBLE PRECISION array Input/Output On entry: if INIT = 0, then Y need not be set. If INIT /= 0, then the array Y must contain an initial approximation to the solution such that Y(j,i) contains an approximation to y (x ) , i=1,2,...,NP ; j=1,2,...,n. j i On exit: the approximate solution z (x ) satisfying (5) on j i the final mesh, that is Y(j,i)=z (x ) , i=1,2,...,NP ; j=1,2,...,n, j i where NP is the number of points in the final mesh. If an error has occurred then Y contains the latest approximation to the solution. The remaining columns of Y are not used. 10: IY -- INTEGER Input On entry: the first dimension of the array Y as declared in the (sub)program from which D02RAF is called. Constraint: IY >= N. 11: ABT(N) -- DOUBLE PRECISION array Output On exit: ABT(i), for i=1,2,...,n, holds the largest estimated error (in magnitude) of the ith component of the solution over all mesh points. 12: FCN -- SUBROUTINE, supplied by the user. External Procedure FCN must evaluate the functions f (i.e., the derivatives i y' ) at a general point x for a given value of (epsilon), i the continuation parameter (see Section 3). Its specification is: SUBROUTINE FCN (X, EPS, Y, F, N) INTEGER N DOUBLE PRECISION X, EPS, Y(N), F(N) 1: X -- DOUBLE PRECISION Input On entry: the value of the argument x. 2: EPS -- DOUBLE PRECISION Input On entry: the value of the continuation parameter, (epsilon). This is 1 if continuation is not being used. 3: Y(N) -- DOUBLE PRECISION array Input On entry: the value of the argument y , for i i=1,2,...,n. 4: F(N) -- DOUBLE PRECISION array Output On exit: the values of f , for i=1,2,...,n. i 5: N -- INTEGER Input On entry: the number of equations. FCN must be declared as EXTERNAL in the (sub)program from which D02RAF is called. Parameters denoted as Input must not be changed by this procedure. 13: G -- SUBROUTINE, supplied by the user. External Procedure G must evaluate the boundary conditions in equation (3) and place them in the array BC. Its specification is: SUBROUTINE G (EPS, YA, YB, BC, N) INTEGER N DOUBLE PRECISION EPS, YA(N), YB(N), BC(N) 1: EPS -- DOUBLE PRECISION Input On entry: the value of the continuation parameter, (epsilon). This is 1 if continuation is not being used. 2: YA(N) -- DOUBLE PRECISION array Input On entry: the value y (a), for i=1,2,...,n. i 3: YB(N) -- DOUBLE PRECISION array Input On entry: the value y (b), for i=1,2,...,n. i 4: BC(N) -- DOUBLE PRECISION array Output On exit: the values g (y(a),y(b),(epsilon)), for i i=1,2,...,n. These must be ordered as follows: (i) first, the conditions involving only y(a) (see NUMBEG description above); (ii) next, the NUMMIX coupled conditions involving both y(a) and y(b) (see NUMMIX description above); and, (iii) finally, the conditions involving only y(b) (N- NUMBEG-NUMMIX). 5: N -- INTEGER Input On entry: the number of equations, n. G must be declared as EXTERNAL in the (sub)program from which D02RAF is called. Parameters denoted as Input must not be changed by this procedure. 14: IJAC -- INTEGER Input On entry: indicates whether or not the user is supplying Jacobian evaluation routines. If IJAC /= 0 then the user must supply routines JACOBF and JACOBG and also, when continuation is used, routines JACEPS and JACGEP. If IJAC = 0 numerical differentiation is used to calculate the Jacobian and the routines D02GAZ, D02GAY, D02GAZ and D02GAX ( ddf ) ( i) JACOBF must evaluate the Jacobian ( ----) for i,j=1,2,...,n, ( ddy ) ( j) given x and y , for j=1,2,...,n. j Its specification is: SUBROUTINE JACOBF (X, EPS, Y, F, N) INTEGER N DOUBLE PRECISION X, EPS, Y(N), F(N,N) 1: X -- DOUBLE PRECISION Input On entry: the value of the argument x. 2: EPS -- DOUBLE PRECISION Input On entry: the value of the continuation parameter (epsilon). This is 1 if continuation is not being used. 3: Y(N) -- DOUBLE PRECISION array Input On entry: the value of the argument y , for i i=1,2,...,n. 4: F(N,N) -- DOUBLE PRECISION array Output ddf i On exit: F(i,j) must be set to the value of ----, ddy j evaluated at the point (x,y), for i,j=1,2,...,n. 5: N -- INTEGER Input On entry: the number of equations, n. JACOBF must be declared as EXTERNAL in the (sub)program from which D02RAF is called. Parameters denoted as Input must not be changed by this procedure. 16: JACOBG -- SUBROUTINE, supplied by the user. External Procedure ( ddg ) ( ddg ) ( i ) ( i ) JACOBG must evaluate the Jacobians ( -------) and ( -------) ( ddy (a)) ( ddy (b)) ( j ) ( j ) The ordering of the rows of AJ and BJ must correspond to the ordering of the boundary conditions described in the specification of subroutine G above. Its specification is: SUBROUTINE JACOBG (EPS, YA, YB, AJ, BJ, N) INTEGER N DOUBLE PRECISION EPS, YA(N), YB(N), AJ(N,N), BJ 1 (N,N) 1: EPS -- DOUBLE PRECISION Input On entry: the value of the continuation parameter, (epsilon). This is 1 if continuation is not being used. 2: YA(N) -- DOUBLE PRECISION array Input On entry: the value y (a), for i=1,2,...,n. i 3: YB(N) -- DOUBLE PRECISION array Input On entry: the value y (b), for i=1,2,...,n. i 4: AJ(N,N) -- DOUBLE PRECISION array Output ddg i On exit: AJ(i,j) must be set to the value -------, ddy (a) j for i,j=1,2,...,n. 5: BJ(N,N) -- DOUBLE PRECISION array Output ddg i On exit: BJ(i,j) must be set to the value -------, ddy (b) j for i,j=1,2...,n. 6: N -- INTEGER Input On entry: the number of equations, n. JACOBG must be declared as EXTERNAL in the (sub)program from which D02RAF is called. Parameters denoted as Input must not be changed by this procedure. 17: DELEPS -- DOUBLE PRECISION Input/Output On entry: DELEPS must be given a value which specifies whether continuation is required. If DELEPS <= 0.0 or DELEPS >= 1.0 then it is assumed that continuation is not required. If 0.0 < DELEPS < 1.0 then it is assumed that continuation is required unless DELEPS < square root of machine precision when an error exit is taken. DELEPS is used as the increment (epsilon) -(epsilon) (see (4)) and the choice DELEPS = 0.1 2 1 is recommended. On exit: an overestimate of the increment (epsilon) -(epsilon) (in fact the value of the increment p p-1 which would have been tried if the restriction (epsilon) =1 p had not been imposed). If continuation was not requested then DELEPS = 0.0. If continuation is not requested then the parameters JACEPS and JACGEP may be replaced by dummy actual parameters in the call to D02RAF. (D02GAZ and D02GAX respectively may be used as the dummy parameters.) 18: JACEPS -- SUBROUTINE, supplied by the user. External Procedure ddf i JACEPS must evaluate the derivative ----------- given x and dd(epsilon) y if continuation is being used. Its specification is: SUBROUTINE JACEPS (X, EPS, Y, F, N) INTEGER N DOUBLE PRECISION X, EPS, Y(N), F(N) 1: X -- DOUBLE PRECISION Input On entry: the value of the argument x. 2: EPS -- DOUBLE PRECISION Input On entry: the value of the continuation parameter, (epsilon). 3: Y(N) -- DOUBLE PRECISION array Input On entry: the solution values y at the point x, for i i=1,2,...,n. 4: F(N) -- DOUBLE PRECISION array Output ddf i On exit: F(i) must contain the value ----------- at dd(epsilon) the point (x,y), for i=1,2,...,n. 5: N -- INTEGER Input On entry: the number of equations, n. JACEPS must be declared as EXTERNAL in the (sub)program from which D02RAF is called. Parameters denoted as Input must not be changed by this procedure. 19: JACGEP -- SUBROUTINE, supplied by the user. External Procedure ddg i JACGEP must evaluate the derivatives ----------- if dd(epsilon) continuation is being used. Its specification is: SUBROUTINE JACGEP (EPS, YA, YB, BCEP, N) INTEGER N DOUBLE PRECISION EPS, YA(N), YB(N), BCEP(N) 1: EPS -- DOUBLE PRECISION Input On entry: the value of the continuation parameter, (epsilon). 2: YA(N) -- DOUBLE PRECISION array Input On entry: the value of y (a), for i=1,2,...,n. i 3: YB(N) -- DOUBLE PRECISION array Input On entry: the value of y (b), for i=1,2,...,n. i 4: BCEP(N) -- DOUBLE PRECISION array Output On exit: BCEP(i) must contain the value of ddg i -----------, for i=1,2,...,n. dd(epsilon) 5: N -- INTEGER Input On entry: the number of equations, n. JACGEP must be declared as EXTERNAL in the (sub)program from which D02RAF is called. Parameters denoted as Input must not be changed by this procedure. 20: WORK(LWORK) -- DOUBLE PRECISION array Workspace 21: LWORK -- INTEGER Input On entry: the dimension of the array WORK as declared in the (sub)program from which D02RAF is called. 2 2 Constraint: LWORK>=MNP*(3N +6N+2)+4N +3N. 22: IWORK(LIWORK) -- INTEGER array Workspace 23: LIWORK -- INTEGER Input On entry: the dimension of the array IWORK as declared in the (sub)program from which D02RAF is called. Constraints: LIWORK>=MNP*(2*N+1)+N, if IJAC /= 0, 2 LIWORK>=MNP*(2*N+1)+N +4*N+2, if IJAC = 0. 24: IFAIL -- INTEGER Input/Output For this routine, the normal use of IFAIL is extended to control the printing of error and warning messages as well as specifying hard or soft failure (see the Essential Introduction). Before entry, IFAIL must be set to a value with the decimal expansion cba, where each of the decimal digits c, b and a must have a value of 0 or 1. a=0 specifies hard failure, otherwise soft failure; b=0 suppresses error messages, otherwise error messages will be printed (see Section 6); c=0 suppresses warning messages, otherwise warning messages will be printed (see Section 6). The recommended value for inexperienced users is 110 (i.e., hard failure with all messages printed). Unless the routine detects an error (see Section 6), IFAIL contains 0 on exit. 6. Error Indicators and Warnings Errors detected by the routine: For each error, an explanatory error message is output on the current error message unit (as defined by X04AAF), unless suppressed by the value of IFAIL on entry. IFAIL= 1 One or more of the parameters N, MNP, NP, NUMBEG, NUMMIX, TOL, DELEPS, LWORK or LIWORK has been incorrectly set, or X (1) >= X(NP) or the mesh points X(i) are not in strictly ascending order. IFAIL= 2 A finer mesh is required for the accuracy requested; that is MNP is not large enough. This error exit normally occurs when the problem being solved is difficult (for example, there is a boundary layer) and high accuracy is requested. A poor initial choice of mesh points will make this error exit more likely. IFAIL= 3 The Newton iteration has failed to converge. There are several possible causes for this error: (i) faulty coding in one of the Jacobian calculation routines; (ii) if IJAC = 0 then inaccurate Jacobians may have been calculated numerically (this is a very unlikely cause); or, (iii) a poor initial mesh or initial approximate solution has been selected either by the user or by default or there are not enough points in the initial mesh. Possibly, the user should try the continuation facility. IFAIL= 4 The Newton iteration has reached round-off error level. It could be however that the answer returned is satisfactory. The error is likely to occur if too high an accuracy is requested. IFAIL= 5 The Jacobian calculated by JACOBG (or the equivalent matrix calculated by numerical differentiation) is singular. This may occur due to faulty coding of JACOBG or, in some circumstances, to a zero initial choice of approximate solution (such as is chosen when INIT = 0). IFAIL= 6 There is no dependence on (epsilon) when continuation is being used. This can be due to faulty coding of JACEPS or JACGEP or, in some circumstances, to a zero initial choice of approximate solution (such as is chosen when INIT = 0). IFAIL= 7 DELEPS is required to be less than machine precision for continuation to proceed. It is likely that either the problem (3) has no solution for some value near the current value of (epsilon) (see the advisory print out from D02RAF) or that the problem is so difficult that even with continuation it is unlikely to be solved using this routine. If the latter cause is suspected then using more mesh points initially may help. IFAIL= 8 Indicates that a serious error has occurred in a call to D02RAF. Check all array subscripts and subroutine parameter lists in calls to D02RAF. Seek expert help. IFAIL= 9 Indicates that a serious error has occurred in a call to D02RAR. Check all array subscripts and subroutine parameter lists in calls to D02RAF. Seek expert help. 7. Accuracy The solution returned by the routine will be accurate to the user's tolerance as defined by the relation (5) except in extreme circumstances. The final error estimate over the whole mesh for each component is given in the array ABT. If too many points are specified in the initial mesh, the solution may be more accurate than requested and the error may not be approximately equidistributed. 8. Further Comments There are too many factors present to quantify the timing. The time taken by the routine is negligible only on very simple problems. The user is strongly recommended to set IFAIL to obtain self- explanatory error messages, and also monitoring information about the course of the computation. In the case where the user wishes to solve a sequence of similar problems, the use of the final mesh and solution from one case as the initial mesh is strongly recommended for the next. 9. Example We solve the differential equation 2 y'''=-yy''-2(epsilon)(1-y' ) with (epsilon)=1 and boundary conditions y(0)=y'(0)=0, y'(10)=1 to an accuracy specified by TOL=1.0E-4. The continuation facility is used with the continuation parameter (epsilon) introduced as in the differential equation above and with DELEPS = 0.1 initially. (The continuation facility is not needed for this problem and is used here for illustration.) The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXd03}{NAG On-line Documentation: d03} \beginscroll \begin{verbatim} D03(3NAG) Foundation Library (12/10/92) D03(3NAG) D03 -- Partial Differential Equations Introduction -- D03 Chapter D03 Partial Differential Equations 1. Scope of the Chapter This chapter is concerned with the solution of partial differential equations. 2. Background to the Problems The definition of a partial differential equation problem includes not only the equation itself but also the domain of interest and appropriate subsidiary conditions. Indeed, partial differential equations are usually classified as elliptic, hyperbolic or parabolic according to the form of the equation and the form of the subsidiary conditions which must be assigned to produce a well-posed problem. Ultimately it is hoped that this chapter will contain routines for the solution of equations of each of these types together with automatic mesh generation routines and other utility routines particular to the solution of partial differential equations. The routines in this chapter will often call upon routines from the Linear Algebra Chapter F04 -- Simultaneous Linear Equations. The classification of partial differential equations is easily described in the case of linear equations of the second order in two independent variables, i.e., equations of the form au +2bu +cu +du +eu +fu+g=0, (1) xx xy yy x y where a, b, c, d, e, f and g are functions of x and y only. Equation (1) is called elliptic, hyperbolic or parabolic 2 according as ac-b is positive, negative or zero. Useful definitions of the concepts of elliptic, hyperbolic and parabolic character can also be given for differential equations in more than two independent variables, for systems and for nonlinear differential equations. For elliptic equations, of which Laplace's equation u +u =0 (2) xx yy is the simplest example of second order, the subsidiary conditions take the form of boundary conditions, i.e., conditions which provide information about the solution at all points of a closed boundary. For example, if equation (2) holds in a plane domain D bounded by a contour C, a solution u may be sought subject to the condition u=f on C, (3) where f is a given function. The condition (3) is known as a Dirichlet boundary condition. Equally common is the Neumann boundary condition u'=g on C, (4) which is one form of a more general condition u'+fu=g on C, (5) where u' denotes the derivative of u normal to the contour C and f and g are given functions. Provided that f and g satisfy certain restrictions, condition (5) yields a well-posed boundary value problem for Laplace's equation. In the case of the Neumann problem, one further piece of information, e.g. the value of u at a particular point, is necessary for uniqueness of the solution. Boundary conditions similar to the above are applicable to more general second order elliptic equations, whilst two such conditions are required for equations of fourth order. For hyperbolic equations, the wave equation u -u =0 (6) tt xx is the simplest example of second order. It is equivalent to a first order system u -v =0, v -u =0. (7) t x t x The subsidiary conditions may take the form of initial conditions, i.e., conditions which provide information about the solution at points on a suitable open boundary. For example, if equation (6) is satisfied for t>0, a solution u may be sought such that u(x,0)=f(x), u (x,0)=g(x), (8) t where f and g are given functions. This is an example of an initial value problem, sometimes known as Cauchy's problem. For parabolic equations, of which the heat conduction equation u -u =0 (9) t xx is the simplest example, the subsidiary conditions always include some of initial type and may also include some of boundary type. For example, if equation (9) is satisfied for t>0 and 00. (11) This is an example of a mixed initial/boundary value problem. For all types of partial differential equations, finite difference methods (Mitchell and Griffiths [5]) and finite element methods (Wait and Mitchell [9]) are the most common means of solution and such methods obviously feature prominently either in this chapter or in the companion NAG Finite Element Library. Many of the utility routines in this chapter are concerned with the solution of the large sparse systems of equations which arise from the finite difference and finite element methods. Alternative methods of solution are often suitable for special classes of problems. For example, the method of characteristics is the most common for hyperbolic equations involving time and one space dimension (Smith [7]). The method of lines (see Mikhlin and Smolitsky [4]) may be used to reduce a parabolic equation to a (stiff) system of ordinary differential equations, which may be solved by means of routines from Chapter D02 -- Ordinary Differential Equations. Similarly, integral equation or boundary element methods (Jaswon and Symm [3]) are frequently used for elliptic equations. Typically, in the latter case, the solution of a boundary value problem is represented in terms of certain boundary functions by an integral expression which satisfies the differential equation throughout the relevant domain. The boundary functions are obtained by applying the given boundary conditions to this representation. Implementation of this method necessitates discretisation of only the boundary of the domain, the dimensionality of the problem thus being effectively reduced by one. The boundary conditions yield a full system of simultaneous equations, as opposed to the sparse systems yielded by the finite difference and finite element methods, but the full system is usually of much lower order. Solution of this system yields the boundary functions, from which the solution of the problem may be obtained, by quadrature, as and where required. 2.1. References [1] Ames W F (1977) Nonlinear Partial Differential Equations in Engineering. Academic Press (2nd Edition). [2] Berzins M (1990) Developments in the NAG Library Software for Parabolic Equations. Scientific Software Systems. (ed J C Mason and M G Cox) Chapman and Hall. 59--72. [3] Jaswon M A and Symm G T (1977) Integral Equation Methods in Potential Theory and Elastostatics. Academic Press. [4] Mikhlin S G and Smolitsky K L (1967) Approximate Methods for the Solution of Differential and Integral Equations. Elsevier. [5] Mitchell A R and Griffiths D F (1980) The Finite Difference Method in Partial Differential Equations. Wiley. [6] Richtmyer R D and Morton K W (1967) Difference Methods for Initial-value Problems. Interscience (2nd Edition). [7] Smith G D (1985) Numerical Solution of Partial Differential Equations: Finite Difference Methods. Oxford University Press (3rd Edition). [8] Swarztrauber P N and Sweet R A (1979) Efficient Fortran Subprograms for the Solution of Separable Elliptic Partial Differential Equations. ACM Trans. Math. Softw. 5 352--364. [9] Wait R and Mitchell A R (1985) Finite Element Analysis and Application. Wiley. 3. Recommendations on Choice and Use of Routines The choice of routine will depend first of all upon the type of partial differential equation to be solved. At present no special allowances are made for problems with boundary singularities such as may arise at corners of domains or at points where boundary conditions change. For such problems results should be treated with caution. Users may wish to construct their own partial differential equation solution software for problems not solvable by the routines described in Sections 3.1 to 3.4 below. In such cases users can employ appropriate routines from the Linear Algebra Chapters to solve the resulting linear systems; see Section 3.5 for further details. 3.1. Elliptic Equations The routine D03EDF solves a system of seven-point difference equations in a rectangular grid (in two dimensions), using the multigrid iterative method. The equations are supplied by the user, and the seven-point form allows cross-derivative terms to be represented (see Mitchell and Griffiths [5]). The method is particularly efficient for large systems of equations with diagonal dominance. The routine D03EEF discretises a second-order equation on a two- dimensional rectangular region using finite differences and a seven-point molecule. The routine allows for cross-derivative terms, Dirichlet, Neumann or mixed boundary conditions, and either central or upwind differences. The resulting seven- diagonal difference equations are in a form suitable for passing directly to the multigrid routine D03EDF, although other solution methods could easily be used. The routine D03FAF, based on the routine HW3CRT from FISHPACK (Swarztrauber and Sweet [8]), solves the Helmholtz equation in a three-dimensional cuboidal region, with any combination of Dirichlet, Neumann or periodic boundary conditions. The method used is based on the fast Fourier transform algorithm, and is likely to be particularly efficient on vector-processing machines. 3.2. Hyperbolic Equations There are no routines available yet for the solution of these equations. 3.3. Parabolic Equations There are no routines available yet for the solution of these equations. But problems in two space dimensions plus time may be treated as a succession of elliptic equations [1], [6] using appropriate D03E- routines or one may use codes from the NAG Finite Element Library. 3.4. Utility Routines There are no utility routines available yet, but routines are available in the Linear Algebra Chapters for the direct and iterative solution of linear equations. Here we point to some of the routines that may be of use in solving the linear systems that arise from finite difference or finite element approximations to partial differential equation solutions. Chapters F01 and F04 should be consulted for further information and for the routine documents. Decision trees for the solution of linear systems are given in Section 3.5 of the F04 Chapter Introduction. The following routines allow the direct solution of symmetric positive-definite systems: Band F04ACF Variable band F01MCF and F04MCF (skyline) Tridiagonal F04FAF Sparse F01MAF* and F04MAF (* the parameter DROPTL should be set to zero for F01MAF to give a direct solution) and the following routines allow the iterative solution of symmetric positive-definite systems: Sparse (incomplete F01MAF and F04MBF Cholesky) Sparse (conjugate F04MBF gradient) The latter routine above allows the user to supply a pre- conditioner and also allows the solution of indefinite symmetric systems. The following routines allow the direct solution of unsymmetric systems: Band F01LBF and F04LDF Almost block- F01LHF and F04LHF diagonal Tridiagonal F01LEF and F04LEF or F04EAF Sparse F01BRF (and F01BSF) and F04AXF and the following routine allows the iterative solution of unsymmetric systems: Sparse F04QAF The above routine allows the user to supply a pre-conditioner and also allows the solution of least-squares systems. 3.5. Index Elliptic equations equations on rectangular grid (seven-point 2-D D03EDF molecule) discretisation on rectangular grid (seven-point 2-D D03EEF molecule) Helmholtz's equation in three dimensions D03FAF D03 -- Partial Differential Equations Contents -- D03 Chapter D03 Partial Differential Equations D03EDF Elliptic PDE, solution of finite difference equations by a multigrid technique D03EEF Discretize a 2nd order elliptic PDE on a rectangle D03FAF Elliptic PDE, Helmholtz equation, 3-D Cartesian co- ordinates \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXd03edf}{NAG On-line Documentation: d03edf} \beginscroll \begin{verbatim} D03EDF(3NAG) Foundation Library (12/10/92) D03EDF(3NAG) D03 -- Partial Differential Equations D03EDF D03EDF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose D03EDF solves seven-diagonal systems of linear equations which arise from the discretization of an elliptic partial differential equation on a rectangular region. This routine uses a multigrid technique. 2. Specification SUBROUTINE D03EDF (NGX, NGY, LDA, A, RHS, UB, MAXIT, ACC, 1 US, U, IOUT, NUMIT, IFAIL) INTEGER NGX, NGY, LDA, MAXIT, IOUT, NUMIT, IFAIL DOUBLE PRECISION A(LDA,7), RHS(LDA), UB(NGX*NGY), ACC, US 1 (LDA), U(LDA) 3. Description D03EDF solves, by multigrid iteration, the seven-point scheme 6 7 A u +A u i,j i-1,j+1 i,j i,j+1 3 4 5 +A u +A u +A u i,j i-1,j i,j ij i,j i+1,j 1 2 +A u +A u =f , i,j i,j-1 i,j i+1,j-1 ij i=1,2,...,n ; j=1,2,...,n , x y which arises from the discretization of an elliptic partial differential equation of the form (alpha)(x,y)U +(beta)(x,y)U +(gamma)(x,y)U +(delta)(x,y)U xx xy yy x +(epsilon)(x,y)U +(phi)(x,y)U=(psi)(x,y) y and its boundary conditions, defined on a rectangular region. This we write in matrix form as Au=f The algorithm is described in separate reports by Wesseling [2], [3] and McCarthy [1]. Systems of linear equations, matching the seven-point stencil defined above, are solved by a multigrid iteration. An initial estimate of the solution must be provided by the user. A zero guess may be supplied if no better approximation is available. A 'smoother' based on incomplete Crout decomposition is used to eliminate the high frequency components of the error. A restriction operator is then used to map the system on to a sequence of coarser grids. The errors are then smoothed and prolongated (mapped onto successively finer grids). When the finest cycle is reached, the approximation to the solution is corrected. The cycle is repeated for MAXIT iterations or until the required accuracy, ACC, is reached. D03EDF will automatically determine the number l of possible coarse grids, 'levels' of the multigrid scheme, for a particular problem. In other words, D03EDF determines the maximum integer l so that n and n can be expressed in the form x y l-1 l-1 n =m2 +1, n =n2 +1, with m>=2 and n>=2. x y It should be noted that the rate of convergence improves significantly with the number of levels used (see McCarthy [1]), so that n and n should be carefully chosen so that n -1 and x y x l n -1 have factors of the form 2 , with l as large as possible. y For good convergence the integer l should be at least 2. D03EDF has been found to be robust in application, but being an iterative method the problem of divergence can arise. For a strictly diagonally dominant matrix A 4 k ij -- ij |A |> > |A | -- k/=4 no such problem is foreseen. The diagonal dominance of A is not a necessary condition, but should this condition be strongly violated then divergence may occur. The quickest test is to try the routine. 4. References [1] McCarthy G J (1983) Investigation into the Multigrid Code MGD1. Report AERE-R 10889. Harwell. [2] Wesseling P (1982) MGD1 - A Robust and Efficient Multigrid Method. Multigrid Methods. Lecture Notes in Mathematics. 960 Springer-Verlag. 614--630. [3] Wesseling P (1982) Theoretical Aspects of a Multigrid Method. SIAM J. Sci. Statist. Comput. 3 387--407. 5. Parameters 1: NGX -- INTEGER Input On entry: the number of interior grid points in the x- direction, n . NGX-1 should preferably be divisible by as x high a power of 2 as possible. Constraint: NGX >= 3. 2: NGY -- INTEGER Input On entry: the number of interior grid points in the y- direction, n . NGY-1 should preferably be divisible by as y high a power of 2 as possible. Constraint: NGY >= 3. 3: LDA -- INTEGER Input On entry: the first dimension of the array A as declared in the (sub)program from which D03EDF is called, which must also be a lower bound for the dimensions of the arrays RHS, US and U. It is always sufficient to set LDA>=(4*(NGX+1)*(NGY+1))/3, but slightly smaller values may be permitted, depending on the values of NGX and NGY. If on entry, LDA is too small, an error message gives the minimum permitted value. (LDA must be large enough to allow space for the coarse-grid approximations). 4: A(LDA,7) -- DOUBLE PRECISION array Input/Output k On entry: A(i+(j-1)*NGX,k) must be set to A , for i = ij 1,2,...,NGX; j = 1,2,...,NGY and k = 1,2,...,7. On exit: A is overwritten. 5: RHS(LDA) -- DOUBLE PRECISION array Input/Output On entry: RHS(i+(j-1)*NGX) must be set to f , for i = ij 1,2,...,NGX; j = 1,2,...,NGY. On exit: the first NGX*NGY elements are unchanged and the rest of the array is used as workspace. 6: UB(NGX*NGY) -- DOUBLE PRECISION array Input/Output On entry: UB(i+(j-1)*NGX) must be set to the initial estimate for the solution u . On exit: the corresponding ij component of the residual r=f-Au. 7: MAXIT -- INTEGER Input On entry: the maximum permitted number of multigrid iterations. If MAXIT = 0, no multigrid iterations are performed, but the coarse-grid approximations and incomplete Crout decompositions are computed, and may be output if IOUT is set accordingly. Constraint: MAXIT >= 0. 8: ACC -- DOUBLE PRECISION Input On entry: the required tolerance for convergence of the residual 2-norm: / NGX*NGY / -- 2 ||r|| = / > (r ) 2 / -- k \/ k=1 where r=f-Au and u is the computed solution. Note that the norm is not scaled by the number of equations. The routine will stop after fewer than MAXIT iterations if the residual 2-norm is less than the specified tolerance. (If MAXIT > 0, at least one iteration is always performed.) If on entry ACC = 0.0, then the machine precision is used as a default value for the tolerance; if ACC > 0.0, but ACC is less than the machine precision, then the routine will stop when the residual 2-norm is less than the machine precision and IFAIL will be set to 4. Constraint: ACC >= 0.0. 9: US(LDA) -- DOUBLE PRECISION array Output On exit: the residual 2-norm, stored in element US(1). 10: U(LDA) -- DOUBLE PRECISION array Output On exit: the computed solution u is returned in U(i+(j- ij 1)*NGX), for i = 1,2,...,NGX; j = 1,2,...,NGY. 11: IOUT -- INTEGER Input On entry: controls the output of printed information to the advisory message unit as returned by X04ABF: IOUT = 0 No output. IOUT = 1 The solution u , for i = 1,2,...,NGX; j = 1,2,... ij ,NGY. IOUT = 2 The residual 2-norm after each iteration, with the reduction factor over the previous iteration. IOUT = 3 As for IOUT = 1 and IOUT = 2. IOUT = 4 As for IOUT = 3, plus the final residual (as returned in UB). IOUT = 5 As for IOUT = 4, plus the initial elements of A and RHS. IOUT = 6 As for IOUT = 5, plus the Galerkin coarse grid approximations. IOUT = 7 As for IOUT = 6, plus the incomplete Crout decompositions. IOUT = 8 As for IOUT = 7, plus the residual after each iteration. The elements A(p,k), the Galerkin coarse grid approximations and the incomplete Crout decompositions are output in the format: Y-index = j X-index = i A(p,1) A(p,2) A(p,3) A(p,4) A(p,5) A(p,6) A(p,7) where p=i+(j-1)*NGX, i = 1,2,...,NGX and j = 1,2,... ,NGY. The vectors U(p), UB(p), RHS(p) are output in matrix form with NGY rows and NGX columns. Where NGX > 10, the NGX values for a given j-value are produced in rows of 10. Values of IOUT > 4 may yield considerable amounts of output. Constraint: 0 <= IOUT <= 8. 12: NUMIT -- INTEGER Output On exit: the number of iterations performed. 13: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. On exit: IFAIL = 0 unless the routine detects an error (see Section 6). 6. Error Indicators and Warnings Errors detected by the routine: If on entry IFAIL = 0 or -1, explanatory error messages are output on the current error message unit (as defined by X04AAF). IFAIL= 1 On entry NGX < 3, or NGY < 3, or LDA is too small, or ACC < 0.0, or MAXIT < 0, or IOUT < 0, or IOUT > 8. IFAIL= 2 MAXIT iterations have been performed with the residual 2- norm decreasing at each iteration but the residual 2-norm has not been reduced to less than the specified tolerance (see ACC). Examine the progress of the iteration by setting IOUT >= 2. IFAIL= 3 As for IFAIL = 2, except that at one or more iterations the residual 2-norm did not decrease. It is likely that the method fails to converge for the given matrix A. IFAIL= 4 On entry ACC is less than the machine precision. The routine terminated because the residual norm is less than the machine precision. 7. Accuracy See ACC (Section 5). 8. Further Comments The rate of convergence of this routine is strongly dependent upon the number of levels, l, in the multigrid scheme, and thus the choice of NGX and NGY is very important. The user is advised to experiment with different values of NGX and NGY to see the effect they have on the rate of convergence; for example, using a 6 value such as NGX = 65 (=2 +1) followed by NGX = 64 (for which l = 1). 9. Example The program solves the elliptic partial differential equation U -(alpha)U +U =-4, (alpha)=1.7 xx xy yy on the unit square 0<=x,y<=1, with boundary conditions {x=0, (0<=y<=1) U=0 on {y=0, (0<=x<=1) U=1 on x=1, 0<=y<=1. {y=1, (0<=x<=1) For the equation to be elliptic, (alpha) must be less than 2. The equation is discretized on a square grid with mesh spacing h in both directions using the following approximations: Please see figure in printed Reference Manual 1 U ~= --(U -2U +U ) xx 2 E O W h 1 U ~= --(U -2U +U ) yy 2 N O S h 1 U ~= ---(U -U +U -2U +U -U +U ). xy 2 N NW E O W SE S 2h Thus the following equations are solved: 1 1 -(alpha)u + (1- -(alpha))u 2 i-1,j+1 2 i,j+1 1 1 +(1- -(alpha))u + (-4+(alpha))u + (1- -(alpha))u 2 i-1,j ij 2 i+1,j 1 1 + (1- -(alpha))u + -(alpha)u 2 i,j-1 2 i+1,j-1 2 =-4h The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXd03eef}{NAG On-line Documentation: d03eef} \beginscroll \begin{verbatim} D03EEF(3NAG) Foundation Library (12/10/92) D03EEF(3NAG) D03 -- Partial Differential Equations D03EEF D03EEF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose D03EEF discretizes a second order elliptic partial differential equation (PDE) on a rectangular region. 2. Specification SUBROUTINE D03EEF (XMIN, XMAX, YMIN, YMAX, PDEF, BNDY, 1 NGX, NGY, LDA, A, RHS, SCHEME, IFAIL) INTEGER NGX, NGY, LDA, IFAIL DOUBLE PRECISION XMIN, XMAX, YMIN, YMAX, A(LDA,7), RHS(LDA) CHARACTER*1 SCHEME EXTERNAL PDEF, BNDY 3. Description D03EEF discretizes a second order linear elliptic partial differential equation of the form 2 2 2 dd U dd U dd U (alpha)(x,y) ----+(beta)(x,y) ------+(gamma)(x,y) ---- 2 ddxddy 2 dd x dd y ddU ddU +(delta)(x,y) ---+(epsilon)(x,y) ---+(phi)(x,y)U=(psi)(x,y) (1) ddx ddy on a rectangular region x <=x<=x A B y <=y<=y A B subject to boundary conditions of the form ddU a(x,y)U+b(x,y) ---=c(x,y) ddn ddU where --- denotes the outward pointing normal derivative on the ddn boundary. Equation (1) is said to be elliptic if 2 4(alpha)(x,y)(gamma)(x,y)>=((beta)(x,y)) for all points in the rectangular region. The linear equations produced are in a form suitable for passing directly to the multigrid routine D03EDF. The equation is discretized on a rectangular grid, with n grid x points in the x-direction and n grid points in the y-direction. y The grid spacing used is therefore h =(x -x )/(n -1) x B A x h =(y -y )/(n -1) y B A y and the co-ordinates of the grid points (x ,y ) are i j x =x +(i-1)h , i=1,2,...,n , i A x x y =y +(j-1)h , j=1,2,...,n . j A y y At each grid point (x ,y ) six neighbouring grid points are used i j to approximate the partial differential equation, so that the equation is discretized on the following seven-point stencil: Please see figure in printed Reference Manual For convenience the approximation u to the exact solution ij U(x ,y ) is denoted by u , and the neighbouring approximations i j O are labelled according to points of the compass as shown. Where numerical labels for the seven points are required, these are also shown above. The following approximations are used for the second derivatives: 2 dd U 1 ----~= --(u -2u +u ) 2 2 E O W dd x h x 2 dd U 1 ----~= --(u -2u +u ) 2 2 N O S dd y h y 2 dd U 1 ------~= -----(u -u +u -2u +u -u +u ). ddxddy 2h h N NW E O W SE S x y Two possible schemes may be used to approximate the first derivatives: Central Differences ddU 1 ---~= ---(u -u ) ddx 2h W E x ddU 1 ---~= ---(u -u ) ddy 2h N S y Upwind Differences ddU 1 ---~= --(u -u )if (delta)(x,y)>0 ddx h E O x ddU 1 ---~= --(u -u )if (delta)(x,y)<0 ddx h O W x ddU 1 ---~= --(u -u )if (epsilon)(x,y)>0 ddy h N O y ddU 1 ---~= --(u -u )if (epsilon)(x,y)<0. ddy h O S y Central differences are more accurate than upwind differences, but upwind differences may lead to a more diagonally dominant matrix for those problems where the coefficients of the first derivatives are significantly larger than the coefficients of the second derivatives. The approximations used for the first derivatives may be written in a more compact form as follows: ddU 1 ---~= ---((k -1)u -2k u +(k +1)u ) ddx 2h ( x W x O x E) x ddU 1 ---~= ---((k -1)u -2k u +(k +1)u ) ddy 2h ( y S y O y N) y where k =sign (delta) and k =sign (epsilon) for upwind x y differences, and k =k =0 for central differences. x y At all points in the rectangular domain, including the boundary, the coefficients in the partial differential equation are evaluated by calling the user-supplied subroutine PDEF, and applying the approximations. This leads to a seven-diagonal system of linear equations of the form: 6 7 A u +A u ij i-1,j+1 ij i,j+1 3 4 5 +A u +A u +A u ij i-1,j ij ij ij i+1,j 1 2 +A u +A u =f , i=1,2,...,n ;j=1,2,...,n , ij i,j-1 ij i+1,j-1 ij x y where the coefficients are given by 1 1 1 1 A =(beta)(x ,y )----- +(gamma)(x ,y )-- +(epsilon)(x ,y )---(k -1) ij i j 2h h i j 2 i j 2h y x y h y y 2 1 A =-(beta)(x ,y ) ----- ij i j 2h h x y 3 1 1 1 A =(alpha)(x ,y ) --+(beta)(x ,y ) -----+(delta)(x ,y ) ---(k -1) ij i j 2 i j 2h h i j 2h x h x y x x 4 2 1 2 A =-(alpha)(x ,y ) -- -(beta)(x ,y ) ---- -(gamma)(x ,y ) -- ij i j 2 i j h h i j 2 h x y h x y k k y y -(delta)(x ,y ) -- -(epsilon)(x ,y ) -- -(phi)(x ,y ) i j h i j h i j x y 5 1 1 1 A =(alpha)(x ,y ) --+(beta)(x ,y ) -----+(delta)(x ,y ) ---(k +1) ij i j 2 i j 2h h i j 2h x h x y x x 6 1 A =-(beta)(x ,y ) ----- ij i j 2h h x y 7 1 1 1 A =(beta)(x ,y )-----+(gamma)(x ,y )--+(epsilon)(x ,y )---(k +1) ij i j 2h h i j 2 i j 2h y x y h y y f =(psi)(x ,y ) ij i j These equations then have to be modified to take account of the boundary conditions. These may be Dirichlet (where the solution is given), Neumann (where the derivative of the solution is given), or mixed (where a linear combination of solution and derivative is given). If the boundary conditions are Dirichlet, there are an infinity of possible equations which may be applied: (mu)u =(mu)f , (mu)/=0. (2) ij ij If D03EDF is used to solve the discretized equations, it turns out that the choice of (mu) can have a dramatic effect on the rate of convergence, and the obvious choice (mu)=1 is not the best. Some choices may even cause the multigrid method to fail altogether. In practice it has been found that a value of the same order as the other diagonal elements of the matrix is best, and the following value has been found to work well in practice: ( { 2 2 } 4 ) (mu)=min (-{ --+ --},A ). ij( { 2 2} ij) ( { h h } ) ( { x y} ) If the boundary conditions are either mixed or Neumann (i.e., B /= 0 on return from the user-supplied subroutine BNDY), then one of the points in the seven-point stencil lies outside the domain. In this case the normal derivative in the boundary conditions is used to eliminate the 'fictitious' point, u : outside ddU 1 ---~= --(u -u ). (3) ddn 2h outside inside It should be noted that if the boundary conditions are Neumann and (phi)(x,y)==0, then there is no unique solution. The routine returns with IFAIL = 5 in this case, and the seven-diagonal matrix is singular. The four corners are treated separately. The user-supplied subroutine BNDY is called twice, once along each of the edges meeting at the corner. If both boundary conditions at this point are Dirichlet and the prescribed solution values agree, then this value is used in an equation of the form (2). If the prescribed solution is discontinuous at the corner, then the average of the two values is used. If one boundary condition is Dirichlet and the other is mixed, then the value prescribed by the Dirichlet condition is used in an equation of the form given above. Finally, if both conditions are mixed or Neumann, then two ' fictitious' points are eliminated using two equations of the form (3). It is possible that equations for which the solution is known at all points on the boundary, have coefficients which are not defined on the boundary. Since this routine calls the user- supplied subroutine PDEF at all points in the domain, including boundary points, arithmetic errors may occur in the user's routine PDEF which this routine cannot trap. If the user has an equation with Dirichlet boundary conditions (i.e., B = 0 at all points on the boundary), but with PDE coefficients which are singular on the boundary, then D03EDF could be called directly only using interior grid points with the user's own discretization. After the equations have been set up as described above, they are checked for diagonal dominance. That is to say, 4 -- k |A | > > |A |, i=1,2,...,n ; j=1,2,...,n . ij -- ij x y k/=4 If this condition is not satisfied then the routine returns with IFAIL = 6. The multigrid routine D03EDF may still converge in this case, but if the coefficients of the first derivatives in the partial differential equation are large compared with the coefficients of the second derivative, the user should consider using upwind differences (SCHEME = 'U'). Since this routine is designed primarily for use with D03EDF, this document should be read in conjunction with the document for that routine. 4. References [1] Wesseling P (1982) MGD1 - A Robust and Efficient Multigrid Method. Multigrid Methods. Lecture Notes in Mathematics. 960 Springer-Verlag. 614--630. 5. Parameters 1: XMIN -- DOUBLE PRECISION Input 2: XMAX -- DOUBLE PRECISION Input On entry: the lower and upper x co-ordinates of the rectangular region respectively, x and x . Constraint: XMIN A B < XMAX. 3: YMIN -- DOUBLE PRECISION Input 4: YMAX -- DOUBLE PRECISION Input On entry: the lower and upper y co-ordinates of the rectangular region respectively, y and y . Constraint: YMIN A B < YMAX. 5: PDEF -- SUBROUTINE, supplied by the user. External Procedure PDEF must evaluate the functions (alpha)(x,y), (beta)(x,y), (gamma)(x,y), (delta)(x,y), (epsilon)(x,y), (phi)(x,y) and (psi)(x,y) which define the equation at a general point (x,y). Its specification is: SUBROUTINE PDEF (X, Y, ALPHA, BETA, GAMMA, 1 DELTA, EPSLON, PHI, PSI) DOUBLE PRECISION X, Y, ALPHA, BETA, GAMMA, DELTA, 1 EPSLON, PHI, PSI 1: X -- DOUBLE PRECISION Input 2: Y -- DOUBLE PRECISION Input On entry: the x and y co-ordinates of the point at which the coefficients of the partial differential equation are to be evaluated. 8 3: ALPHA -- DOUBLE PRECISION Output 4: BETA -- DOUBLE PRECISION Output 5: GAMMA -- DOUBLE PRECISION Output 6: DELTA -- DOUBLE PRECISION Output 7: EPSLON -- DOUBLE PRECISION Output 8: PHI -- DOUBLE PRECISION Output 9: PSI -- DOUBLE PRECISION Output On exit: ALPHA, BETA, GAMMA, DELTA, EPSLON, PHI and PSI must be set to the values of (alpha)(x,y), (beta)(x,y), (gamma)(x,y), (delta)(x,y), (epsilon)(x,y), (phi)(x,y) and (psi)(x,y) respectively at the point specified by X and Y. PDEF must be declared as EXTERNAL in the (sub)program from which D03EEF is called. Parameters denoted as Input must not be changed by this procedure. 6: BNDY -- SUBROUTINE, supplied by the user. External Procedure BNDY must evaluate the functions a(x,y), b(x,y), and c(x,y) involved in the boundary conditions. Its specification is: SUBROUTINE BNDY (X, Y, A, B, C, IBND) INTEGER IBND DOUBLE PRECISION X, Y, A, B, C 1: X -- DOUBLE PRECISION Input 2: Y -- DOUBLE PRECISION Input On entry: the x and y co-ordinates of the point at which the boundary conditions are to be evaluated. 3: A -- DOUBLE PRECISION Output 4: B -- DOUBLE PRECISION Output 5: C -- DOUBLE PRECISION Output On exit: A, B and C must be set to the values of the functions appearing in the boundary conditions. 6: IBND -- INTEGER Input On entry: specifies on which boundary the point (X,Y) lies. IBND = 0, 1, 2 or 3 according as the point lies on the bottom, right, top or left boundary. BNDY must be declared as EXTERNAL in the (sub)program from which D03EEF is called. Parameters denoted as Input must not be changed by this procedure. 7: NGX -- INTEGER Input 8: NGY -- INTEGER Input On entry: the number of interior grid points in the x- and y -directions respectively, n and n . If the seven-diagonal x y equations are to be solved by D03EDF, then NGX-1 and NGY-1 should preferably be divisible by as high a power of 2 as possible. Constraint: NGX >= 3, NGY >= 3. 9: LDA -- INTEGER Input On entry: the first dimension of the array A as declared in the (sub)program from which D03EEF is called. Constraint: if only the seven-diagonal equations are required, then LDA >= NGX*NGY. If a call to this routine is to be followed by a call to D03EDF to solve the seven- diagonal linear equations, LDA >= (4*(NGX+1)*(NGY+1))/3. Note: this routine only checks the former condition. D03EDF, if called, will check the latter condition. 10: A(LDA,7) -- DOUBLE PRECISION array Output On exit: A(i,j), for i=1,2,...,NGX*NGY; j = 1,2,...,7, contains the seven-diagonal linear equations produced by the discretization described above. If LDA > NGX*NGY, the remaining elements are not referenced by the routine, but if LDA >= (4*(NGX+1)*(NGY+1))/3 then the array A can be passed directly to D03EDF, where these elements are used as workspace. 11: RHS(LDA) -- DOUBLE PRECISION array Output On exit: the first NGX*NGY elements contain the right-hand sides of the seven-diagonal linear equations produced by the discretization described above. If LDA > NGX*NGY, the remaining elements are not referenced by the routine, but if LDA >= (4*(NGY+1)*(NGY+1))/3 then the array RHS can be passed directly to D03EDF, where these elements are used as workspace. 12: SCHEME -- CHARACTER*1 Input On entry: the type of approximation to be used for the first derivatives which occur in the partial differential equation. If SCHEME = 'C', then central differences are used. If SCHEME = 'U', then upwind differences are used. Constraint: SCHEME = 'C' or 'U'. Note: generally speaking, if at least one of the coefficients multiplying the first derivatives (DELTA or EPSLON as returned by PDEF) are large compared with the coefficients multiplying the second derivatives, then upwind differences may be more appropriate. Upwind differences are less accurate than central differences, but may result in more rapid convergence for strongly convective equations. The easiest test is to try both schemes. 13: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. Users who are unfamiliar with this parameter should refer to the Essential Introduction for details. On exit: IFAIL = 0 unless the routine detects an error or gives a warning (see Section 6). For this routine, because the values of output parameters may be useful even if IFAIL /=0 on exit, users are recommended to set IFAIL to -1 before entry. It is then essential to test the value of IFAIL on exit. 6. Error Indicators and Warnings Errors or warnings specified by the routine: If on entry IFAIL = 0 or -1, explanatory error messages are output on the current error message unit (as defined by X04AAF). IFAIL= 1 On entry XMIN >= XMAX, or YMIN >= YMAX, or NGX < 3, or NGY < 3, or LDA < NGX*NGY, or SCHEME is not one of 'C' or 'U'. IFAIL= 2 At some point on the boundary there is a derivative in the boundary conditions (B /= 0 on return from a BNDY) and there 2 dd U is a non-zero coefficient of the mixed derivative ------ ddxddy (BETA /= 0 on return from PDEF). IFAIL= 3 A null boundary has been specified, i.e., at some point both A and B are zero on return from a call to BNDY. IFAIL= 4 2 The equation is not elliptic, i.e., 4*ALPHA*GAMMA= 5. 4: LBDCND -- INTEGER Input On entry: indicates the type of boundary conditions at x = XS and x = XF. LBDCND = 0 if the solution is periodic in x, i.e., u(XS,y,z)=u(XF,y,z). LBDCND = 1 if the solution is specified at x = XS and x = XF. LBDCND = 2 if the solution is specified at x = XS and the derivative of the solution with respect to x is specified at x = XF. LBDCND = 3 if the derivative of the solution with respect to x is specified at x = XS and x = XF. LBDCND = 4 if the derivative of the solution with respect to x is specified at x = XS and the solution is specified at x = XF. Constraint: 0 <= LBDCND <= 4. 5: BDXS(MDIMF,N+1) -- DOUBLE PRECISION array Input On entry: the values of the derivative of the solution with respect to x at x = XS. When LBDCND = 3 or 4, BDXS (j,k)=u (XS,y ,z ), for j=1,2,...,M+1; k=1,2,...,N+1. x j k When LBDCND has any other value, BDXS is not referenced. 6: BDXF(MDIMF,N+1) -- DOUBLE PRECISION array Input On entry: the values of the derivative of the solution with respect to x at x = XF. When LBDCND = 2 or 3, BDXF(j,k) = u (XF,y ,z ), for j=1,2,...,M+1; k=1,2,...,N+1. x i k When LBDCND has any other value, BDXF is not referenced. 7: YS -- DOUBLE PRECISION Input On entry: the lower bound of the range of y, i.e., YS <= y <= YF. Constraint: YS < YF. 8: YF -- DOUBLE PRECISION Input On entry: the upper bound of the range of y, i.e., YS <= y <= YF. Constraint: YS < YF. 9: M -- INTEGER Input On entry: the number of panels into which the interval (YS,YF) is subdivided. Hence, there will be M+1 grid points in the y-direction given by y = YS + (j-1)*(delta)y for j j=1,2,...,M+1, where (delta)y=(YF-YS)/M is the panel width. Constraint: M >= 5. 10: MBDCND -- INTEGER Input On entry: indicates the type of boundary conditions at y = YS and y = YF. MBDCND = 0 if the solution is periodic in y, i.e., u(x,YF,z)=u(x,YS,z). MBDCND = 1 if the solution is specified at y = YS and y = YF. MBDCND = 2 if the solution is specified at y = YS and the derivative of the solution with respect to y is specified at y = YF. MBDCND = 3 if the derivative of the solution with respect to y is specified at y = YS and y = YF. MBDCND = 4 if the derivative of the solution with respect to y is specified at y = YS and the solution is specified at y = YF. Constraint: 0 <= MBDCND <= 4. 11: BDYS(LDIMF,N+1) -- DOUBLE PRECISION array Input On entry: the values of the derivative of the solution with respect to y at y = YS. When MBDCND = 3 or 4, BDYS (i,k)=u (x ,y ,z ), for i=1,2,...,L+1; k=1,2,...,N+1. y i s k When MBDCND has any other value, BDYS is not referenced. 12: BDYF(LDIMF,N+1) -- DOUBLE PRECISION array Input On entry: the values of the derivative of the solution with respect to y at y = YF. When MBDCND = 2 or 3, BDYF (i,k)=u (x ,YF,z ), for i=1,2,...,L+1; k=1,2,...,N+1. y i k When MBDCND has any other value, BDYF is not referenced. 13: ZS -- DOUBLE PRECISION Input On entry: the lower bound of the range of z, i.e., ZS <= z <= ZF. Constraint: ZS < ZF. 14: ZF -- DOUBLE PRECISION Input On entry: the upper bound of the range of z, i.e., ZS <= z <= ZF. Constraint: ZS < ZF. 15: N -- INTEGER Input On entry: the number of panels into which the interval (ZS,ZF) is subdivided. Hence, there will be N+1 grid points in the z-direction given by z =ZS+(k-1)*(delta)z, for k k=1,2,...,N+1, where (delta)z=(ZF-ZS)/N is the panel width. Constraint: N >= 5. 16: NBDCND -- INTEGER Input On entry: specifies the type of boundary conditions at z = ZS and z = ZF. NBDCND = 0 if the solution is periodic in z, i.e., u(x,y,ZF)=u(x,y,ZS). NBDCND = 1 if the solution is specified at z = ZS and z = ZF. NBDCND = 2 if the solution is specified at z = ZS and the derivative of the solution with respect to z is specified at z = ZF. NBDCND = 3 if the derivative of the solution with respect to z is specified at z = ZS and z = ZF. NBDCND = 4 if the derivative of the solution with respect to z is specified at z = ZS and the solution is specified at z = ZF. Constraint: 0 <= NBDCND <= 4. 17: BDZS(LDIMF,M+1) -- DOUBLE PRECISION array Input On entry: the values of the derivative of the solution with respect to z at z = ZS. When NBDCND = 3 or 4, BDZS (i,j)=u (x ,y ,ZS)=u(x,y,z), for i=1,2,...,L+1; z i j j=1,2,...,M+1. When NBDCND has any other value, BDZS is not referenced. 18: BDZF(LDIMF,M+1) -- DOUBLE PRECISION array Input On entry: the values of the derivative of the solution with respect to z at z = ZF. When NBDCND = 2 or 3, BDZF (i,j)=u (x ,y ,ZF)=u(x,y,z), for i=1,2,...,L+1; z i j j=1,2,...,M+1. When NBDCND has any other value, BDZF is not referenced. 19: LAMBDA -- DOUBLE PRECISION Input On entry: the constant (lambda) in the Helmholtz equation. For certain positive values of (lambda) a solution to the differential equation may not exist, and close to these values the solution of the discretized problem will be extremely ill-conditioned. If (lambda)>0, then D03FAF will set IFAIL to 3, but will still attempt to find a solution. However, since in general the values of (lambda) for which no solution exists cannot be predicted a priori, the user is advised to treat any results computed with (lambda)>0 with great caution. 20: LDIMF -- INTEGER Input On entry: the first dimension of the arrays F, BDYS, BDYF, BDZS and BDZF as declared in the (sub)program from which D03FAF is called. Constraint: LDIMF >= L + 1. 21: MDIMF -- INTEGER Input On entry: the second dimension of the array F and the first dimension of the arrays BDXS and BDXF as declared in the (sub)program from which D03FAF is called. Constraint: MDIMF >= M + 1. 22: F(LDIMF,MDIMF,N+1) -- DOUBLE PRECISION array Input/Output On entry: the values of the right-side of the Helmholtz equation and boundary values (if any). F(i,j,k)=f(x ,y ,z ) i=2,3,...,L, j=2,3,...,M and k i j k =2,3,...,N. On the boundaries F is defined by LBDCND F(1,j,k) F(L+1,j,k) 0 f(XS,y ,z )f(XS,y ,z ) j k j k 1 u(XS,y ,z )u(XF,y ,z ) j k j k 2 u(XS,y ,z )f(XF,y ,z ) j=1,2,...,M+1 j k j k 3 f(XS,y ,z )f(XF,y ,z ) k=1,2,...,N+1 j k j k 4 f(XS,y ,z )u(XF,y ,z ) j k j k MBDCND F(i,1,k) F( i,M+1,k) 0 f(x ,YS,z )f(x ,YS,z ) i k i k 1 u(x ,YS,z )u(x ,YF,z ) i k i k 2 u(x ,YS,z )f(x ,YF,z ) i=1,2,...,L+1 i k i k 3 f(x ,YS,z )f(x ,YF,z ) k=1,2,...,N+1 i k i k 4 f(x ,YS,z )u(x ,YF,z ) i k i k NBDCND F(i,j,1) F(i,j,N+1 ) 0 f(x ,y ,ZS)f(x ,y ,ZS) i j i j 1 u(x ,y ,ZS)u(x ,y ,ZF) i j i j 2 u(x ,y ,ZS)f(x ,y ,ZF) i=1,2,...,L+1 i j i j 3 f(x ,y ,ZS)f(x ,y ,ZF) j=1,2,...,M+1 i j i j 4 f(x ,y ,ZS)u(x ,y ,ZF) i j i j Note: if the table calls for both the solution u and the right-hand side f on a boundary, then the solution must be specified. On exit: F contains the solution u(i,j,k) of the finite difference approximation for the grid point ( x ,y ,z ) for i=1,2,...,L+1, j=1,2,...,M+1 and k=1,2,...,N+1. i j k 23: PERTRB -- DOUBLE PRECISION Output On exit: PERTRB = 0, unless a solution to Poisson's equation ((lambda)=0) is required with a combination of periodic or derivative boundary conditions (LBDCND, MBDCND and NBDCND = 0 or 3). In this case a solution may not exist. PERTRB is a constant, calculated and subtracted from the array F, which ensures that a solution exists. D03FAF then computes this solution, which is a least-squares solution to the original approximation. This solution is not unique and is unnormalised. The value of PERTRB should be small compared to the right-hand side F, otherwise a solution has been obtained to an essentially different problem. This comparison should always be made to insure that a meaningful solution has been obtained. 24: W(LWRK) -- DOUBLE PRECISION array Workspace 25: LWRK -- INTEGER Input On entry: the dimension of the array W as declared in the (sub)program from which D03FAF is called. LRWK>=2*(N+1)*max(L,M)+3*L+3*M+4*N+6 is an upper bound on the required size of W. If LWRK is too small, the routine exits with IFAIL = 2, and if on entry IFAIL = 0 or IFAIL = - 1, a message is output giving the exact value of LWRK required to solve the current problem. 26: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. Users who are unfamiliar with this parameter should refer to the Essential Introduction for details. On exit: IFAIL = 0 unless the routine detects an error or gives a warning (see Section 6). For this routine, because the values of output parameters may be useful even if IFAIL /=0 on exit, users are recommended to set IFAIL to -1 before entry. It is then essential to test the value of IFAIL on exit. 6. Error Indicators and Warnings Errors or warnings specified by the routine: If on entry IFAIL = 0 or -1, explanatory error messages are output on the current error message unit (as defined by X04AAF). IFAIL= 1 On entry XS >= XF, or L < 5, or LBDCND < 0, or LBDCND > 4, or YS >= YF, or M < 5, or MBDCND < 0, or MBDCND > 4, or ZS >= ZF, or N < 5, or NBDCND < 0, or NBDCND > 4, or LDIMF < L + 1 > 0, or MDIMF < M + 1. IFAIL= 2 On entry LWRK is too small. IFAIL= 3 On entry (lambda) > 0. 7. Accuracy None. 8. Further Comments The execution time is roughly proportional to L*M*N*(log L+log M+5), but also depends on input parameters 2 2 LBDCND and MBDCND. 9. Example The example solves the Helmholz equation 2 2 2 dd u dd u dd u ----+ ----+ ----+(lambda)u=f(x,y,z) 2 2 2 ddx ddy ddz (pi) for (x,y,z)is in [0,1]*[0,2(pi)]*[0, ----] where (lambda)=-2, and 2 f(x,y,z) is derived from the exact solution 4 u(x,y,z)=x sin(y)cos(z). The equation is subject to the following boundary conditions, again derived from the exact solution given above. u(0,y,z) and u(1,y,z) are prescribed (i.e., LBDCND = 1). u(x,0,z)=u(x,2(pi),z) (i.e., MBDCND = 0). (pi) u(x,y,0) and u (x,y, ----) are prescribed (i.e. NBDCND = 2). x 2 The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page}