arXiv:1607.08269v3 [math.CA] 24 Feb 2017
Noname manuscript No.
(will be inserted by the editor)
Computation of asymptotic expansions of turning point
problems via Cauchy’s integral formula: Bessel functions
T. M. Dunster · A. Gil · J. Segura
Received: date / Accepted: date
Abstract Linear second order differ ential equations having a large real pa-
rameter and turning point in the co mplex plane are considered. Classical
asymptotic expansions fo r solutions involve the Airy function and its deriva-
tive, along with two infinite series, the c oefficients of which are usually difficult
to compute. By considering the series as asymptotic expansions for two ex-
plicitly defined analytic functions, Cauchy’s integral formula is employed to
compute the coe fficient functions to high order of accuracy. The method em-
ploys a certain exponential form of Liouville-Green expansions for solutions of
the differential equation, as well as for the Airy function. We illustrate the use
The authors acknowledge support from Ministerio de Econom´ıa y Competitividad, project
MTM2015-67142-P (MINECO/FEDER, UE). A.G. and J.S. acknowledge support from Min-
isterio de Econom´ıa y Competitividad, project MTM2012-34787. A.G. acknowledges the
Fulbright/MEC Program f or support during her stay at SDSU. J. S. acknowledges the Sal-
vador de Madariaga Program for support during his stay at SDSU.
T. M. Dunster
Department of Mathematics and Statistics
San Diego State University. 5500 Campanile Dri ve San Diego, CA, USA.
E-mail: mdunster@mail.sdsu.edu
A. Gil
Departamento de Matem´atica Aplicada y CC. de la Computaci´on.
ETSI Caminos. Universidad de Cantabria. 39005-Santander, Spain.
and
Department of Mathematics and Statistics
San Diego State University. 5500 Campanile Dri ve San Diego, CA, USA.
E-mail: amparo.gil@unican.es
J. Segura
Departamento de Matem´aticas, Estad´ıstica y Computaci´on.
Universidad de Cantabria. 39005-Santander, Spain.
and
Department of Mathematics and Statistics
San Diego State University. 5500 Campanile Dri ve San Diego, CA, USA.
E-mail: segurajj@unican.es
2 T.M. Dunster, A. Gil, J. Segura
of the method with the high accuracy c omputation of Airy-type expansions of
Bessel functions of complex argument.
Keywords Turning point problems · Asymptotic expansions · Bessel
functions · Numerical computation
Mathematics Subject Classification (2010) MSC 34E05 · 34E20 ·
33C10 · 33F05
1 Introdu ction
In this paper we study linear second order differe ntial equations having a
simple turning point. Spe c ific ally, we consider the differential equation
d
2
w/dz
2
=
u
2
f(z) + g(z)
w, (1.1)
where u is positive and large, f(z) has a simple zero (turning point) at z = z
0
(say), and f (z) and g(z) are analytic in an unbounded domain containing the
turning po int.
This is a classical problem, with applications to numerous special functions.
To obtain asymptotic solutions, the Liouville tra nsformation
2
3
ζ
3/2
= ±
Z
z
z
0
f
1/2
(t)dt, W = ζ
−1/4
f
1/4
(z)w, (1.2)
is applied, where either sign in front of the integral c an be chosen. As a result
we transform (1.1) to the form
d
2
W/dζ
2
=
u
2
ζ + ψ(ζ)
W, (1.3)
where
ψ(ζ) =
5
16ζ
2
+
4f(z)f
′′
(z) −5f
′2
(z)
ζ
16f
3
(z)
+
ζg(z)
f(z)
. (1.4)
The lower integration limit in (1.2) ensures that the turning point z = z
0
of
(1.1) is mapped to the turning point ζ = 0 of (1.3). Throughout this paper we
shall assume that this turning point is bounded away from any other turning
points or sing ularities of (1.1), equivalently ψ(ζ) is analytic for 0 ≤ |ζ| < R
for some positive R which is independent of u.
When the turning point z
0
is real and f (z) is re al on a real interval around
z
0
, the sign in (1.2) is usually chosen in such a way that the new variable ζ is
real when z is real and in a neighborhood of the turning point.
From [9, Chap. 11, Theorem 9.1] we obta in solutions having the following
asymptotic expansions in terms of Airy functions
W
2n+1,j
(u, ζ) = Ai
j
u
2/3
ζ
n
X
s=0
A
s
(ζ)
u
2s
+
Ai
′
j
u
2/3
ζ
u
4/3
n−1
X
s=0
B
s
(ζ)
u
2s
+ ε
2n+1,j
(u, ζ),
(1.5)
Computation of asymptotic solutions of turning point problems: B essel functions 3
0
−1
2π/3
S
S
S
1
Fig. 1 The sectors S
j
in the complex plane
for j = 0, ±1. Here Ai
j
(u
2/3
ζ) = Ai(u
2/3
ζe
−2πij/3
), which are the Airy func-
tions that are recessive in the sectors S
j
:=
ζ : |ar g(ζe
−2πij/3
)| ≤ π/3
(Fig.
1); se e [10, §9.2(iii)]. Also, note that Ai
′
j
(z) = dAi
j
(z)/dz = e
−2πij/3
Ai
′
(ze
−2πij/3
).
From Olver’s explicit error bounds we have
ε
2n+1,j
(u, ζ) = env
n
Ai
j
u
2/3
ζ
o
O
u
−2n−1
,
as u → ∞, for ζ lying in certain domains described in [9, Chap. 11, §9], and
which we assume to be unbounded. For a definition of the envelope function
env for Airy functions, see [10, §2 .8(iii)].
In (1.5) A
0
(ζ) is an arbitrar y non-zero constant (typically taken to be 1),
and for s = 0, 1, 2, ···, the other coefficients satisfy the recursion relations
B
s
(ζ) =
1
2ζ
1/2
Z
ζ
0
{ψ(t)A
s
(t) − A
′′
s
(t)}
dt
t
1/2
, (1.6)
and
A
s+1
(ζ) = −
1
2
B
′
s
(ζ) +
1
2
Z
ψ(ζ)B
s
(ζ)dζ. (1.7)
We remark that the lower integration limit in (1.6) must be 0 in order for
each B
s
(ζ) to be analytic at ζ = 0, whe reas in (1.7) there is no restriction in
the choice of integration cons tant. This will be of significance to us below.
In general, these coefficients are difficult to compute, primarily due to the
requirement of repeated integrations. They also show cancellations near the
turning point. For complex ζ close to 0 one can compute these coefficients in
a numerically stable way by considering power s eries expansions for the coeffi-
cients, as done in [1,2], where they are expanded in powers of ω =
√
1 − z
2
. For
4 T.M. Dunster, A. Gil, J. Segura
other computational approaches to compute the c oefficients, and in par ticular
for real values of ζ, see [12].
The purpose of this paper is to provide a more simple means of computing
a larg e number of these terms. We shall employ Cauchy’s integral formula to
do so, and our results will be valid for real a nd complex ζ lying in a bounded
(but not necessarily small) domain containing the turning point ζ = 0. O ur
approach can potentially be extended to other situations, including the cases
of simple poles [9, Chap. 12], and coalescing turning points.
We illustrate the use of the method with the high accuracy computation
of Airy-type expansions of Bessel functions of complex argument.
2 General method
We first present Liouville-Green expansions for solutions of (1.1), a certain
form of which will be required for our method. These only involve elementary
(exp onential) functions, but are not valid at the turning point. The appropriate
Liouville-Green transforma tion is given by [9, Chap. 10, §2], namely
ξ =
2
3
ζ
3/2
, V = f
1/4
(z)w. (2.1)
With these, equation (1.1) is transformed to equation
d
2
V/dξ
2
=
u
2
+ φ(ξ)
V, (2.2)
where
φ(ξ) =
4f(z)f
′′
(z) −5f
′2
(z)
16f
3
(z)
+
g(z)
f(z)
. (2.3)
The branch for the first of (2 .1) will be dependent on the solutions under
consideration, as described below. Note, as ζ completes one circuit about the
turning point ζ = 0, the variable ξ correspondingly cros ses more than one
Riemann sheet.
It turns out that solutions where asymptotic expansions appear inside ex-
ponentials are more convenient for our purposes. Specifically, from [9, Chap.
10, Ex. 2.1] we have solutions
V
±
n
(u, ξ) = exp
(
±uξ +
n−1
X
s=1
(±1)
s
E
s
(ξ)
u
s
)
+ ε
±
n
(u, ξ). (2.4)
In these, the co efficients are given by
E
s
(ξ) =
Z
F
s
(ξ)dξ (s = 1, 2, 3, ···), (2.5)
where
F
1
(ξ) =
1
2
φ(ξ), F
2
(ξ) = −
1
4
φ
′
(ξ), (2.6)
Computation of asymptotic solutions of turning point problems: B essel functions 5
and
F
s+1
(ξ) = −
1
2
F
′
s
(ξ) −
1
2
s−1
X
j=1
F
j
(ξ)F
s−j
(ξ) (s = 2, 3, ···). (2.7)
Primes are derivatives with respect to ξ. The integration constants in (2.5)
will be disc ussed below, and we find that for the odd c oefficients E
2j+1
(ξ)
(j = 0, 1, 2, ···) they must be suitably chosen.
Explicit error bounds for ε
±
n
(u, ξ), which verify the asymptotic validity of
the ex pansions (2.4), are given in [5]. In particular, for arbitrary δ > 0, under
certain conditions on ψ(ξ), we have that ε
±
n
(u, ξ) = e
±uξ
O(u
−n
) as ξ → ∞ in
certain domains Ξ
±
as described in [9, Chap. 10, §3]. These are the same as
those for the corresponding asymptotic solutions of the more common form
V ∼ exp{±uξ}
∞
X
s=0
(±1)
s
A
s
(ξ)u
−s
. (2.8)
It is the relation (2.7) that is the reason why the e xpansions (2.4) are nu-
merically advantageous: the coefficients F
s
(ξ) can a ll be determined explicitly
without resorting to integration. Furthermore, from (2.5) we observe that only
one integration is r e quired (numerical or explicit) to evaluate each E
s
(ξ), as
opposed to repeated integrals fo r computing the coefficients A
s
(ξ) in (2.8).
Remarkably, it turns out that integration is not required to evaluate the
even terms E
2j
(ξ) (j = 1, 2, 3, ···). To see this, consider the Wronskian of the
solutions V
±
n
(u, ξ) given by (2.4). Since this is a constant (by Abel’s theorem)
we infer that
u +
∞
X
j=0
F
2j+1
(ξ)
u
2j+1
exp
2
∞
X
j=1
E
2j
(ξ)
u
2j
∼ constant,
which, on taking loga rithms, yields
∞
X
j=1
E
2j
(ξ)
u
2j
∼ −
1
2
ln
1 +
∞
X
j=0
F
2j+1
(ξ)
u
2j+2
+ constant. (2.9)
We then asymptotically expand the RHS of this relation in inverse powers
of u
2
, and equate the coefficients of both sides. As a result, we find that (to
within an arbitrary additive constant in each instance)
E
2
(ξ) = −
1
2
F
1
(ξ), E
4
(ξ) =
1
4
F
2
1
(ξ) −
1
2
F
3
(ξ), (2.10)
and so on. In particular, the even coefficients E
2j
(ξ) (j = 1 , 2, 3, ···) are
explicitly given in terms of F
2k+1
(ξ) (k = 0, 1, 2, ··· , j −1) (which in turn are
given by (2.6) and (2.7)).
At this stage we consider Liouville- Green solutions of (1.3). Comparing
(1.2), (2.1) and (2.4) we obtain three asymptotic solutions W
j
(u, ζ) (j = 0, ±1)
![Fig. 5 Relative errors in the comparison of the function values obtained with the Airy-type expansion with n terms for H (1) ν (ν(1 + 0.1i)) against those obtained with Amos’ algorithm [1]. In the evaluation of the coefficients the Airy functions are also computed with [1]](/figures/fig-5-relative-errors-in-the-comparison-of-the-function-2q9usd57.webp)
![Fig. 6 Relative errors in the comparison of the function values obtained with the Airytype expansion for H (1) ν (ν(1 + 0.1i)) against those obtained with Amos’ algorithm [1]. The asymptotic expansion of the coefficients is considered over the Cauchy contour. We take n = 14 = 2m in Eqs. (5.4) and (5.5).](/figures/fig-6-relative-errors-in-the-comparison-of-the-function-11e7auhj.webp)
