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

Abstract

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 P, and the number of algebraic operations scales roughly linearly with P when P 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 (1)−s2(d2dz2+1−ν+−ν−zddz+ν+ν−z2)ψ(z)+1z∑n=0Nvnznψ(z)=0, is solved numerically to very high precision. The evaluation point z 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), and optionally ψ′(z), is evaluated at the point z by executing the recursion (2)Am+1(z)=s−2(m+1+ν−ν+)(m+1+ν−ν−)∑n=0NVn(z)Am−n(z),(3)ψ(m+1)(z)=ψ(m)(z)+Am+1(z), to sufficiently large m. Here ν is either ν+ or ν−, and Vn(z)=vnzn+1. The recursion is initialized by (4)A−n(z)=δn0zν,for n=0,1,…,N(5)ψ(0)(z)=A0(z).Restrictions: No solution is computed if z=0, or s=0, or if ν=ν− (assuming Reν+≥Reν−) with ν+−ν− an integer, except when ν+−ν−=1 and v0=0 (i.e. when z 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=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.