Legendre numerical Kent State

Electronic Transactions on Numerical Analysis. Volume 9, 1999, pp. 137-146.Copyright 1999, Kent State University. ISSN 1068-9613.ETNA Kent State University [email protected] EVALUATION OFASSOCIATED LEGENDRE FUNCTIONSOFFTHECUTAND PARABOLIC CYLINDER FUNCTIONS JAVIER SEGURAyzANDAMPARO GILyx Abstract. We review a set of algorithms to evaluate associated Legendre functions off the cut; in particular, we consider prolate spheroidal, oblate spheroidal andtoroidal harmonics. Asimilar schemecanbeapplied tootherfam-ilies of special functions like Bessel and parabolic cylinder functions; we will describe the corresponding algorithmfor the evaluation of parabolic cylinder functions. Key words. computation of special functions, Legendre functions, parabolic cylinder functions. AMSsubject classifications. 65D20, 33-04, 33C05, 33A70. 1. Introduction. The evaluation of associated Legendre (ALF) and parabolic cylinder functions (PCF) is a matter of relevance because these functions appear in the solution of Dirichlet problems in different geometries [12]. Then, they show up in a vast number ofapplications [12, 8, 9] in different fields such as, for instance, lattice field theory[5], ther- monuclearfusion[16],biology[10]orcristallography[18]. Recently,a seriesof codesto evaluateALF[8, 9, 24]and PCF[23] havebeendeveloped, fillinga considerablegapin numericallibraries. For associated Legendre functions off the cut there was no available routine; only Gautschi [6], in 1965, presented a set of algorithms in ALGOL60 to evaluate them. Our approach is similar to Gautschi’s: Legendre functions off the cut satisfy three term recur- rencerelations, beingone ofthe independentsolutionsa minimalsolution[28, 7]. However, ourcodehassomeimportantdifferencesfromGautschi’swhich,in fact,allowsittobemoreaccurateandvalidfora largerrangeoftheparameters[24]. Other examples of families of real functions of real variable satisfying three term re- currences with a minimal solution are Bessel and Modified Bessel functions and parabolic cylinderfunctions. Bessel functionshavebeenbroadlydiscussedintheliteratureandmanyalgorithmswith differentcharacteristicsexist[19,2,27,26,22]. Buttherewasaconsiderablelackofnumer- icalalgorithmsforPCF: therewasonlyonepublishedprogramtoevaluatePCFs[25],whichas we discussed [23], has serious problems. Different approachesto the evaluation of PCFs canbefoundin [11, 20,15,21]. ALFs,PCFsandBesselfunctionsareclassicalinthesensethatallstandardbooksonspe- cialfunctions[1,12,28]devoteatleastachaptertothem. Inaddition,thetaskofdeveloping numericalmethodstoevaluatetheclassicalspecialfunctionshasgainedrenewedinterestdue totheongoingprogramtorevisetheAbramowitz &StegunHandbookonMathematicalfunc- tions [14]. A comprehensivenumericallibrary to generate valuesfor all functionsdescribed insuchrevisedversionisintendedto bebuilt. 2. Legendreandparaboliccylinderfunctions: definitionand properties. 2.1. Associated Legendre functions off the cut. The associated Legendre functions P m (z)andQm (z)[1] aresolutionsofthedifferentialequation Received November 1, 1998. Accepted for publicaton December 1, 1999. Recommended by F. Marcell´ an. yInstituto de Bioingenier´ ıa, Universidad Miguel Hern´ andez, Edificio La Galia, 03202-Elche (Alicante), Spain. z([email protected]) x([email protected]) 137 ETNA Kent State University [email protected] 138 Legendre functions and parabolic cylinder functions (1−z2)u00−2zu0+ (+1 )−m2 1−z2 u=0; (2.1) where,in mostpracticalsituations, misa nonnegativeinteger. FromnowonwewillconsiderassociatedLegendrefunctionswith zoutsidetheinterval [−1;1], that is, associated Legendre functions off the cut. For the evaluation of associated Legendrefunctionsonthecut,see[17]. TherecurrencerelationssatisfiedbytheassociatedLegendrefunctionsoffthecut(ALF), bothoverdegrees ’s andorders m’s are (−m+1 )Pm +1(z)−(2+1 )zPm (z)+(+m)Pm −1(z)=0; (2.2) Pm+1 (z)+2mz (z2−1)1=2Pm (z)−(−m+1 ) (+m)Pm−1 =0; (2.3) wherethesamerelationsapplyforthe Q’s. TheWronskianrelationbetween P’s and Q’s is W(Pm (z);Qm (z)) =Γ(+m+1 ) Γ(−m+1 )(−1)m 1−z2: (2.4) From which follow two useful relations between consecutive degrees (eq.(5)) and orders (eq.(6)) Pm (z)Qm −1(z)−Pm −1(z)Qm (z)=Γ(+m) Γ(−m+1 )(−1)m; (2.5) Pm (z)Qm+1 (z)−Pm+1 (z)Qm (z)=Γ(+m+1 ) Γ(−m+1 )(−1)m p z2−1: (2.6) For half-integer degrees n−1=2;n=0;1;2; :::and real arguments zx>1, the functions fPm n−1=2(x);Qm n−1=2(x)gare called toroidal harmonics .W h e n is an integer n=0;1;2; :::and for real arguments x>1, the functions fPm n(x);Qm n(x)gare calledprolate spheroidal harmonics , while, for purely imaginary arguments, the functions fPm n(ix);Qm n(ix)gwithx>0areknownas oblatespheroidalharmonics . Both prolate spheroidaland toroidalharmonicsare real functionsof the real variable x. Theoblatespheroidalharmonics fPm n(ix);Qm n(ix)gcanberealorimaginaryvaluedforreal x;however,thenewsetoffunctions fRm n(x);Tm n(x)gx>0;n0definedby Rm n(x)=exp(−in 2)Pm n(ix); Tm n(x)=iexp(in 2)Qm n(ix)(2.7) arerealfunctionsoftherealvariable x,andmoreconvenientfornumericalevaluation. From nowon,wewillreferto Rm n(x)andTm n(x)asoblatespheroidalharmonics (OSH)ofthefirst andsecondkindsrespectively. Reference[9]isthefirstonetoprovideanumericalalgorithm tocomputeOSHs. ETNA Kent State University [email protected] J.Segura and A.Gil 139 2.2. Parabolic cylinder functions. The parabolic cylinder functions V(a; x)and U(a; x)[1] aresolutionsofthe differentialequation y00−(a+1 4x2)y=0: (2.8) TheVsandUssatisfy thefollowingrecurrencerelations: V(a+1;x)=xV(a; x)+(a−1=2)V(a−1;x); (2.9) U(a−1;x)=xU(a; x)+(a+1=2)U(a+1;x): (2.10) TheWronskianrelationbetween V’s and U’s is WfU(a; x);V(a; x)g=p 2= (2.11) fromwhichit followsthat (a−1=2)U(a; x)V(a−1;x)+U(a−1;x)V(a; x)=r 2 : (2.12) 3. Recurrencerelationsandstability. BothassociatedLegendre fP;Qgandparabolic cylinder functions fU;Vghave two common characteristics: they satisfy three-term recur- rencerelationsandoneofthesolutionsisminimal. A threetermrecurrencerelation yk+1+akyk+bkyk−1=0 (3.1) issaid toadmitaminimalsolutionwhenthereexisttwolinearlyindependentsolutions y# k,y" ksuchthat lim k!1y# k y" k=0 ; (3.2) the solution y# kis called minimal solution(whichis unique)while y" kis adominantsolution. The recurrence relation should be applied backwards to evaluate the minimal solution, and neverforward,sinceanysmallroundingerrorwouldintroduceadominantcomponent. Onthe otherhand,therecurrencerelationhastobeappliedforwardtocalculatedominantsolutions. ImportantresultsareprovidedbyPerron’s[28,19]andPincherle’s[3,19]theorems: Per- ron’s theorem helps in studying the stability of recurrences and the existence of a minimal solution. On the other hand, Pincherle’s theorem guarantees the existence of a continued fraction(CF)fortheratioofconsecutiveminimalsolutions y# k=y# k−1andgivesaprescription to estimate the speed of convergenceof the resulting CF; in case the recurrence(3.1)admits minimalsolutionPincherle’stheoremstatesthat theratiocanbeevaluatedintheform: ETNA Kent State University [email protected] 140 Legendre functions and parabolic cylinder functions y# k=y# k−1=−bk ak−bk+1 ak+1−::: (3.3) It is easy to check that Qm is the minimal solution of the three term recurrencerelation (2.2) while Pm is a dominant solution. Correspondingly, for oblate spheroidal harmonics, theminimalsolutionis Tm nandthedominantoneis Rm n. Thecharacterofthefunctions Qm , Pm changes if one consider the recurrence relation given by eq.(2.3): Pm is the minimal solution while Qm is a dominant one. Notice that, in case is integer, because Pm n=0 whenm>nandn;mare integers, the CF (3.3) for the ratio Pm n−1=Pm−1 nbecomes a finite continuedfraction. For prolate and oblate spheroidal harmonics we always assume nm; the recurrence overnstartingwith n=mandn=m+1willbeenoughfortheirevaluation. However,the evaluationofTH, needsbothrecurrences(over mandn). In the case of parabolic cylinderfunctions, one can easily establish the character of U’s as the minimal solution of recurrence (2.10) and the character of Vs as a dominant one of (2.9). We have focused our attention on integer and half-integer values of the order aand non-negativearguments xwhicharethe casesofgreatestapplicability. 4. Numerical evaluation of ALF and PCF. The numerical evaluation of PSH, OSH andPCFfollowa similarschemeandthe procedurecanbedescribedintermsofa singleba- sic algorithm. The main differencesare in the evaluation of the starting values to “feed” therecurrences, the study of the convergence of the continued fraction (and substitution when- ever it converges slowly) and the handling of possible numerical overflows. For issues of convergenceof the CFs and controlof overflows, we refer to [8, 9, 23, 24]. We describe the basic algorithm and the evaluation of the starting values. Also, we will explicitly show the resultingalgorithmforOSH. The algorithms for the evaluation of TH are considerably more involved. In this case, both recurrences have to be combined in the algorithm. We will present one of the threealgorithmsdescribedinref. [24] 4.1. Basic algorithm. The main ingredients of our algorithms are the character of the functionsasminimalordominantsolutionsofathreetermrecurrencerelationandtheWron-skian relatingbothsolutions. Essentially the procedurecan be describedas follows: Givena threetermrecurrencerelation y k+1+akyk+bkyk−1=0;k1 (4.1) withy# ktheminimalsolutionand y" ka dominantone,andconsidering y# ky" k−1+ck(x)y# k−1y" k=dk(x) (4.2) the Wronskian relating both solutions, the following steps are considered to evaluate the set fy" k;y# k;k=0;1; :::Kg: /circlecopyrtEvaluate y" 0;y" 1. /circlecopyrtUse forwardrecurrencetoobtaintheset fy" 0;y" 1; :::; y" Kg. /circlecopyrtCombine y# K=y# K−1=−bKaK−bK+1aK+1−:::with the Wronskian relation (4.2) and y" K;y" K−1toget y# K;y# K−1. ETNA Kent State University [email protected] J.Segura and A.Gil 141 /circlecopyrtUse backwardrecurrencetoobtain fy# K;y# K−1; :::; y# 0g. Some interesting features of the algorithm described are: first, unlike Miller’s method, no renormalization has to be carried out. This is important in order to have a good control of accuracyandoverflows. Second,bothdominantandminimalsolutionscan be obtainedat the same time; this feature is of interest when both solutions are needed, as happens whensolvingDirichletproblems. The basic ingredient in Gautschi’s codes for PSH and TH was also the application of recurrence relations. However, although the underlying theory is the same as in Gautschi’s codes,ourapproachleadstoalgorithmsvalidforlargerrangesoftheparameters,muchfaster whenseveralorders/degreesareneeded,andwithhigherprecision[24]. On the other hand, surprisingly, recurrence relations where rarely used for PCF and the associated continuedfractionand, usefulas it is, was notconsideredin [11, 20, 15, 21]. Our code for PCFs, as we discussed in [23], solves the problems of Taubmann’s code [25] and enlargesconsiderablytherangesofparameters. In principle, we only need to evaluate the two starting values for the recurrences. How- ever,onealso needstotake careofpossiblebadconvergenceoftheCFs (takingintoaccount Pincherle’stheorem)andtoreplacetheCFbyseriesorasymptoticexpansionswhenneeded. Formoredetailssee[23, 24]. Letus now summarizehowthe evaluationof the starting valuesis carriedin eachof the casesdescribed. 4.2. Evaluation of the starting values for the recurrences. To feed the recurrence relationsweneedtwostartingvalues,whichareevaluatedasfollows: /circlecopyrtProlateandoblatespheroidalharmonics: For prolate and oblate spheroidal harmonics a closed expression can be found for theinitial values: P m m(x)=( 2 m−1)!!(x2−1)m=2;Pm m+1(x)=x(2m+1 )Pm m(x); (4.3) Rm m(x)=( 2 m−1)!!(x2+1 )m=2;Rm m+1(x)=x(2m+1 )Rm m(x): (4.4) /circlecopyrtToroidalharmonics: Inthiscaseweusetherelationof Q0 −1=2andQ1 −1=2withtheellipticintegrals Eand K: Q0 −1=2(x)=p 2=(x+1 )K(p 2=(x+1 ) ); Q1 −1=2(x)=−E(p 2=(x+1 ) )=p 2(x−1);(4.5) andweevaluate EandKbymeansoftheCarlson’sduplicationtheorem[4]. /circlecopyrtParaboliccylinderfunctionsofintegerorder a0: For paraboliccylinderfunctions V(a; x)with integervaluesofthe parameter a,w e considertherelationof V(0;x)andV(1;x)with themodifiedBesselfunctions I V(0;x)=px 2/parenleftbig I−1=4(x2=4) +I1=4(x2=4)
(4.6) ETNA Kent State University [email protected] 142 Legendre functions and parabolic cylinder functions V(1;x)=x3=2 4/parenleftbig I−1=4(x2=4) +I1=4(x2=4) +I−3=4(x2=4) +I3=4(x2=4)
:(4.7) To evaluate the Bessel functions we have followed the scheme of reference [19], complementedwithanasymptoticexpansion[23]forlarge x. /circlecopyrtParaboliccylinderfunctionsofhalf-integerorder a1=2: In the half-integer case, the expressions for the initial parabolic cylinder functionsaresimpler: V(1=2;x)=r 2 ex2=4;V(3=2;x)=r 2 xex2=4: (4.8) 4.3. An explicit example: oblate spheroidal harmonics. As an example of the basic algorithm,we showthecorrespondingtotheevaluationofoblatespheroidalharmonics: Letrn(x)=Rm m+n(x)andtn(x)=Tm m+n(x). The following steps are followed to evaluatetheset frn;tn;n=0; :::; Ng: /circlecopyrtEvaluate r0(x)=( 2 m−1)!!(x2+1 )m=2>0andr1(x)=x(2m+1 )r0(x)0. /circlecopyrtApplytherecurrencerelation rn+1=1 n+1[(2n+2m+1 )xrn(x)+(n+2m)rn−1(x)]0 (forward)uptoa maximumdegree n=N. /circlecopyrtUse the Wronskianrelation,combinedwith the CF for HN(x)=tN(x)=tN−1(x)(con- vergentfor x>0)toobtain tN−1=(2m+N−1)! N!(−1)m 1 rN(x)+HN(x)rN−1(x); tN(x)=tN−1(x)HN(x): /circlecopyrtUsing tN;tN−1asstartingvalues,therecurrencerelation tn−1(x)=1 (n+2m)[(n+1 )tn+1(x)+( 2 n+2m+1 )xtn(x)] isappliedbackwards. Taking into accountthat Rm n(x)0;Tm n(x)08x0, one can see that no subtrac- tionstakeplaceinapplyingtheforwardandthebackwardrecurrencesfor Rm n(x)andTm n(x) respectively. Then,nosignificantroundofferrorsareexpectedtooccur. 4.4. A more involved example: toroidal harmonics.. The algorithm for toroidal har- monics(TH)deservesaseparateanalysis. Themaindifficulty,comparedwithOSHandPSH, concernsthestarting pointofthe algorithm. ForTHwedonothaveclosedformexpressions like(4.3),(4.4),whichallowedtheevaluationofOSHandPSHforfixed musingonlyrecur- rence(2.2). ForTHoneneedstouse recurrence(2.2)combinedwith(2.3). ETNA Kent State University [email protected] J.Segura and A.Gil 143 Ascommented,the Q’sareminimalandthe P’sdominantforrecursionoverthedegree n while,forrecursionovertheorder m,theP’saretheminimalsolutionandthe Q’sdominant. This“dual” behaviortogetherwith the two associated CF’s and the two Wronskianrelations (2.5) and (2.6) makes it possible to reach any order mor degree nfrom two starting and consecutivevalues. Using this fact, the algorithmfortoroidalharmonicscan be summarized asfollows: Theset fPm n−1=2;Qm n−1=2g,n=0;1; :::; N +1,m=0;1; :::; Mcanbegeneratedfrom: a)m-recurrence(basicalgorithm): Startingfrom Q0 −1=2andQ1 −1=2,generate fPm −1=2;Qm −1=2;0mMg (forlarge xbetteruse series for PM −1=2insteadoftheCF). b)evaluate PM +1=2(andPM−1 +1=2): QM −1=2,H=QM +1=2=QM −1=2!QM 1=2 QM −1=2,QM 1=2,PM −1=2!PM 1=2(fromtheWronskian). c)m-recurrence(backward): PM 1=2,PM−1 1=2!Pm 1=2,0mM Thena)+c)give Pm 1=2with0mM. d)n-recurrence(forward): Pm −1=2,Pm +1=20mM!Pm n−1=20mM,0nN. e)CF+ Wronskiantoget Q0 N+1=2,Q0 N−1=2fromP0 N+1=2,P0 N−1=2. Andsimilarlyweget Q1 N+1=2,Q1 N−1=2. f)m-recurrence(forward): Q0 N1=2,Q1 N1=2!Qm N1=2,0mM. g)n-recurrence(backward): Qm N+1=2,Qm N−1=2,0mM!Qm n1=2,0mM,0nN. This algorithm for toroidal harmonics evaluates and stores in each run first and second kindtoroidalharmonics. 4.5. NumericaltestsandCPUtimes. Inallcases,ouralgorithmshavebeenextensively tested in orderto controlthe accuracyand CPU times [8, 9, 24]. In the case of slow conver- gence of the CF, we have replaced it with series or asymptotic expansions. For paraboliccylinder functions of integer orders a, we have compared our code with other existing code byTaubmann[25],concludingthatourcode[23] clearlysupersedesit. In double precision arithmetic, the codes for PSH and OSH were shown to reach an accuracy of 10 −15in their ranges of validity. For PCF and TH the accuracy was better than 10−12. In tables 1, 2, 3 and 4 we show the CPU time spent on a HP715/100 computer for ourroutinestoevaluateprolatespheroidalharmonics(DPROH),oblatespheroidalharmonics (DOBLH),toroidalharmonics(DTORH3)and paraboliccylinderfunctionsof integerorders (DINPCF),respectively. Routine DOBLH uses the algorithm explicitly shown and DPROH use a similar one. Both routines evaluate, for a fixed order m, the first and second kind corresponding ALF’s of orders n=0;1; :::; N(withNchosen at will). Routine DINPCF also uses our basic algorithmtoevaluateintegerorderPCF’softhefirstandsecondkindsoforders n=0; :::; N. RoutineDTORH3usesthealgorithmdescribedinsection4.4. ETNA Kent State University [email protected] 144 Legendre functions and parabolic cylinder functions xMNMax 103CPU-t 103CPU-t N=NMax N=1 0 1:01 54393 10:87s 0:27s 501983 5:06s 0:16s 1:151411 3:51s 0:14s 50 709 1:88s 0:15s 10:5208 0:54s 0:07s 50 92 0:26s 0:12s 1000:5 79 0:20s 0:06s 50 14 0:06s 0:11s Table 1. Subroutine DPROH. CPU times (in 10−3s) for several values of xand M. The demanded precision is EPS= 10−15.NMaxaccounts for the maximum order that can be reached for an overflow 1.d+280. xMNMax 103CPU-t 103CPU-t N=NMax N=1 0 0:01560803 153:57s 4:11s 5015472 46:73s 4:76s 0:156211 15:48s 0:46s 502651 6:87s 0:50s 1:5712 1:77s 0:10s 50 365 1:00s 0:15s 10:5208 0:51s 0:07s 50 92 0:32s 0:12s 1000:5 79 0:23s 0:06s 50 14 0:13s 0:12s Table 2. Subroutine DOBLH. CPU times (in 10−3s) for several values of xand M. The demanded precision is EPS= 10−15.NMaxaccounts for the maximum order that can be reached for an overflow 1.d+280. x 102CPU-t. M50 ;N50 1:1 1:41s 10: 1:43s 100: 1:41s 1000: 1:41s Table 3. Subroutine DTORH3. CPU times (in 1/100 s) in evaluating fPm n−1=2;Qm n−1=2gfor several values of x,m=0 ;1; :::; 50and n=0 ;1; :::; 50. The demanded precision isEPS= 10−12. ETNA Kent State University [email protected] J.Segura and A.Gil 145 NMax 103CPU-t NMax 103CPU-t 103CPU-t xMODE=0 N=NMaxMODE=1 N=NMax N=1 0 0:1 276 0:54s 276 0:54s 0:09s 1:0 271 1:63s 271 1:67s 0:12s 2:0 265 0:91s 265 0:88s 0:28s 10: 222 0:46s 230 0:48s 0:13s 1000: 93 0:17s 0:06s Table 4. Subroutine DINPCF. CPU times (in 10−3s) for several values of xand N. EPS= 10−15. NMaxaccounts for the maximum order thatcan be reached for an overflow 1.d+280. 5. Conclusions. A set of algorithms to evaluate oblate and prolate spheroidal harmon- ics, toroidal harmonics and parabolic cylinder functions of integer and half-integer orders, have been described. These functions appear in a large variety of fields. Prolate spheroidalandoblatespheroidalharmonicsappearinthesolutionofthepotentialproblemsfordomains bounded by spheroids while toroidal harmonics appear in domains bounded by tori. On the other hand, parabolic cylinder functionsof integer and half-integerorders are used in statis-tical thermodynamics, lattice field theory, etc. In spite of their importance, there were very fewcodesinthenumericallibrariestoevaluatethem. Ouralgorithmsandtheresultingcodes fill thisgap. Acknowledgments. The authorswish to thank the Departamentode F´ ısica Te´orica (U. Valencia)theuseoftheircomputerfacilities. Theauthorsalsowouldliketoacknowledgethe hospitalityofCWI (Amsterdam)wherethisworkwasconcluded. REFERENCES [1] M. A BRAMOWITZ &I. STEGUN(eds.),Handbook of Mathematical Functions , Dover Publications, Inc. (1972). [2] D.E.A MOS,Algorithm644 : AportablepackageforBesselfunctionsofacomplexargumentandnonnegative order, ACM Trans. Math. Software, 12 (1986), p. 265. [3] C. B REZINSKI (ed.),Continued Fractions and Pad´ e approximants , North Holland (1990), p. 142. [4] B.C. C ARLSON,E . M .N OTIS, Algorithm 577 : Algorithms for incomplete elliptic integrals ,A C MT r a n s . Math. Software, vol. 7, (1981), p. 398. [5] J. L. DELYRA,S.K.F OONG,T .E.G ALLIVAN ,Finite Lattice Systems with True Critical Behaviour , Phys. Rev., D46 (1992), p. 1643. [6] W. G AUTSCHI , Algorithm 259 : Legendre functions for arguments larger than one , Communications of the ACM, vol. 8, no. 8 (1965), p. 488. [7] W. G AUSTCHI ,Computational aspects of three-term recurrence relations , SIAM Rev.,9 (1967), p.24. [8] A. G IL,J.SEGURA,Evaluation of Legendre functions ofargument greater than one , Comput. Phys. Comm., 105 (1997), p. 273. [9] A. G IL,J.SEGURA,A code to evaluate Prolate and Oblate Spheroidal Harmonics , Comput. Phys. Comm., 108 (1998), p. 267. [10] S. K UYUCAK,M.HOYLES,S.H.C HUNG,Analytical solutions of Poisson’s equation for realistic geometri- cal shapes of membrane ion channels , Biophys. J., 74 (1998), p. 22. [11] W.P. L ATHAM,R . W .R EDDING,On the calculation of the parabolic cylinder functions , J. Computational Phys.,16 (1974), p.66. [12] N.N.L EBEDEV,Special functions &their applications , Dover Publications, Inc. (1972). [13] D. W. L OZIER AND F. W. J. O LVER,Numerical Evaluation of Special Functions , in Mathematics of Com- putation 1943-1993: A Half Century of Computational Mathematics, Walter Gautschi, ed., ProceedingofSymposia in Applied Mathematics, 48 (British Columbia U., 1993). [14] D. W. L OZIER,Towards a Revised NBS Handbook of Mathematical Functions , Preprint NISTIR 6072. http://www.nist.gov/ DigitalMathLib. [15] G. M AINO,E.MENAPACE ,A.VENTURA,Computation of parabolic cylinder functions by means of a Tri- comi expansion . J. Comput. Phys.,40 (1981), p. 294. ETNA Kent State University [email protected] 146 Legendre functions and parabolic cylinder functions [16] B.P H.VANMILLIGEN,A.L´OPEZFRAGUAS,Expansion of vacuum magnetic fields in toroidal harmonics . Comput. Phys. Comm., 81 (1994), p. 74. [17] F.W.J.O LVER,J.M.S MITH,Associated Legendre functions on the cut . J. Comput. Phys., 51 (1983), p 502. [18] N.S. P ANNU,R .J .R EAD,Improved Structure Refinement Through Maximum Likehood , Acta Cryst., A52 (1996), p. 659. [19] W.H. P RESS,S.A.T EUKOLSKY ,W .T.V ETTERLINGAND B.P.FLANNERY ,Numerical Recipes in Fortran , Cambridge University Press, 1992, Chap. 6. [20] R.W. R EDDING,W . P .L ATHAM,On the calculation of the parabolic cylinder functions. II. The function V(a; x). J. Comput. Phys., 20 (1976), p. 256. [21] Z. S CHULTEN ,R . G .G ORDON, D.G.M. A NDERSON ,A numerical algorithm for the evaluation of Weber parabolic cylinder functions U(a; x),V(a; x), and W(a; x), J. Comput. Phys.,42 (1981), p. 213. [22] J. S EGURA,P .FERN´ANDEZ DE CORDOBA,YU.L. RATIS,A code to evaluate modified Bessel functions based on the continued fraction method , Comput. Phys. Comm.,105 (1997), p. 263. [23] J. S EGURA,A.GIL,Parabolic Cylinder Functions of integer and half-integer orders for non-negative argu- ments. Comput. Phys. Comm., 115 (1998), p. 69. [24] J. S EGURA,A.GIL,Evaluation of Toroidal Harmonics , submitted for publication in Comput. Phys. Comm. [25] G. T AUBMANN ,ParabolicCylinder Functions U(n,x)fornatural nandpositive x , Comput. Phys.Comm.,69 (1992), p. 415. [26] I.J.T HOMPSON ,A.R.B ARNETT,CoulombandBesselfunctionsofcomplexargumentsandorder ,J.Comput. Phys.,64 (1986), p.490. [27] I.J. T HOMPSON ,A.R.B ARNETT,Modified Bessel functions I(z),K(z)of real order and complex argu- ment, to selected accuracy , Comput. Phys. Comm.,47 (1987), p. 245. [28] N.M.T EMME,Special functions: Anintroduction totheClassical FunctionsofMathematical Physics ,Wiley, 1996.