An improved algorithm and a Fortran 90 module for computing the conical function [formula omitted

Abstract

In this paper we describe an algorithm and a Fortran 90 module ( Conical) for the computation of the conical function P − 1 2 + i τ m ( x ) for x > − 1 , m ⩾ 0 , τ > 0 . These functions appear in the solution of Dirichlet problems for domains bounded by cones; because of this, they are involved in a large number of applications in engineering and physics. In the Fortran 90 module, the admissible parameter ranges for computing the conical functions in standard IEEE double precision arithmetic are restricted to ( x , m , τ ) ∈ ( − 1 , 1 ) × [ 0 , 40 ] × [ 0 , 100 ] and ( x , m , τ ) ∈ ( 1 , 100 ) × [ 0 , 100 ] × [ 0 , 100 ] . Based on tests of the three-term recurrence relation satisfied by these functions and direct comparison with Maple, we claim a relative accuracy close to 10 − 12 in the full parameter range, although a mild loss of accuracy can be found at some points of the oscillatory region of the conical functions. The relative accuracy increases to 10 − 13 – 10 − 14 in the region of the monotonic regime of the functions where integral representations are computed ( − 1 < x < 0 ). Program summary Program title: Conical Catalogue identifier: AELD_v1_0 Program summary URL: http://cpc.cs.qub.ac.uk/summaries/AELD_v1_0.html Program obtainable from: CPC Program Library, Queenʼs University, Belfast, N. Ireland Licensing provisions: Standard CPC licence, http://cpc.cs.qub.ac.uk/licence/licence.html No. of lines in distributed program, including test data, etc.: 5387 No. of bytes in distributed program, including test data, etc.: 49 615 Distribution format: tar.gz Programming language: Fortran 90 Computer: Any supporting a FORTRAN compiler Operating system: Any supporting a FORTRAN compiler RAM: A few MB Classification: 4.7 Nature of problem: Conical functions appear in a large number of applications because these functions are the natural function basis for solving Dirichlet problems bounded by conical domains. Also, they are the Kernel of the Mehler–Fock transform. Solution method: The algorithm uses different methods of computation depending on the range of parameters: asymptotic expansions, quadrature methods and recurrence relations. Restrictions: In order to avoid underflow/overflow problems, the admissible parameter ranges for computing the conical functions in standard IEEE double precision arithmetic are restricted to ( x , m , τ ) ∈ ( − 1 , 1 ) × [ 0 , 40 ] × [ 0 , 100 ] and ( x , m , τ ) ∈ ( 1 , 100 ) × [ 0 , 100 ] × [ 0 , 100 ] . Additional comments: The module Conical uses a Fortran 90 version of the routine dkia (developed by the authors) for computing the modified Bessel functions K i a ( x ) and its derivative. This routine is included in the distribution file and is also available at http://toms.calgo.org. Running time: Depending on the parameter range: when numerical quadrature is used (for x < 0 ), the algorithm is 10–20 times slower than the computations made using asymptotic expansions + recurrence relations.