Computing Airy funciona a través de su serie Taylor.

Question

He estado tratando de encontrar buenos métodos para numéricamente computación funciones de Airy, y he encontrado que para los pequeños argumentos es conveniente el uso de sus series de Taylor de las expansiones:

$$\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$ términos de la recurrencia de los grandes argumentos, parece que hay algunas numérica de problemas.

---

Mis preguntas son las siguientes:

• Puede este algoritmo se mejoró en algunos de manera obvia para una mejor convergencia y estabilidad, teniendo en cuenta la posibilidad de redondear/cancelación de los errores?

• Y ¿cuál es la mejor manera para calcular estas funciones para las grandes argumentos? Puede que yo ya uso el asintótica expansiones aparece en algunas fuentes para $|z|>7$?

---

El código R hice para este algoritmo:

#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))

---

Actualización:

Quería señalar una mejora evidente. En lugar de calcular los coeficientes, es mucho mejor utilizar atrás recurrencia. Por ejemplo, para la primera parte de la serie:

$$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$ debe ser precalculadas).

Answer 1

Esto es demasiado largo para un comentario:

A partir de los gráficos https://dlmf.nist.gov/9.3.i se puede ver que para el real $x$ sólo $\mathrm{Bi(x)}$ grandes $x>0$ puede ser evaluado con el poder de la serie de una manera estable. De lo contrario, la oscilatorio o de manera exponencial pequeñas funciones no son adecuadas para la serie de implementación. Las alternativas dependen de las funciones disponibles:

Para mi biblioteca me encontré con que la serie es mejor sólo para un pequeño $x$intervalo de $-2 \le x \le 1$ para $\mathrm{Ai}(x)$ e $-1 \le x \le 1$ para $\mathrm{Bi}(x)$ (a pesar de que suma los positivos y los negativos de los términos por separado), de lo contrario, el uso de la relación con las funciones de Bessel. Para $x>1$ la $\mathrm{Ai}$ función se evalúa como (véase el comunicado de Prensa et al, recetas Numérica, Ch. 6.7) con $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$$ y

$$\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$$

Para grandes negativo $x$ las funciones de Airy tienen expansiones asintóticas oscilante con $\cos(z + \pi/4)$ o $\sin(z + \pi/4),$ por lo tanto la fase de información se convierte en totalmente confiable para $x < -(2/\epsilon)^{2/3}$ (sobre $-0.43\cdot 10^{11}$ IEEE precisión doble) y el error relativo se incrementa fuertemente para $x$ menos que el cuadrado de la raíz ($\approx -200000$).