High precision series solutions of differential equations: Ordinary and regular singular points of second order ODEs

A subroutine for a very-high-precision numerical solution of a class of ordinary differential equations is provided. For a given evaluation point and equation parameters the memory requirement scales linearly with precision PP, and the number of algebraic operations scales roughly linearly with PP when PP becomes sufficiently large. We discuss results from extensive tests of the code, and how one, for a given evaluation point and equation parameters, may estimate precision loss and computing time in advance.Program summaryProgram title: seriesSolveOde1Catalogue identifier: AEMW_v1_0Program summary URL: http://cpc.cs.qub.ac.uk/summaries/AEMW_v1_0.htmlProgram obtainable from: CPC Program Library, Queen’s University, Belfast, N. IrelandLicensing provisions: Standard CPC licence, http://cpc.cs.qub.ac.uk/licence/licence.htmlNo. of lines in distributed program, including test data, etc.: 991No. of bytes in distributed program, including test data, etc.: 488116Distribution format: tar.gzProgramming language: C++Computer: PC’s or higher performance computers.Operating system: Linux and MacOSRAM: Few to many megabytes (problem dependent).Classification: 2.7, 4.3External routines: CLN — Class Library for Numbers [1] built with the GNU MP library [2], and GSL — GNU Scientific Library [3] (only for time measurements).Nature of problem:The differential equation equation(1)−s2(d2dz2+1−ν+−ν−zddz+ν+ν−z2)ψ(z)+1z∑n=0Nvnznψ(z)=0, is solved numerically to very high precision. The evaluation point zz and some or all of the equation parameters may be complex numbers; some or all of them may be represented exactly in terms of rational numbers.Solution method  : The solution ψ(z)ψ(z), and optionally ψ′(z)ψ′(z), is evaluated at the point zz by executing the recursion equation(2)Am+1(z)=s−2(m+1+ν−ν+)(m+1+ν−ν−)∑n=0NVn(z)Am−n(z),equation(3)ψ(m+1)(z)=ψ(m)(z)+Am+1(z),ψ(m+1)(z)=ψ(m)(z)+Am+1(z), to sufficiently large mm. Here νν is either ν+ν+ or ν−ν−, and Vn(z)=vnzn+1. The recursion is initialized by equation(4)A−n(z)=δn0zν,for n=0,1,…,Nequation(5)ψ(0)(z)=A0(z).ψ(0)(z)=A0(z).Restrictions:   No solution is computed if z=0z=0, or s=0s=0, or if ν=ν−ν=ν− (assuming Reν+≥Reν−) with ν+−ν−ν+−ν− an integer, except when ν+−ν−=1ν+−ν−=1 and v0=0v0=0 (i.e. when zz is an ordinary point for z−ν−ψ(z)).Additional comments: The code of the main algorithm is in the file seriesSolveOde1.cc, which “#include” the file checkForBreakOde1.cc. These routines, and the programs using them, must “#include” the file seriesSolveOde1.cc.Running time:   On a Linux PC that is a few years old, at y=10 to an accuracy of P=200P=200 decimal digits, evaluating the ground state wavefunction of the anharmonic oscillator (with the eigenvalue known in advance); (cf. Eq. (6)) takes about 2 ms, and about 40 min at an accuracy of P=100000 decimal digits.References:[1] B. Haible and R.B. Kreckel, CLN — Class Library for Numbers, http://www.ginac.de/CLN/[2] T. Granlund and collaborators, GMP — The GNU Multiple Precision Arithmetic Library, http://gmplib.org/[3] M. Galassi et al., GNU Scientific Library Reference Manual (3rd Ed.), ISBN 0954612078., http://www.gnu.org/software/gsl/