How to do a fast numerical computation of an oscillatory integral including $\operatorname{HeunC}$ function using Mathematica?

Question

I am trying to numerically compute the following integral in Mathematica

$$\int_{1}^{1000} dx \,(x+2)e^{-2Ia(x+2)}\operatorname{HeunC}[-4Ib,-4Ib,1+4Ib,1,-4Ib,-x/2]$$

where $\operatorname{HeunC}$ is the confluent Heun function. Mathematica can numerically compute this integral for the known values of $a$ and $b$. There are two problems here: When I increase the value of $a$ and $b$, the integral becomes highly oscillatory and it takes a long time to be computed. Also, even for the small value of $a$ and $b$, Mathematica spends a long time to compute this integral when the upper bound is $1000$, but it is pretty fast when I decrease the upper bound. Could you please help me with this problem?

Answer 1

It looks like HeunC was introduced in Mathematica only recently and may have implementation issues.

For shortness let's introduce $$ F(x) = \operatorname{HeunC}(-4ib, -4ib, 1+4ib, 1,-4ib; -x/2) $$

Evaluating it with MachinePrecision gives completely wrong results (obtained with https://www.wolframcloud.com/):

F[b_,x_] = HeunC[-4I b,-4I b,1+4I b,1,-4I b,-x/2];
F[10,500.]
-2.333438484634705*^4219+4.27023838819976*^4218 I
F[10,500.`1000]
0.00398375247818748....-0.0000498007848027014.... I
F[1,500.]
-1.84935*^270-1.79429*^266 I
F[1,500.`1000]
0.0039527529500429...-0.00049881082215... I

Plotting the function for small $x$ suggests that its absolute value monotonicaly decreases and never exceeds 1, so values with $10^{4219}$ magnitude are obviously wrong. I guess this leads to slow (and simply wrong) numerical integration in your case.

One of the possible solutions is to implement HeunC by yourself. We know that HeunC is a solution to $$ y''(z) + \left(\frac{\gamma}{z} + \frac{\delta}{z-1} + \epsilon\right) y'(z) + \frac{\alpha z - q}{z(z-1)} y(z) = 0\\ y(0) = 1 $$ In your case $$ q = \alpha = \epsilon = -4ib\\ \gamma = 1 + 4ib\\ \delta = 1\\ z = -\frac{x}{2} $$ Rewriting the equation in terms of $F$, $b$ and $x$ gives: $$ F''(x) + \left( \frac{1+4ib}{x} + \frac{1}{x+2} + 2ib \right) F'(x) + \frac{2ib}{x} F(x) = 0\\ F'(0) = 1 $$

The problem is that you cannot just take this equation and ask Mathematica to solve it numerically. This is because

• this second order equation has only one initial condition,
• the point $x = 0$ is a singularity.

The former is a consequence of the latter. The usual solution of the problem is to step a bit away from the $x=0$ and to start numerical integration from $x=x_0 > 0$. Anyway, to start integration from $x = x_0$ one needs initial conditions, both for $F(x_0)$ and $F'(x_0)$.

One possibility is to ask Mathematica for them, since we know that $F(x)$ may be expressed in terms of HeunC. Another approach is to use series expansion. The series expansion can be derived from the expansion for HeunC, or may be obtained quite easily from the equation itself (Frobenius method).

Assume that $F(x) = \sum_{n=0}^\infty c_n x^n$. We already know that $c_0 = 1$. Multiplying the equation with $x$ and taking limit when $x \to 0$ yields $$ (1 + 4 ib) F'(0) + 2 i b F(0) = 0 \implies c_1 = F'(0) = -\frac{2ib}{1 + 4ib} $$ Now let's multiply the equation with $x(x+2)$: $$ (x^2+2x) F''(x) + \left( 2ib x^2 + (2 + 8ib) x + 2 + 8 ib \right) F'(x) + (2ib x + 4) F(x) = 0 $$ Let's now plug $F(x) = \sum_n c_n x^n$ and collect the coefficients before $x^n$ term: $$ n(n-1) c_n + 2(n+1) n c_{n+1} + 2ib (n-1) c_{n-1} + (2 + 8ib) n c_n + (2 + 8ib) (n+1) c_{n+1} + 2ib c_{n-1} + 4ib c_n = 0 $$ $$ \left[2(n+1) n + (2 + 8ib) (n+1)\right]c_{n+1} + \left[n(n-1) + (2 + 8ib) n + 4ib\right] c_n + \left[2ib(n-1) + 2ib\right] c_{n-1} = 0,\\ 2 (n+1) (1 + 4ib + n) c_{n+1} + \left[n(n+1) + 4ib(1+2n)\right] c_n + 2ib n c_{n-1} = 0. $$ This is a three-terms recurrence for $c_n$ from which we may obtain as many $c_n$ as we need. $$ c_{n+1} = -\frac{\left[n(n+1) + 4ib(1+2n)\right] c_n + 2ib n c_{n-1}}{2 (n+1) (1 + 4ib + n)} $$ From the theory of the confluent Heun's function it is known that the radius of convergence of its series at $z=0$ is $R = 1$. Since we use $x/2$ as its argument, the radius of convergence of $\sum_{n=0}^\infty c_n x^n$ should be equal to $2$. Numerical experiment agrees with it:

Here's a plot of $\frac{|c_n|}{|c_{n+1}|}$ which illustrate that the limit for $n \to \infty$ is 2. We may take $x_0 = 1$ and expect that series $\sum_{n=0}^\infty c_n x_0^n$ converge fast enough for plactical usage. Estimating truncation error with the first neglected term gives us that $M=50$ terms are enough for double precision (16 digits).

Now we can obtain $F(x)$ numerically by solving the problem $$ F''(x) + \left( \frac{1+4ib}{x} + \frac{1}{x+2} + 2ib \right) F'(x) + \frac{2ib}{x} F(x) = 0, \quad x \in [1, 1000]\\ F(1) = \sum_{n=0}^{50} c_n\\ F'(1) = \sum_{n=1}^{50} n c_n $$

To minimize computational problems arising from numerical integration of an (piecewise) numerical solution I suggest evaluating integral in the same ODE problem. This is achieved by expanding the ODE problem with an auxiliary equation $$ J'(x) = (x+2)e^{-2Ia(x+2)} F(x)\\ J(1) = 0. $$ It is clear that $$ J(x) = \int_1^x (t+2)e^{-2Ia(t+2)} F(t) dt $$ and we're interested in $J(1000)$ only. Combining the problems into single ODE guarantees that ODE integration grid will be the same as quadrature integration grid. This cannot be achieved by separate NDSolve + NIntegrate calls.

Here's an implementation of this method in a notebook. But be warned - WorkingPrecision stongly affects the result. Using MachinePrecision is very fast, but the relative error is about $10^{-5}$. Using MachinePrecision=20 improves the result by an order of magnitude or two, but timing is much worse.