Computing Airy functions through their Taylor series

Question

I've been trying to find good methods for numerically computing Airy functions, and I found that for small arguments it's convenient to use their Taylor series expansions:

$$\text{Ai} (z)=\frac{1}{3^{2/3} \pi} \sum_{n=0}^{\infty} \frac{\Gamma(1/3(n+1))}{n!}(3^{1/3}z)^n \sin[(2(n+1)\pi)/3]$$

$$\text{Bi} (z)=\frac{1}{3^{1/6} \pi} \sum_{n=0}^{\infty} \frac{\Gamma(1/3(n+1))}{n!}(3^{1/3}z)^n |\sin[(2(n+1)\pi)/3]|$$

While these expansions are inconvenient for direct implementation, I've modified them in the following way. Let's define some quantities:

$$C_0=\frac{\Gamma(1/3)}{3^{1/3}} \\ D_0=\Gamma(2/3)$$

And recurrence relations:

$$C_k=\frac{C_{k-1}}{3k(3k-1)} \\ D_k=\frac{D_{k-1}}{3k(3k-2)}$$

Then the general terms will be:

$$A_k=z^{3k} \left(C_k-\frac{z}{3k+1}D_k \right) \\ B_k=z^{3k} \left(C_k+\frac{z}{3k+1}D_k \right)$$

And finally:

$$\text{Ai}(z)= \frac{3^{1/6}}{2 \pi}\sum_{k=0}^\infty A_k$$

$$\text{Bi}(z)= \frac{3^{2/3}}{2 \pi}\sum_{k=0}^\infty B_k$$

This seems to work great for (approximately) $|z|<7$ (including complex values) using around $20-25$ terms of the recurrence, for the larger arguments there seem to be some numerical issues.

---

My questions are as follows:

• Can this algorithm be improved in some obvious way for better convergence/stability, taking into account possible rounding/cancellation errors?

• And what is the best way to compute these functions for larger arguments? Can I already use the asymptotic expansions listed in some sources for $|z|>7$?

---

The R code I made for this algorithm:

#Airy functions using series
AiBi <- function(z){
g13 <- 2.678938534707747633655693; #Gamma(1/3)
g23 <- 1.354117939426400416945288; #Gamma(2/3)
Km <- 25; #Number of terms used
C <- 1:Km;
D <- C;
A <- C;
B <- C;
C[1] <- g13/3^(1/3);
D[1] <- g23;
A[1] <- C[1]-D[1]*z;
B[1] <- C[1]+D[1]*z;
k <- 1;
while (k < Km) {k <- k+1;
        C[k] <- C[k-1]/(3*k-3)/(3*k-4);
        D[k] <- D[k-1]/(3*k-3)/(3*k-5);
        A[k] <- z^(3*k-3)*(C[k]-D[k]*z/(3*k-2));
        B[k] <- z^(3*k-3)*(C[k]+D[k]*z/(3*k-2))
};
c(sum(A)*3^(1/6),sum(B)*3^(2/3))/2/pi};
AiBi(pi/exp(1))

---

Update:

Wanted to point out an obvious improvement. Instead of computing the coefficients, it's much better to use backwards recurrence. For example, for the first part of the series:

$$f(z) = \sum_{k=0}^\infty C_k z^{3k}$$

Pick some $K_m$ large enough, then define (for the sake of simplicity I changed $C_0=1$):

$$f_0=1 \\ f_{k'}=1+\frac{z^3}{3(K_m-k')(3(K_m-k')-1)}f_{k'-1}, \qquad k'=0,1,2,\dots,K_m-1$$

This should make the numerical computation faster and more stable (obviously, $z^3$ should be precomputed).

Answer 1

This is far too long for a comment:

From the graphics https://dlmf.nist.gov/9.3.i you can see that for real $x$ only $\mathrm{Bi(x)}$ for larger $x>0$ may be evaluated with the power series in a stable way. Otherwise the oscillatory or exponentially small functions are not suited for the series implementation. The alternatives depend on the functions available:

For my library I found that the series are better only for a small $x$ intervall $-2 \le x \le 1$ for $\mathrm{Ai}(x)$ and $-1 \le x \le 1$ for $\mathrm{Bi}(x)$ (although I sum the negative and positive terms separately), otherwise I use the relation to the Bessel functions. For $x>1$ the $\mathrm{Ai}$ function is evaluated as (see Press et al, Numerical recipes, Ch. 6.7) with $z=(2/3)|x|^{3/2}$ $$\mathrm{Ai}(x) = \frac{1}{\pi} \sqrt{\frac{x}{3}}\:K_{1/3}(z), \quad x > 1$$ $$\mathrm{Ai(x)}=\tfrac{1}{2} \sqrt{-x}\left( J_{1/3}(z) - \frac{1}{\sqrt{3}} Y_{1/3}(z)\right), \quad x < -2$$ and

$$\mathrm{Bi}(x)=\sqrt{x}\left( \frac{2}{\sqrt{3}} I_{1/3}(z) + \frac{1}{\pi} K_{1/3}(z)\right), \quad x > 1$$ $$\mathrm{Bi}(x)=-\tfrac{1}{2} \sqrt{-x}\left( \frac{1}{\sqrt{3}} J_{1/3}(z) + Y_{1/3}(z)\right), \quad x<-1$$

For large negative $x$ the Airy functions have asymptotic expansions oscillating with $\cos(z + \pi/4)$ or $\sin(z + \pi/4),$ therefore the phase information becomes totally unreliable for $x < -(2/\epsilon)^{2/3}$ (about $-0.43\cdot 10^{11}$ for IEEE double precision) and the relative error increases strongly for $x$ less than the square root ($\approx -200000$).