# Airy Functions (and Their First Derivative) Compute the Airy functions $Ai$ or $Bi$ or their first derivatives, $d/dz Ai(z)$ and $d/dz Bi(z)$. The Airy functions are solutions of the differential equation [REMOVE_ME]$$ w'' = z w [REMOVE_ME_2]$$ for $w(z)$, and are related to each other and to the (modified) Bessel functions via (many identities, see [https://dlmf.nist.gov/9.6](https://dlmf.nist.gov/9.6)), e.g., if c("$%$", "$\n$", "$ zta := (2/3) z * sqrt(z) = 2/3 z^(3/2)$"), [REMOVE_ME]$$ Ai(z) = \pi^{-1}\sqrt{z/3}K_{1/3}(\zeta) =\frac{1}{3}\sqrt{z}\left(I_{-1/3}(\zeta) - I_{1/3}(\zeta)\right),%Ai(z) = 1/pi sqrt(z/3) K[1/3](zta) =1/3 sqrt(z) * (I[-1/3](zta) - I[1/3](zta)), [REMOVE_ME_2]$$ and [REMOVE_ME]$$ Bi(z) = \sqrt{z/3} \left(I_{-1/3}(\zeta) + I_{1/3}(\zeta)\right).%Bi(z) = sqrt(z/3) * (I[-1/3](zta) + I[1/3](zta)). [REMOVE_ME_2]$$ 1.1 ## Description Compute the Airy functions $Ai$ or $Bi$ or their first derivatives, $d/dz Ai(z)$ and $d/dz Bi(z)$. The Airy functions are solutions of the differential equation $$ w'' = z w $$ for $w(z)$, and are related to each other and to the (modified) Bessel functions via (many identities, see [https://dlmf.nist.gov/9.6](https://dlmf.nist.gov/9.6)), e.g., if c("$%$", "$\n$", "$ zta := (2/3) z * sqrt(z) = 2/3 z^(3/2)$"), $$ Ai(z) = \pi^{-1}\sqrt{z/3}K_{1/3}(\zeta) =\frac{1}{3}\sqrt{z}\left(I_{-1/3}(\zeta) - I_{1/3}(\zeta)\right),%Ai(z) = 1/pi sqrt(z/3) K[1/3](zta) =1/3 sqrt(z) * (I[-1/3](zta) - I[1/3](zta)), $$ and $$ Bi(z) = \sqrt{z/3} \left(I_{-1/3}(\zeta) + I_{1/3}(\zeta)\right).%Bi(z) = sqrt(z/3) * (I[-1/3](zta) + I[1/3](zta)). $$ ```r AiryA(z, deriv = 0, expon.scaled = FALSE, verbose = 0) AiryB(z, deriv = 0, expon.scaled = FALSE, verbose = 0) ``` ## Arguments - `z`: complex or numeric vector. - `deriv`: order of derivative; must be 0 or 1. - `expon.scaled`: logical indicating if the result should be scaled by an exponential factor (typically to avoid under- or over-flow). - `verbose`: integer defaulting to 0, indicating the level of verbosity notably from C code. ## Details By default, when `expon.scaled` is false, `AiryA()` computes the complex Airy function $Ai(z)$ or its derivative $d/dz Ai(z)$ on `deriv=0` or `deriv=1` respectively. When `expon.scaled` is true, it returns $exp(zta)*Ai(z)$ or $exp(zta)* d/dz Ai(z)$, effectively removing the exponential decay in $-pi/3 < arg(z) < pi/3$ and the exponential growth in $pi/3 < abs(arg(z)) < pi$, where $zta=(2/3)*z*sqrt(z)$, and $arg(z) = $`Arg(z)`. While the Airy functions $Ai(z)$ and $d/dz Ai(z)$ are analytic in the whole $z$ plane, the corresponding scaled functions (for `expon.scaled=TRUE`) have a cut along the negative real axis. By default, when `expon.scaled` is false, `AiryB()` computes the complex Airy function $Bi(z)$ or its derivative $d/dz Bi(z)$ on `deriv=0` or `deriv=1` respectively. When `expon.scaled` is true, it returns $exp(-abs(Re(zta)))*Bi(z)$ or c("$%$", "$\n$", "$ exp(-abs(Re(zta)))* dBi(z)/dz$"), to remove the exponential behavior in both the left and right half planes where, as above, $zta=(2/3)*z*sqrt(z)$. ## Returns a complex or numeric vector of the same length (and class) as `z`. ## References see `BesselJ`; notably for many results the Digital Library of Mathematical Functions (DLMF), Chapter 9 **Airy and Related Functions** at [https://dlmf.nist.gov/9](https://dlmf.nist.gov/9). ## Author(s) Donald E. Amos, Sandia National Laboratories, wrote the original fortran code. Martin Maechler did the interface. ## See Also `BesselI` etc; the Hankel functions `Hankel`. The CRAN package [list("Rmpfr")](https://CRAN.R-project.org/package=Rmpfr) has `Ai(x)` for arbitrary precise "mpfr"-numbers `x`. ## Examples ```r ## The AiryA() := Ai() function ------------- curve(AiryA, -20, 100, n=1001) curve(AiryA, -1, 100, n=1011, log="y") -> Aix curve(AiryA(x, expon.scaled=TRUE), -1, 50, n=1001) ## Numerically "proving" the 1st identity above : z <- Aix$x; i <- z > 0; head(z <- z[i <- z > 0]) Aix <- Aix$y[i]; zeta <- 2/3*z*sqrt(z) stopifnot(all.equal(Aix, 1/pi * sqrt(z/3)* BesselK(zeta, nu = 1/3), tol = 4e-15)) # 64b Lnx: 7.9e-16; 32b Win: 1.8e-15 ## This gives many warnings (248 on nb-mm4, F24) about lost accuracy, but on Windows takes ~ 4 sec: curve(AiryA(x, expon.scaled=TRUE), 1, 10000, n=1001, log="xy") ## The AiryB() := Bi() function ------------- curve(AiryB, -20, 2, n=1001); abline(h=0,v=0, col="gray",lty=2) curve(AiryB, -1, 20, n=1001, log = "y") # exponential growth (x > 0) curve(AiryB(x,expon.scaled=TRUE), -1, 20, n=1001) curve(AiryB(x,expon.scaled=TRUE), 1, 10000, n=1001, log="x") ```