On the nonoscillatory phase function for Legendre's differential equation

TL;DR

This paper introduces a nonoscillatory phase function for Legendre's differential equation, enabling more stable and efficient evaluation of Legendre functions of large orders through asymptotic expansion and numerical validation.

Contribution

It presents the first derivation of a nonoscillatory phase function for Legendre's equation and develops an asymptotic expansion for large-order Legendre functions.

Findings

The phase function is nonoscillatory and well-behaved.

The asymptotic expansion improves evaluation accuracy for large orders.

Numerical experiments confirm the effectiveness of the proposed methods.

Abstract

We express a certain complex-valued solution of Legendre's differential equation as the product of an oscillatory exponential function and an integral involving only nonoscillatory elementary functions. By calculating the logarithmic derivative of this solution, we show that Legendre's differential equation admits a nonoscillatory phase function. Moreover, we derive from our expression an asymptotic expansion useful for evaluating Legendre functions of the first and second kinds of large orders, as well as the derivative of the nonoscillatory phase function. Numerical experiments demonstrating the properties of our asymptotic expansion are presented.

Tables

Table 1. Table 1: The relative accuracy of the expansion ( 68 ) as a function of ν 𝜈 \nu and N 𝑁 N . Here, the expansion was evaluated using double precision arithmetic.

$\nu$	$N=2$	$N=3$	$N=4$	$N=5$	$N=6$
${10}^2$	$1.55\times {10}^{-06}$	$5.30\times {10}^{-08}$	$1.48\times {10}^{-09}$	$6.05\times {10}^{-11}$	$1.17\times {10}^{-11}$
${10}^2\pi$	$5.02\times {10}^{-08}$	$5.00\times {10}^{-10}$	$2.84\times {10}^{-12}$	$6.49\times {10}^{-14}$	$5.63\times {10}^{-14}$
${10}^3$	$1.55\times {10}^{-09}$	$4.74\times {10}^{-12}$	$2.09\times {10}^{-13}$	$2.09\times {10}^{-13}$	$2.09\times {10}^{-13}$
${10}^3\pi$	$5.02\times {10}^{-11}$	$1.16\times {10}^{-12}$	$1.16\times {10}^{-12}$	$1.16\times {10}^{-12}$	$1.16\times {10}^{-12}$
${10}^4$	$2.46\times {10}^{-12}$	$1.90\times {10}^{-12}$	$1.90\times {10}^{-12}$	$1.90\times {10}^{-12}$	$1.90\times {10}^{-12}$
${10}^4\pi$	$6.70\times {10}^{-12}$	$6.70\times {10}^{-12}$	$6.70\times {10}^{-12}$	$6.70\times {10}^{-12}$	$6.70\times {10}^{-12}$
${10}^5$	$2.17\times {10}^{-11}$	$2.17\times {10}^{-11}$	$2.17\times {10}^{-11}$	$2.17\times {10}^{-11}$	$2.17\times {10}^{-11}$
${10}^6$	$2.15\times {10}^{-10}$	$2.15\times {10}^{-10}$	$2.15\times {10}^{-10}$	$2.15\times {10}^{-10}$	$2.15\times {10}^{-10}$
${10}^7$	$2.00\times {10}^{-09}$	$2.00\times {10}^{-09}$	$2.00\times {10}^{-09}$	$2.00\times {10}^{-09}$	$2.00\times {10}^{-09}$
${10}^8$	$2.33\times {10}^{-08}$	$2.33\times {10}^{-08}$	$2.33\times {10}^{-08}$	$2.33\times {10}^{-08}$	$2.33\times {10}^{-08}$
${10}^9$	$2.15\times {10}^{-07}$	$2.15\times {10}^{-07}$	$2.15\times {10}^{-07}$	$2.15\times {10}^{-07}$	$2.15\times {10}^{-07}$

Table 2. Table 2: The relative accuracy of the expansion ( 68 ) as a function of ν 𝜈 \nu and N 𝑁 N . Here, the expansion was evaluated using quadruple precision arithmetic.

$\nu$	$N=2$	$N=3$	$N=4$	$N=5$	$N=6$
${10}^2$	$1.55\times {10}^{-06}$	$5.30\times {10}^{-08}$	$1.48\times {10}^{-09}$	$6.05\times {10}^{-11}$	$1.17\times {10}^{-11}$
${10}^2\pi$	$5.02\times {10}^{-08}$	$5.00\times {10}^{-10}$	$2.84\times {10}^{-12}$	$3.28\times {10}^{-14}$	$3.15\times {10}^{-15}$
${10}^3$	$1.55\times {10}^{-09}$	$4.74\times {10}^{-12}$	$7.06\times {10}^{-15}$	$3.94\times {10}^{-16}$	$4.10\times {10}^{-16}$
${10}^3\pi$	$5.02\times {10}^{-11}$	$4.83\times {10}^{-14}$	$3.72\times {10}^{-16}$	$3.58\times {10}^{-16}$	$3.61\times {10}^{-16}$
${10}^4$	$1.55\times {10}^{-12}$	$4.80\times {10}^{-16}$	$4.94\times {10}^{-16}$	$4.70\times {10}^{-16}$	$4.73\times {10}^{-16}$
${10}^4\pi$	$5.02\times {10}^{-14}$	$4.95\times {10}^{-16}$	$5.06\times {10}^{-16}$	$4.84\times {10}^{-16}$	$4.87\times {10}^{-16}$
${10}^5$	$1.53\times {10}^{-15}$	$3.41\times {10}^{-16}$	$3.50\times {10}^{-16}$	$3.28\times {10}^{-16}$	$3.31\times {10}^{-16}$
${10}^6$	$5.10\times {10}^{-16}$	$4.86\times {10}^{-16}$	$4.95\times {10}^{-16}$	$4.73\times {10}^{-16}$	$4.76\times {10}^{-16}$
${10}^7$	$4.85\times {10}^{-16}$	$4.57\times {10}^{-16}$	$4.67\times {10}^{-16}$	$4.42\times {10}^{-16}$	$4.45\times {10}^{-16}$
${10}^8$	$4.43\times {10}^{-16}$	$4.14\times {10}^{-16}$	$4.24\times {10}^{-16}$	$3.98\times {10}^{-16}$	$4.02\times {10}^{-16}$
${10}^9$	$4.91\times {10}^{-16}$	$4.63\times {10}^{-16}$	$4.72\times {10}^{-16}$	$4.48\times {10}^{-16}$	$4.51\times {10}^{-16}$

Table 3. Table 3: A comparison of the average time (in seconds) required to evaluate ( 68 ) for various values of N 𝑁 N and ν 𝜈 \nu with the time requires to evaluate the first 17 17 17 terms of Stieltjes’ expansion ( 42 ). Note that Stieltjes’ formula yields only the value of P ν subscript 𝑃 𝜈 P_{\nu} while ( 68 ) yields both P ν subscript 𝑃 𝜈 P_{\nu} and Q ν subscript 𝑄 𝜈 Q_{\nu} . Moreover, unlike ( 68 ), Stieltjes’ is not accurate for arguments close to the singular points of Legendre’s differential equation.

$\nu$	$N=3$	$N=4$	$N=5$	$N=6$	Stieltjes’ formula
${10}^2$	$2.36\times {10}^{-06}$	$2.93\times {10}^{-06}$	$3.92\times {10}^{-06}$	$4.27\times {10}^{-06}$	$1.57\times {10}^{-06}$
${10}^2\pi$	$2.35\times {10}^{-06}$	$3.02\times {10}^{-06}$	$3.98\times {10}^{-06}$	$4.31\times {10}^{-06}$	$1.71\times {10}^{-06}$
${10}^3$	$2.38\times {10}^{-06}$	$2.98\times {10}^{-06}$	$3.95\times {10}^{-06}$	$4.31\times {10}^{-06}$	$1.58\times {10}^{-06}$
${10}^3\pi$	$2.38\times {10}^{-06}$	$2.96\times {10}^{-06}$	$3.93\times {10}^{-06}$	$4.26\times {10}^{-06}$	$1.64\times {10}^{-06}$
${10}^4$	$2.38\times {10}^{-06}$	$3.00\times {10}^{-06}$	$3.94\times {10}^{-06}$	$4.41\times {10}^{-06}$	$1.55\times {10}^{-06}$
${10}^4\pi$	$2.33\times {10}^{-06}$	$2.97\times {10}^{-06}$	$3.95\times {10}^{-06}$	$4.26\times {10}^{-06}$	$1.56\times {10}^{-06}$
${10}^5$	$2.43\times {10}^{-06}$	$2.95\times {10}^{-06}$	$3.93\times {10}^{-06}$	$4.31\times {10}^{-06}$	$1.63\times {10}^{-06}$
${10}^6$	$2.33\times {10}^{-06}$	$2.94\times {10}^{-06}$	$3.93\times {10}^{-06}$	$4.23\times {10}^{-06}$	$1.55\times {10}^{-06}$
${10}^7$	$2.44\times {10}^{-06}$	$2.98\times {10}^{-06}$	$3.96\times {10}^{-06}$	$4.37\times {10}^{-06}$	$1.54\times {10}^{-06}$
${10}^8$	$2.36\times {10}^{-06}$	$3.01\times {10}^{-06}$	$3.99\times {10}^{-06}$	$4.23\times {10}^{-06}$	$1.61\times {10}^{-06}$
${10}^9$	$2.39\times {10}^{-06}$	$2.99\times {10}^{-06}$	$3.94\times {10}^{-06}$	$4.38\times {10}^{-06}$	$1.67\times {10}^{-06}$

Table 4. Table 4: The relative accuracy with which the derivative of the nonoscillatory phase function for Legendre’s differential equation is evaluated via the expansions of Section  5 . These calculations were performed using double precision arithmetic.

$\nu$	$N=2$	$N=3$	$N=4$	$N=5$	$N=6$
${10}^2$	$4.87\times {10}^{-07}$	$2.07\times {10}^{-08}$	$7.36\times {10}^{-09}$	$7.04\times {10}^{-09}$	$7.03\times {10}^{-09}$
${10}^2\pi$	$1.53\times {10}^{-08}$	$1.64\times {10}^{-10}$	$1.60\times {10}^{-10}$	$1.60\times {10}^{-10}$	$1.60\times {10}^{-10}$
${10}^3$	$4.78\times {10}^{-10}$	$1.22\times {10}^{-12}$	$2.42\times {10}^{-15}$	$1.27\times {10}^{-15}$	$1.41\times {10}^{-15}$
${10}^3\pi$	$1.54\times {10}^{-11}$	$1.29\times {10}^{-14}$	$1.07\times {10}^{-15}$	$1.27\times {10}^{-15}$	$1.36\times {10}^{-15}$
${10}^4$	$4.78\times {10}^{-13}$	$1.36\times {10}^{-15}$	$1.36\times {10}^{-15}$	$1.48\times {10}^{-15}$	$1.36\times {10}^{-15}$
${10}^4\pi$	$1.56\times {10}^{-14}$	$1.08\times {10}^{-15}$	$1.26\times {10}^{-15}$	$1.08\times {10}^{-15}$	$1.34\times {10}^{-15}$
${10}^5$	$1.10\times {10}^{-15}$	$9.65\times {10}^{-16}$	$1.38\times {10}^{-15}$	$1.14\times {10}^{-15}$	$1.20\times {10}^{-15}$
${10}^6$	$1.30\times {10}^{-15}$	$1.30\times {10}^{-15}$	$1.19\times {10}^{-15}$	$1.30\times {10}^{-15}$	$1.33\times {10}^{-15}$
${10}^7$	$1.22\times {10}^{-15}$	$1.22\times {10}^{-15}$	$1.22\times {10}^{-15}$	$1.70\times {10}^{-15}$	$1.62\times {10}^{-15}$
${10}^8$	$1.56\times {10}^{-15}$	$1.42\times {10}^{-15}$	$1.45\times {10}^{-15}$	$1.56\times {10}^{-15}$	$1.42\times {10}^{-15}$
${10}^9$	$1.12\times {10}^{-15}$	$1.34\times {10}^{-15}$	$1.34\times {10}^{-15}$	$1.41\times {10}^{-15}$	$1.38\times {10}^{-15}$

Equations

$${\psi }^{\prime \prime }(z)-\frac{2z}{1-z^2}{\psi }^{\prime }(z)+\frac{\nu (\nu +1)}{1-z^2}\psi (z)=0$$

$$y^{\prime \prime }(x)+q(x)y(x)=0\text{mbox}foralla<x<b,$$

$$\frac{\cos (\alpha (x))}{\sqrt{\left|{\alpha }^{\prime }(x)\right|}}$$

$$\frac{\sin (\alpha (x))}{\sqrt{\left|{\alpha }^{\prime }(x)\right|}}$$

$${\eta }^{\prime \prime }(x)+p_0(x){\eta }^{\prime }(x)+q_0(x)\eta (x)=0\text{mbox}foralla<x<b,$$

$$q(x)=q_0(x)-\frac{1}{4}(p_0(x){)}^2-\frac{1}{2}{p}_0^{\prime }(x).$$

$$y(x)=\frac{\eta (x)}{\sqrt{W(x)}},$$

$$W(x)=\exp \left(-{\int }_a^x{p}_0(u)du\right).$$

$$u(x)=\sqrt{\frac{W(x)}{\left|{\alpha }^{\prime }(x)\right|}}\cos (\alpha (x))$$

$$v(x)=\sqrt{\frac{W(x)}{\left|{\alpha }^{\prime }(x)\right|}}\sin (\alpha (x))$$

$$y^{\prime \prime }(x)+{\lambda }^{2}q(x)y(x)=0\text{mbox}foralla<x<b$$

$$u_1(x)=\frac{1}{\sqrt{\lambda }(q(x){)}^{1/4}}\cos \left(\lambda {\int }_a^x\sqrt{q(u)}du\right)+O\left(\frac{1}{\lambda }\right)$$

$$v_1(x)=\frac{1}{\sqrt{\lambda }(q(x){)}^{1/4}}\sin \left(\lambda {\int }_a^x\sqrt{q(u)}du\right)+O\left(\frac{1}{\lambda }\right).$$

$${\alpha }_1(x)=\lambda {\int }_a^x\sqrt{q(u)}du$$

$$u_{\infty }(x)=\frac{\cos ({\alpha }_{\infty }(x))}{\sqrt{\left|{\alpha }_{\infty }^{\prime }(x)\right|}}+O\left(\exp (-\mu \lambda )\right)$$

$$v_{\infty }(x)=\frac{\sin ({\alpha }_{\infty }(x))}{\sqrt{\left|{\alpha }_{\infty }^{\prime }(x)\right|}}+O\left(\exp (-\mu \lambda )\right).$$

$$y^{\prime \prime }(x)+q(x,\lambda )y(x)=0$$

$$\tan (\alpha (x))=\frac{v(x)}{u(x)}$$

$$\mathrm{cotan}(\alpha (x))=\frac{u(x)}{v(x)}$$

$${\alpha }^{\prime }(x)=\frac{u(x)v^{\prime }(x)-u^{\prime }(x)v(x)}{(u(x){)}^2+(v(x){)}^2}$$

$$\psi (x)=u(x)+iv(x),$$

$$\frac{{\psi }^{\prime }(x)}{\psi (x)}=\frac{u(x)u^{\prime }(x)+v(x)v^{\prime }(x)}{(u(x){)}^2+(v(x){)}^2}+i\frac{u(x)v^{\prime }(x)-u^{\prime }(x)v(x)}{(u(x){)}^2+(v(x){)}^2}$$

$$\text{Hypergeometric}21a,bcz=\sum _{n=0}^{\infty }\frac{{\left(a\right)}_n{\left(b\right)}_n}{{\left(c\right)}_n}\frac{z^n}{\Gamma (n+1)}.$$

$$z(z+1)(z+2)\cdots (z+n-1).$$

$$P_{\nu }(z)=\text{Hypergeometric}21-\nu ,\nu +11\frac{1-z}{2}$$

$$Q_{\nu }(z)=(2z{)}^{-\nu -1}\sqrt{\pi }\frac{\Gamma (\nu +1)}{\Gamma (\nu +3/2)}\text{Hypergeometric}21\frac{\nu }{2}+1,\frac{\nu }{2}+\frac{1}{2}\nu +\frac{3}{2}\frac{1}{z^2}$$

$$Q_{\nu }(x)=\underset{y\to 0^{+}}{\lim }\frac{Q_{\nu }(x+iy)+Q_{\nu }(x-iy)}{2}.$$

$$\underset{y\ge y_0}{\sup }\left|F(x+iy)\right|<p_{y_0}(x)\text{mbox}forallx\in \mathbb{R}.$$

$$\underset{y\to 0^{+}}{\lim }{Q}_{\nu }(x+iy)={\phi }_{\nu }(x),$$

$${\phi }_{\nu }(x)=Q_{\nu }(x)-i\frac{\pi }{2}{P}_{\nu }(x)$$

Full text

1 Introduction

The Legendre functions of degree $\nu$ — that is, the solutions of the second order linear ordinary differential equation

[TABLE]

— appear in numerous contexts in physics and applied mathematics. For instance, they arise when certain partial differential equations are solved via separation of variables, they are often used to represent smooth functions defined on bounded intervals, and their roots are the nodes of Gauss-Legendre quadrature rules.

In this article, we give a constructive proof of the existence of a nonoscillatory phase function for Legendre’s differential equation. A smooth function $\alpha:(a,b)\to\mathbb{R}$ is a phase function for the second order linear ordinary differential equation

[TABLE]

where $q$ is smooth and positive on $(a,b)$ (it can have zeros or singularities at the endpoints of the interval), if $\alpha^{\prime}$ does not vanish on $(a,b)$ and the functions

[TABLE]

and

[TABLE]

form a basis in the space of solutions of (2). When $q$ is positive the solutions of (2) are oscillatory, whereas when $q$ is negative the solutions of (2) are roughly increasing and decreasing exponentials. In the latter case, (3) and (4) are not appropriate mechanisms for representing them, hence the requirement that $q$ is positive. The ostensibly more general second order linear ordinary differential equation

[TABLE]

where the coefficients $p_{0}$, $q_{0}$ are real-valued and smooth on $(a,b)$, can be reduced to the form (2) with the function $q$ given by

[TABLE]

This is accomplished by letting

[TABLE]

where $W$ is defined via

[TABLE]

Assuming that the function $q$ defined by (6) is positive, we say $\alpha$ is a phase function for (5) if it is a phase function for the transformed equation (2). In this case, the functions $u$, $v$ defined via the formulas

[TABLE]

and

[TABLE]

form a basis in the space of solutions of the original equation (5). We note that according to Abel’s identity, the Wronskian of any pair of solutions of (5) is a constant multiple of the function $W$, and that the Wronskian of the pair defined by (9) and (10) is $W$.

It has long been known that certain second order differential equations admit nonoscillatory phase functions, and that a large class of such equations do so in an asymptotic sense. Indeed, under mild conditions on $q$, the second order equation

[TABLE]

admits a basis $\{u_{1},v_{1}\}$ of solutions such that

[TABLE]

and

[TABLE]

This standard result can be found in [15] and [10], among many other sources. Assuming that the coefficient $q$ is nonoscillatory, the phase function

[TABLE]

associated with (12) and (13) obviously is as well. Expressions of this type are generally known as WKB approximations, and higher order analogues of them can be constructed. The resulting formulas, however, are quite unwieldy — indeed, the $n^{th}$ order WKB expansion, which achieves $O\left(\lambda^{-n}\right)$ accuracy, involves complicated combinations of the derivatives of $q$ of orders [math] through $n-1$ and various noninteger powers of $q$. Perhaps more importantly, efficient and accurate numerical methods for their computation are lacking.

It is shown in [13, 4] that under mild conditions on $q$, there exist a positive constant $\mu$, a function $\alpha_{\infty}$ which is roughly as oscillatory as $q$, and a basis $\{u_{\infty},v_{\infty}\}$ in the space of solutions of the differential equation (11) such that

[TABLE]

and

[TABLE]

The function $\alpha_{\infty}$ can be represented using various series expansions (such as piecewise polynomial expansions) the number of terms of which do not depend on the parameter $\lambda$. In order words, nonoscillatory phase functions represent the solutions of equations of the form (11) with $O(\exp(-\mu\lambda)$ accuracy using $O(1)$-term expansions. Moreover, a reliable and efficient numerical algorithm for the computation of $\alpha_{\infty}$ which runs in time independent of $\lambda$ and only requires knowledge of the values of $q$ on the interval $(a,b)$ is introduced in [3]. Much like standard WKB estimates, the results of [13, 4, 3] easily apply to large class of differential equations whose coefficients are allowed to vary with $\lambda$ — that is, equations of the more general form

[TABLE]

— assuming that $q$ satisfy certain innocuous conditions independent of $\lambda$.

The framework of [13, 4, 3] applies to (1) and it can be used to, among other things, evaluate Legendre functions of large orders and their zeros in time independent of degree. However, in cases like Legendre’s differential equation where an exact nonoscillatory phase function exists, it is advantageous to derive formulas for it which are as explicit as possible. To that end, in this article, we derive an integral representation of a particular solution $\psi_{\nu}$ of Legendre’s differential equation. This is done in Section 3, after certain preliminaries are dispensed with in Section 2.8. In Section 4, our integral representation of $\psi_{\nu}$ is used to prove the existence of a nonoscillatory phase function $\alpha_{\nu}$ for Legendre’s differential equation. In Section 5, we derive an asymptotic formula for Legendre functions of large degrees, and for the derivative of the nonoscillatory phase function $\alpha_{\nu}$. Numerical experiments demonstrating the properties of these expansions are discussed in Section 6. We conclude this article with brief remarks in Section 7.

2 Preliminaries

2.1 An elementary observation regarding phase functions

Suppose that $u$, $v$ form a basis in the space of solutions of (5). Then, since $u$ and $v$ cannot simultaneously vanish, (9) and (10) hold for all $x\in(a,b)$ if and only if

[TABLE]

for all $x\in(a,b)$ such that $u(x)\neq 0$ and

[TABLE]

for all $x\in(a,b)$ such that $v(x)\neq 0$. By differentiating (18) and (19) we find that (9) and (10) are satisfied if and only if

[TABLE]

for all $a<x<b$. In particular, if $u$, $v$ form a basis in the space of solutions of (5), then any function $\alpha$ whose derivative satisfies (20) is a phase function for (5).

We will frequently make use of the following elementary theorem:

Theorem 1**.**

If $u$, $v$ are real-valued solutions of (5) and

[TABLE]

then the imaginary part of the logarithmic derivative $\psi^{\prime}(x)/\psi(x)$ of $\psi$ is the derivative of a phase function for (5)

It is easily established by observing that

[TABLE]

and comparing (20) and (22).

2.2 Gauss’ hypergeometric function

For any complex-valued parameters $a$, $b$, and $c$ such that $c\neq 0,-1,-2,\ldots$ we use $\Hypergeometric{2}{1}{a,b}{c}{z}$ to denote Gauss’ hypergeometric function, which is defined by the formula

[TABLE]

Here, $(z)_{n}$ is the Pochhammer symbol

[TABLE]

Since the series (23) converges absolutely for all $|z|<1$, $\Hypergeometric{2}{1}{a,b}{c}{z}$ is an analytic function on the unit disk in the complex plane.

2.3 Legendre functions

We refer to the function

[TABLE]

obtained by employing the Frobenius (series solution) method to construct a solution of (1) which is analytic in a neighborhood of the regular singular point $z=1$ as the Legendre function of the first kind of degree $\nu$. When $\nu$ is an integer, the series (25) terminates and $P_{\nu}$ is entire; otherwise, it can be analytically continued to the cut plane $\mathbb{C}\setminus\left(-\infty,-1\right]$ (see, for instance, Section 2.10 of [8]). Similarly, we refer to the solution

[TABLE]

of (1) as the Legendre function of the second kind of degree $\nu$. It can be analytically continued to the cut plane $\mathbb{C}\setminus\left(-\infty,1\right]$. Following the standard convention (see, for example, [5, 8, 11]), we define $Q_{\nu}(x)$ for $x\in(-1,1)$ via the formula

[TABLE]

It can be verified readily that while $Q_{\nu}$ is not an element of any of the classical Hardy spaces on the upper half of the complex plane $\mathbb{H}$, it is an element of the space $H^{+}$ of functions which are analytic on $\mathbb{H}$ and uniformly bounded by a polynomial on every half plane of the form $\operatorname{Im}(z)\geq y_{0}>0$. That is, a function $F$ analytic on $\mathbb{H}$ is an element of $H^{+}$ if for every $y_{0}>0$ there exists a polynomial $p_{y_{0}}$ such that

[TABLE]

The space $H^{+}$ generalizes the classical Hardy spaces $H^{p}$ and many of the useful properties of Hardy spaces still apply in this setting. For instance, if $1\leq p<\infty$ then the elements of $H^{p}$ are the Fourier transforms of $L^{p}$ functions supported on the half line $[0,\infty)$, while those of $H^{+}$ are the Fourier transforms of distributions supported on $[0,\infty)$. The classical theory of Hardy spaces is discussed in the well-known textbook [14] and the properties of $H^{+}$ are described in some detail in [1].

According to Formula (9) in Section 3.4 of [8] (see also 8.732.5 in [11]),

[TABLE]

where

[TABLE]

for all $-1<x<1$. Using the fact that $Q_{\nu}\in H^{+}$ and the observation that it has no zeros in the upper half of the complex plane (this is established, for instance, in [16] via the argument principle), it can be shown that the logarithmic derivative of the function $\varphi_{\nu}$ is nonoscillatory. It then follows from Theorem 1 that Legendre’s differential equation admits a nonoscillatory phase function. The details of this argument are beyond the scope of this article, but we mention it because of its strong relation to the approach taken here.

2.4 A transformation of Legendre’s differential equation

It is easy to verify that the functions $P_{\nu}(\cos(\theta))$ and $Q_{\nu}(\cos(\theta))$ satisfy the second order differential equation

[TABLE]

on the interval $(0,\pi/2)$. We will, by a slight abuse of terminology, refer to both (1) and (31) as Legendre’s differential equation. In all cases, though, it will be clear from the context which specific form of Legendre’s differential equation is intended.

2.5 Bessel functions

For nonnegative integers $n$, the Bessel function of the first kind of order $n$ is the entire function defined via the series

[TABLE]

while the Bessel function of the second kind of order $n$ is

[TABLE]

where $\psi$ denotes the logarithmic derivative of the gamma function. We use the principal branch of the logarithm in (33), so that $Y_{n}$ is defined on the cut plane $\mathbb{C}\setminus\left(-\infty,0\right]$.

The Hankel function of the first kind of order $n$ is

[TABLE]

and we let $K_{n}$ denote the modified Bessel function of the first kind of order $n$. For complex-valued $z$ such that $-\pi<\arg(z)\leq\pi/2$ , $K_{n}$ is defined in terms of $H_{n}$ via the formula

[TABLE]

It is easily verified that

[TABLE]

2.6 The Lipschitz-Hankel integrals

The formulas

[TABLE]

and

[TABLE]

hold when $0<\theta<\pi/2$ and $\operatorname{Re}(\nu)>-1$. Derivations of (37) and (38) can be found in Chapter 13 of [18]; they can also be found as Formulas 6.628(1) and 6.628(2) in [11].

2.7 The Laplace transforms of $\left(x^{2}+bx\right)^{-1/2}$ and $x^{\nu}$

For complex numbers $z$ and $b$ such that $\operatorname{Re}(z)>0$ and $|\arg(b)|<\pi$,

[TABLE]

(as usual, we take the principal branch of the square root function). A careful derivation of (39) can be found in Section 7.3.4 of [9] and it appears as Formula 3.383(8) in [11]. By combining (39) with (35) we see that

[TABLE]

whenever $\operatorname{Re}(z)>0$, $-\pi<\arg(bz)\leq\pi/2$ and $|\arg(b)|<\pi$.

It is abundantly well known that for any complex numbers $z$ and $\nu$ such that $\operatorname{Re}(\nu)>-1$ and $\operatorname{Re}(z)>0$,

[TABLE]

this identity appears, for instance, as Formula 3.381(4) in [11].

2.8 Steiltjes’ asymptotic formula for Legendre functions

For all positive real $\nu$ and all $0<\theta<\pi/2$,

[TABLE]

where

[TABLE]

[TABLE]

and

[TABLE]

This approximation was introduced by Stieltjes; a derivation of the error bound (45) can be found in Chapter 8 of [17], among many other sources. In [2], an efficient method for computing the coefficients $C_{\nu,k}$ is suggested. The coefficient in the first term is given by

[TABLE]

This ratio of gamma functions can be approximated via a series in powers of $1/\nu$; however, it can be calculated more efficiently by observing that the related function

[TABLE]

admits an expansion in powers of $1/x^{2}$. In particular, the $7$-term expansion

[TABLE]

gives roughly double precision accuracy for all $x>10$ (see the discussion in [2]). The coefficient $C_{\nu,0}$ is related to $\tau$ through the formula

[TABLE]

The subsequent coefficients are obtained through the recurrence relation

[TABLE]

3 An integral representation of a particular solution of Legendre’s equation

For $0<\theta<\pi/2$ and $\operatorname{Re}(\nu)>-1$, we define the function $\psi_{\nu}$ via

[TABLE]

Obviously, $\psi_{\nu}$ is a solution of Legendre’s equation (31) and it is related to the function $\varphi_{\nu}$ defined in Section 2.3 via the formula

[TABLE]

In this section, we derive an integral representation of $\psi_{\nu}$ involving only elementary functions. We prefer $\psi_{\nu}$ over $\varphi_{\nu}$ because it leads to an asymptotic formula involving the Hankel functions of the first kind of order [math] instead of the Hankel function of the second kind of order [math].

From (37), (38) and the definition (34) of the Hankel function of the first kind of order [math], we see that

[TABLE]

whenever $0<\theta<\pi/2$ and $\operatorname{Re}(\nu)>-1$. We rearrange (53) as

[TABLE]

let $\beta=\exp(i\theta)\sin(\theta)$ and introduce the new variable $w=\exp(-i\theta)t$ to obtain

[TABLE]

Now by letting $b=2i$ and $z=\beta w$ in (40), we find that

[TABLE]

Formula (40) holds since $0<\arg(z)<\pi/2$ and $\arg(b)=-\pi/2$. Inserting (56) into (55) yields

[TABLE]

Since $\operatorname{Re}(1+\beta x)>0$, we may use Formula (41) to evaluate the integral with respect to $w$ in (57) whenever $\operatorname{Re}(\nu)>-1$. In this way we conclude that

[TABLE]

where $\beta=\sin(\theta)\exp(i\theta)$, for all $\operatorname{Re}(\nu)>-1$ and $0<\theta<\pi/2$.

4 A nonoscillatory phase function for Legendre’s equation

Formula (58) expresses $\psi_{\nu}$ as the product of an oscillating exponential function and an integral involving only nonoscillatory elementary functions. It follows that $\psi_{\nu}$ gives rise to a nonoscillatory phase function for Legendre’s differential equation. To see this, we first define $\sigma_{\nu}$ via the formula

[TABLE]

so that

[TABLE]

Next, we observe that

[TABLE]

Now, we combine (59), (60) and (61) to obtain

[TABLE]

According to Theorem 1, the imaginary part of $\psi_{\nu}^{\prime}/\psi_{\nu}$ is the derivative of a phase function $\alpha_{\nu}$ for Legendre’s differential equation (31). We conclude form Formula (62) that $\alpha_{\nu}$ and $\alpha_{\nu}^{\prime}$ are nonoscillatory functions. Figure 1 contains plots of the imaginary part of the logarithmic derivative of the function $\psi_{\nu}$ when $\nu=100\pi$ and when $\nu=10^{6}$.

5 An asymptotic formula for Legendre functions of large degrees

In this section, we derive an asymptotic expansion for the Legendre functions of large degrees. Our starting point is the formula

[TABLE]

which is obtained from (58) by introducing the new variable $\tau=\beta x=\sin(\theta)\exp(i\theta)x$. We replace the function

[TABLE]

in (63) with a sum of the form

[TABLE]

where $p=\nu+1$, $q=\sqrt{p}$, and the coefficients $a_{0},a_{1},b_{1},a_{2},b_{2},\ldots,a_{N},b_{N}$ are chosen so the power series expansions of $f$ and $g$ around [math] agree to order $2N$. That is, we require that the system of $2N+1$ linear equations

[TABLE]

in the $2N+1$ variables $a_{0},a_{1},b_{1},a_{2},b_{2},\ldots,a_{N},b_{N}$ be satisfied. By so doing we obtain the approximation

[TABLE]

Now by applying Formula (40) to (67) we conclude that

[TABLE]

where $\beta=\sin(\theta)\exp(i\theta)$, $p=\nu+1$, and $q=\sqrt{p}$. The use of Formula (40) is justified so long as $\operatorname{Re}(p)>N^{2}$ and $0<\theta<\pi/2$ (the second condition ensures that $|\arg(2i\beta)|<\pi$). The linear system (66) can be solved easily using a computer algebra system. Our Mathematica script for doing so appears in an appendix of this paper. A second appendix lists the coefficients when $N=2$, $N=3$, $N=4$, $N=5$ and $N=6$.

Error bounds for the formula (68) can be derived quite easily using standard techniques (the monograph [15], for example, contains many similar estimates). We omit them here, though, because our point is not that (68) is especially accurate, but rather that its form is quite different from most widely used asymptotic expansions for Legendre functions in that it expresses $\psi_{\nu}$ as the product of an oscillatory exponential function and a sum of nonoscillatory functions rather than as a sum of oscillatory functions. We do, however, report on extensive numerical experiments which assess the performance of Formula (68) in Section 6 of this article. That the functions

[TABLE]

appearing in (68) are nonoscillatory is obvious from the formula

[TABLE]

obtained by letting $b=-2i$ in (40). Stieltjes’ formula (42), which approximates $P_{\nu}$ via a sum of cosines which oscillate at frequencies on the order of $\nu$, furnishes an example of the latter type of expansion, as does Olver’s asymptotic expansion

[TABLE]

where $u=n+1/2$. The definition of the coefficients $A_{k}$ and $B_{k}$ as well as a derivation of (71) can be found in Chapter 12 of [15]. Yet another example is given by Dunster’s convergent expansions for Legendre functions [7], which have a form similar to (71).

Among other things, the form of the expansion (68) is conducive to computing the derivative of the nonoscillatory phase function $\alpha_{\nu}$ associated with $\psi_{\nu}$. To see this, we rewrite (68) as

[TABLE]

where the function $S$ is defined by

[TABLE]

and differentiate both sides in order to obtain

[TABLE]

From (36) we see that

[TABLE]

By combining (74) and (68) we can evaluate the logarithmic derivative $\psi_{\nu}^{\prime}/\psi_{\nu}$ of $\psi_{\nu}$ and hence the derivative of a nonoscillatory phase function for Legendre’s differential equation, which is the imaginary part of this ratio.

6 Numerical experiments

In this section, we describe numerical experiments which were conducted to assess the performance of the asymptotic expansions of Section 5. Our code was written in Fortran and compiled with the GNU Fortran Compiler version 5.2.1. The calculations were carried out on a laptop equipped with an Intel i7-5600U processor running at 2.60GHz.

6.1 The accuracy of the expansion (68) as a function of $N$ and $\nu$

We measured the accuracy of the expansion (68) for various values of $N$ and $\nu$. For each pair of values of $N$ and $\nu$ considered, we evaluated (68) at a collection of $1,{,}000$ points on the interval $(0,\pi/2)$. The first $500$ points were drawn at random from the uniform distribution on the interval $(0,\pi/2)$, while the remaining points were constructed by drawing $500$ points from the uniform distribution on the interval $(0,1)$ and applying the mapping $t\to\exp(-36t)$ to each of them. In this way, we ensured that the accuracy of (68) was tested near the singular point of Legendre’s equation which occurs when $\theta=0$.

The results are reported in Tables 1 and 2. Table 1 gives the maximum relative error in the value of $\psi_{\nu}$ which was observed when these calculations were carried out in double precision arithmetic as a function of $\nu$ and $N$, and Table 2 reports the maximum relative error in the value of $\psi_{\nu}$ which was observed when they were performed using quadruple precision arithmetic as a function of $\nu$ and $N$. Reference values were computed by running the algorithm of [3] in quadruple precision arithmetic. Like the asymptotic expansion (68), the algorithm of [3] is capable of achieving high accuracy near the singular point of Legendre’s equation whereas Stieltjes’ expansion (42) and the well-known three term recurrence relations are inaccurate when $\theta$ is near 0. See [2, 12], though, for the derivation of asymptotic expansions of Legendre functions which are accurate for $\theta$ near 0 and [15, 6, 7] for expansions of Legendre functions which are accurate for all $\theta\in(0,\pi/2)$.

The routine we used to evaluate the Hankel function of order [math] achieves roughly double precision accuracy, even when executed using quadruple precision arithmetic. Thus the minimum error achieved when the computations were performed using quadruple precision arithmetic was on the order of $10^{-16}$ (see Table 2).

We also observe that relative accuracy was lost as $\nu$ increases when (68) is evaluated using double precision arithmetic. This loss of precision is unsurprising and consistent with the condition number of the evaluation of the highly oscillatory functions $P_{\nu}$ and $Q_{\nu}$. Indeed, evaluating $P_{\nu}(x)$ and $Q_{\nu}(x)$ is analogous to calculating the values of $\cos(\nu x)$ and $\sin(\nu x)$. If $\cos(\nu x)$ and $\sin(\nu x)$ can be evaluated with accuracy on the order of $\epsilon$, then so can the value of $\mbox{mod}(\nu,1)$. This clearly limits the precision with which $\cos(\nu x)$ and $\sin(\nu x)$ can be evaluated using finite precision arithmetic. A similar argument applies to $P_{\nu}(x)$ and $Q_{\nu}(x)$. The only exception is when the value of $\mbox{mod}(\nu,1)$ is known to high precision (for instance, when $\nu$ is an integer). Then $\cos(\nu x)$ and $\sin(\nu x)$ can be calculated with comparable precision, as can $P_{\nu}$ and $Q_{\nu}$ (via the well-known three term recurrence relations, for instance).

6.2 The accuracy of the expansion (68) as a function of $\theta$

In order to measure the accuracy of the expansion (68) as a function of $\theta$, we sampled a collection of $1,{,}000$ points on the interval $(0,\pi/2)$ using the same method as in Section 6.1 and then evaluated (68) at each of the chosen values of $\theta$ with $\nu=1,{,}000$ and $N=2$, $N=3$ and $N=4$. The base-$10$ logarithm of the relative error in each of the resulting values was computed and a plot showing these errors appears in Figure 2. Once again, reference values were computed by running the algorithm of [3] in quadruple precision arithmetic. We observe that the error in Formula (68) is remarkably uniform as a function of $\theta$ — it is nearly constant until $\theta$ nears the singular point at [math], at which point it decreases slightly.

6.3 The speed of the expansion (68) as function of $N$

Next, we measured the time required to evaluate $\psi_{\nu}$ using the expansion (68). In particular, for several pairs of values of $N$ and $\nu$, we evaluated (68) at a collection of $1,{,}000$ points drawn from the uniform distribution on the interval $(0,\pi/2)$. We also applied the same procedure to Stieltjes’ expansion (42) with $M=16$. Table 3 gives the average time required to evaluate (68) and (42).

We note that while the time required to evaluate (68) is slightly larger than the time required to evaluate Stieltjes’ formula (42), the expansion of Section 5 gives both the value of $P_{\nu}$ and that of $Q_{\nu}$ while Stieltjes’ formula gives only the value of $P_{\nu}$. Moreover, unlike Stieltjes’ formula, our asymptotic approximation is accurate for $\theta$ near [math].

6.4 Evaluation of the derivative of a phase function for Legendre’s equation

In this last experiment, we measured the accuracy achieved when Formulas (68) and (74) are combined in order to evaluate the derivative of the nonoscillatory phase function $\alpha_{\nu}$ for Legendre’s equation associated with the function $\psi_{\nu}$. More specifically, for each of several pairs of values of $N$ and $\nu$, we evaluated (68) and (74) at a collection of $1,{,}000$ points in the interval $(0,\pi/2)$ and used these quantities to calculate $\alpha_{\nu}^{\prime}$ at each of these points. The points were chosen as in Section 6.1. These experiments were performed using double precision arithmetic. The obtained values of $\alpha_{\nu}^{\prime}$ were compared with reference values computed by running the algorithm of [3] using quadruple precision arithmetic.

The results are reported in Table 4. For each pair of values of $\nu$ and $N$, it lists the maximum relative error in $\alpha_{\nu}^{\prime}$ which was observed. We note that, unlike the experiments of Section 6.1, near double precision accuracy was obtained by performing the calculations in double precision arithmetic. This is not surprising since the condition number of the evaluation of the nonoscillatory function $\alpha_{\nu}^{\prime}$ is small and not dependent on $\nu$.

7 Conclusions

Nonoscillatory phase functions are powerful analytic and numerical tools. Among other things, explicit formulas for them can be used to efficiently and accurately evaluate special functions, their zeros and to apply special function transforms.

Here, we factored a particular solution of Legendre’s differential equation as the product of an oscillatory exponential function and an integral involving only nonoscillatory elementary functions. By so doing, we proved the existence of a nonoscillatory phase function and derived an asymptotic formula of an unusual type. Specifically, our formula represents the oscillatory functions $P_{\nu}$ and $Q_{\nu}$ as the product of an oscillatory exponential function and a sum of nonoscillatory functions.

We will report on the use of the results of the paper to apply the Legendre transform rapidly and on generalizations of this work to the case of associated Legendre functions, prolate spheroidal wave functions and other related special functions at a later date.

8 Acknowledgments

James Bremer was supported by National Science Foundation grant DMS-1418723. Vladimir Rokhlin was supported by the Office of Naval Research under contract N00014-16-1-2123 and by the Air Force Scientific Research Office under contract FA9550-16-1-0175.

Appendix A Mathematica source code

Below is a Mathematica script for generating the coefficients $a_{0},a_{1},b_{1},\ldots,a_{n},b_{n}$ in the expansion (68).

(**************************************************************************************) ClearAll["Global‘"]

(* Construct the (2n+1)-term expansion ) n=3; f[t_]=1/(1+t)^(q^2); g[t_]=a0Exp[-q^2*t]+Sum[a[k]Exp[-(q^2+kq)*t]+b[k]Exp[-(q^2-kq)*t],{k,1,n}];

eqs={}; vars={}; h0[t_]=f[t]-g[t]; For[k=0,k<=2*n,k++,eqs=Simplify[Join[eqs,{h0[0]==0}]];h0[t_]=h0’[t]]; For[k=1,k<=n,k++,vars=Join[vars,{a[k],b[k]}]]; vars=Join[vars,{a0}]; sol=Simplify[Solve[eqs,vars][[1]],p>0]

(* Write the TEX version of the expansion to the file tmp2 *) OpenWrite["tmp2"]; WriteString["tmp2","\begin{equation}\n\begin{aligned}\n"]; WriteString["tmp2","a_0 &="ToString[TeXForm[a0/.sol]],", \\ \n"]; For[k=1,k<=n,k++, WriteString["tmp2","a_",k," &="ToString[TeXForm[a[k]/.sol]],", \\ \n"]]; For[k=1,k<=n-1,k++, WriteString["tmp2","b_",k," &="ToString[TeXForm[b[k]/.sol]],", \\ \n"]]; WriteString["tmp2","b_",k," &="ToString[TeXForm[b[k]/.sol]],".\n"]; WriteString["tmp2","\end{aligned}\n\end{equation}"]; Close["tmp2"];

(* Write the FORTRAN version of the expansion to tmp ) OpenWrite["tmp"]; WriteString["tmp","a0 = ",ToString[FortranForm[N[Expand[a0/.sol],36]]],"\n"]; For[k=1,k<=n,k++, WriteString["tmp","a(",k,") ="ToString[FortranForm[N[Expand[a[k]/.sol],36]]],"\n"]]; For[k=1,k<=n-1,k++, WriteString["tmp","b(",k,") ="ToString[FortranForm[N[Expand[b[k]/.sol],36]]],"\n"]]; WriteString["tmp","b(",k,") ="ToString[FortranForm[N[Expand[b[k]/.sol],36]]],"\n"]; Close["tmp"]; (**************************************************************************************)

Appendix B The coefficients in the expansion (68)

When $N=2$, the coefficients are

[TABLE]

When $N=3$, they are

[TABLE]

When $N=4$:

[TABLE]

When $N=5$:

[TABLE]

When $N=6$:

[TABLE]