Donald E. Amos

Algorithm 644: A portable package for Bessel functions of a complex argument and nonnegative order

This algorithm is a package of subroutines for Computing Bessel functions <italic>H<subscrpt>v</subscrpt></italic><supscrpt>(1)</supscrpt>(<italic>z</italic>), <italic>H<subscrpt>v</subscrpt></italic><supscrpt>(2)</supscrpt>(<italic>z</italic>), <italic>I<subscrpt>v</subscrpt>(z)</italic>, <italic>J<subscrpt>v</subscrpt>(z)</italic>, <italic>K<subscrpt>v</subscrpt>(z)</italic>, <italic>Y<subscrpt>v</subscrpt>(z)</italic> and Airy functions Ai(<italic>z</italic>), Ai′(<italic>z</italic>), Bi(<italic>z</italic>), Bi′(<italic>z</italic>) for orders<italic>v</italic>≥0 and complex <italic>z</italic> in −π<arg <italic>z</italic>≤π. Eight callable subroutines and their double-precision counterparts are provided. Exponential scaling and sequence generation are auxiliary options.

CDC 6600 Subroutines IBESS and JBESS for Bessel Functions I υ (x) and J υ (x), x≥0,υ≥0

Subroutines IBESS and JBESS for L(x) and J,(x), x >_ O, p > O, are presented which implement , in similar fashions, the power series, the asymptotic expansion for x ~ ~, and the uniform asymptotic expansion for ~-* ¢~ as basic elements. Where the asymptotic expansion for x-~ o0 is used, sequences are started by two evaluations of the expansion followed by backward recursion in IBESS and forward recursion in JBESS for p < x. Two evaluations of the series also start a backward reeursive sequence where x _< ½. Except for very large z and v, where it is more economical to perform two evaluations of a basic formula, other sequences are started by the Miller algorithm normalized by one evaluation of either the power series or the expansion for ~-~ ~o. Recursion is used for pairs (~, x) not covered by one of the basic formulas. A scaling option to remove the exponential growth of I,(x) is also provided together with appropriate diagnostics for underflow or overflow in both IBESS and JBESS.

Solution of an initial-boundary value problem for heat conduction in a parallelepiped by time partitioning

Abstract An initial-boundary value problem for transient heat conduction in a rectangular parallelepiped is studied. Solutions for the temperature and heat flux are represented as integrals involving the Greens function (GF), the initial and boundary data, and volumetric energy generation. Use of the usual GF obtained by separation of variables leads to slowly convergent series. To circumvent this difficulty, the dummy time interval of integration is partitioned into a short time and a long time subintervals where the GFs are approximated by their small and large time representations. This paper deals with the analysis and implementation of this time partitioning method.

A remark on Algorithm 644: “A portable package for Bessel functions of a complex argument and nonnegative order”

Algorithm 644 computes all major Bessel functions of a complex argument and of nonnegative order. Single-precision routine CBESY and double-precision routine ZBESY are modified to reduce the computation time for the Y Bessel function by approximately 25% over a wide range of arguments and orders. Quick check (driver) programs that exercise the package are also modified to make tests more meaningful over larger regions of the complex plane.

Algorithm 556: Exponential Integrals [S13]

DESCRIPTION The Fortran subroutine EXPINT given here is an implementation of [1]. EX-PINT has four machine-dependent parameters XCUT, XLIM, ETOL, and EU-LER which are set into DATA statements. XCUT is a breakpoint such that for x _ XCUT, the series is evaluated and for x > XCUT, the Miller algorithm is applied. XLIM is the approximate underflow limit for e-=, x _ 0, and ETOL is nominally set to 1.E-D where D is the number of base 10 digits in a word. would be appropriate for IBM single-precision arithmetic. The two choices for XCUT reflect the fact that there is a loss of up to two digits on 1 < x _< 2 with the series evaluation. This loss can be tolerated on longer word length machines, but not on shorter word length machines. Maximum accuracy can always be achieved with XCUT = 1. However, for longer word lengths where D = 14, the reduction in computation by moving XCUT from 1 to 2, achieved at the expense of two digits of accuracy, seems to be a worthwhile trade-off, and D ffi 12 reflects this modification for CDC machines. D = 12 is also more consistent with the accuracy attainable from E X P (-X) near the underflow limit X = 667. While the subroutine EXPINT is almost portable, the function DIGAM, which computes the psi function at integer arguments, is supplied as a CDC 6600-7600 Fortran function. Modifications for other machines can be made easily from eq. Algorithms • 421 T h e convergence of eq. (2) in [1] is so rapid that the m = n-1 term, which requires ~(n), is reached for only small values of n. A table of 100 values suffices for virtually all single-precision implementations of E X P I N T. An initialization step to generate a higher precision table might be appropriate if multiple precision is anticipated. For CDC single precision or IBM double precision, 36 values suffice for n as high as 10 le and relative errors of 10-14. was used to check out the routine on parameter ranges were also made at tolerances T O L as low as 10-9. The quadrature~could be relied upon to give the requested accuracy down to T O L-10-1~ over the full exponent range of the CDC 6600. For large x, argument reduction in computing e-~ will result in decreased accuracy, up to three …

Algorithms 683: a portable FORTRAN subroutine for exponential integrals of a complex argument

The algorithm computes exponential integrals <italic>E<subscrpt>n</subscrpt></italic>(<italic>z</italic>) for integer orders <italic>n</italic> ≥ 1 and complex <italic>z</italic> in -π < arg <italic>z</italic> ≤ π. Both single and double precision subroutines are provided. Exponential scaling and stable sequence generation are auxiliary options.

Remark on algorithm 644

Algorithm 644 computes all major Bessel functions of a complex argument and of nonnegative order. Single-precision routine CBESY and double-precision routine ZBESY are modified to reduce the comput...

Computation of exponential integrals of a complex argument

Previous work on exponential integrals of a real argument has produced a computational algorithm which implements backward recurrence on a three-term recurrence relation (Miller algorithm). The process on which the algorithm is based involves the truncation of an infinite sequence with a corresponding analysis of the truncation error. This error is estimated by means of (asymptotic) formulas, which are not only applicable to the real line, but also extendible to the complex plane. This fact makes the algorithm extendible to the complex plane also. However, the rate of convergence decreases sharply when the complex argument is close to the negative real axis. In order to make the algorithm more efficient, analytic continuation into a strip about the negative real axis is carried out by limited use of power series. The analytic details needed for a portable computational algorithm are presented.

Diffusion Penetration Time for Transient Heat Conduction

The time duration for processes involving transient thermal diffusion can be a critical piece of information related to thermal processes in engineering applications. Analytical solutions must be used to calculate these types of time durations because the boundary conditions in such cases can be effectively like semi-infinite conditions. This research involves an investigation into analytical solutions for six geometries, including one-dimensional cases for Cartesian, cylindrical, and spherical coordinates. The fifth case involves a heated surface on the inside of a hole bored through an infinite body, which is a one-dimensional problem in cylindrical coordinates. The sixth case involves two-dimensional conduction from a point heat source on the surface of a slab subjected to insulated boundary conditions elsewhere. The mathematical modeling for this case is done in cylindrical coordinates. For each geometric configuration, a relationship is developed to determine the time required for a temperature rise ...

Remark on “Algorithm 556: Exponential Integrals”

The following changes in the FORTRAN program rescale the normalization of the Miller algorithm so that premature underflow without a diagnostic does not occur. Reference [1, Sec. 6], shows the calculation of XLIM, a machine-dependent parameter related to the exponent range of the machine. Single-and double-precision portable versions of Algorithm 556 (with this modification included) are also provided with the code described in [1].