How to implement a numerically efficient Airy Zeta Function

Question

Define the Airy zeta function for $n=2,3,\dots$, by $$ Z(n) := \sum_r \frac{1}{r^n}. $$ where the sum is over the zeros $r$ of the Airy function $\operatorname{Ai}$.

In Mathematica the $\operatorname{Ai}$ function is implemented as AiryAi, and the zeros of this function is implemented as AiryAiZero.

I have tried to calculate the values of $Z$ using the following: Z[n_] := Sum[1/AiryAiZero[k]^n, {k, 1, Infinity}], and then N[Z[2]] for example.

Sadly, in Mathematica 9.0 it gives $0.499$; however, the correct result is $0.531457$. For larger $n$ values $Z$ is correct; however, for only a few digits, even if I modify $MaxExtraPrecision.

The Airy zeta function MathWorld page gives a closed-form of $Z$, but then you have to implement $n$th derivatives of different Airy-related functions.

How do I implement an efficient $Z$ function?

Answer 1

I'll see if I can come up with a better routine later, but this should be satisfactory for the time being:

SetAttributes[airyZeta, Listable];
airyZeta[n_Integer?NonNegative] :=
π (SeriesCoefficient[AiryBi[\[FormalZ]]/AiryAi[\[FormalZ]],
                     {\[FormalZ], 0, n - 1}]/(3^(2/3) Gamma[1/3]^2) + 
   SeriesCoefficient[AiryAi[\[FormalZ]] AiryBi[\[FormalZ]],
                     {\[FormalZ], 0, n - 2}]/(n - 1) - 
   Sum[SeriesCoefficient[AiryBi[\[FormalZ]]/AiryAi[\[FormalZ]],
                         {\[FormalZ], 0, n - j - 1}]
       SeriesCoefficient[AiryAi[\[FormalZ]]^2, {\[FormalZ], 0, j - 1}]/j,
       {j, n - 1}, Method -> "Procedural"])

where I used SeriesCoefficient[] to evaluate the needed derivatives.

Test:

N[airyZeta[Range[2, 6]], 20]
   {0.53145723196099945287, -0.11256176121511457943, 0.039443078421238584544,
    -0.015533659376623159601, 0.0063892694802911830860}

---

By replacing some of the terms with their closed forms, I managed to produce a slightly faster airyZeta[]:

SetAttributes[airyZeta, Listable];
airyZeta[n_Integer?NonNegative] :=
π SeriesCoefficient[AiryBi[\[FormalZ]]/AiryAi[\[FormalZ]],
                    {\[FormalZ], 0, n - 1}]/(3^(2/3) Gamma[1/3]^2) +
(-1)^n 2^((2 n - 6)/3) 3^((2 n - 9)/6) Gamma[(2 n - 3)/6] Sin[π n/3]/
(Sqrt[π] (n - 1)!) -
Sum[SeriesCoefficient[AiryBi[\[FormalZ]]/AiryAi[\[FormalZ]],
                      {\[FormalZ], 0, n - j - 1}]
    Gamma[(2 j - 1)/6] Cos[(2 j - 1) π/3] 2^((2 j - 4)/3) 3^((2 j - 7)/6)/j!,
    {j, n - 1}, Method -> "Procedural"]/Sqrt[π]

(There is actually an expression for $\left.\dfrac{\mathrm d^k}{\mathrm dz^k}\dfrac{\operatorname{Bi}(z)}{\operatorname{Ai}(z)}\right|_{z=0}$ in terms of a sum of matrix Bell polynomials (BellY[]), but the performance was not too different from the original to justify the replacement.)

Answer 2

Here is a version that produces exact results and memoizes the derivatives and polynomials for efficiency:

airyZeta`dAA[n_ /; IntegerQ[n] && n >= 0] := 
 airyZeta`dAA[n] = 
  If[n == 0, AiryAi[airyZeta`z]^2, 
   D[airyZeta`dAA[n - 1], airyZeta`z]]
airyZeta`dAB[n_ /; IntegerQ[n] && n >= 0] := 
 airyZeta`dAB[n] = 
  If[n == 0, AiryAi[airyZeta`z] AiryBi[airyZeta`z], 
   D[airyZeta`dAB[n - 1], airyZeta`z]]
airyZeta`dBoA[n_ /; IntegerQ[n] && n >= 0] := 
 airyZeta`dBoA[n] = 
  If[n == 0, AiryBi[airyZeta`z]/AiryAi[airyZeta`z], 
   D[airyZeta`dBoA[n - 1], airyZeta`z]]
airyZeta`poly[n_ /; IntegerQ[n] && n >= 2] := 
 airyZeta`poly[
   n] = π/
       Gamma[n] (airyZeta`dBoA[n - 1]/(3^(2/3) Gamma[1/3]^2) + 
         airyZeta`dAB[n - 2] - 
         Sum[Binomial[n - 1, j] airyZeta`dBoA[n - 1 - j] airyZeta`dAA[
            j - 1], {j, 1, n - 1}]) /. airyZeta`z -> 0 /. 
     Gamma[2/3] -> (2 π)/(Sqrt[3] Gamma[1/3]) /. 
    Gamma[1/3] -> Sqrt[(2 π)/(3^(1/6) airyZeta`x)] // Expand
airyZeta[n_ /; IntegerQ[n] && n >= 2] := 
 airyZeta[n] = 
  airyZeta`poly[n] /. airyZeta`x -> (2 π)/(3^(1/6) Gamma[1/3]^2)
airyZetaN[n_ /; IntegerQ[n] && n >= 2] := 
 airyZetaN[n] = 
  airyZeta`poly[n] /. 
   airyZeta`x -> N[(2 π)/(3^(1/6) Gamma[1/3]^2)]

It generates the Airy Zeta polynomials, e.g.:

$$\begin{array}{c} x^2 \\ x^3-\frac{1}{2} \\ x^4-\frac{x}{3} \\ x^5-\frac{5 x^2}{12} \\ x^6-\frac{x^3}{2}+\frac{1}{20} \\ x^7-\frac{7 x^4}{12}+\frac{13 x}{180} \\ x^8-\frac{2 x^5}{3}+\frac{139 x^2}{1260} \\ x^9-\frac{3 x^6}{4}+\frac{87 x^3}{560}-\frac{1}{160} \\ x^{10}-\frac{5 x^7}{6}+\frac{209 x^4}{1008}-\frac{17 x}{1296} \\ x^{11}-\frac{11 x^8}{12}+\frac{671 x^5}{2520}-\frac{2167 x^2}{90720} \\ x^{12}-x^9+\frac{93 x^6}{280}-\frac{13 x^3}{336}+\frac{7}{8800} \\ x^{13}-\frac{13 x^{10}}{12}+\frac{2041 x^7}{5040}-\frac{10543 x^4}{181440}+\frac{9301 x}{4276800} \\ x^{14}-\frac{7 x^{11}}{6}+\frac{349 x^8}{720}-\frac{67 x^5}{810}+\frac{1797097 x^2}{389188800} \\ x^{15}-\frac{5 x^{12}}{4}+\frac{4 x^9}{7}-\frac{19 x^6}{168}+\frac{56909 x^3}{6726720}-\frac{1}{9856} \\ x^{16}-\frac{4 x^{13}}{3}+\frac{419 x^{10}}{630}-\frac{1699 x^7}{11340}+\frac{7688249 x^4}{544864320}-\frac{30671 x}{89812800} \\ x^{17}-\frac{17 x^{14}}{12}+\frac{3859 x^{11}}{5040}-\frac{1003 x^8}{5184}+\frac{240005881 x^5}{10897286400}-\frac{620143 x^2}{747242496} \\ x^{18}-\frac{3 x^{15}}{2}+\frac{489 x^{12}}{560}-\frac{137 x^9}{560}+\frac{1466711 x^6}{44844800}-\frac{30439 x^3}{17937920}+\frac{2169}{167552000} \\ x^{19}-\frac{19 x^{16}}{12}+\frac{2489 x^{13}}{2520}-\frac{27569 x^{10}}{90720}+\frac{509000671 x^7}{10897286400}-\frac{406671269 x^4}{130767436800}+\frac{56881351 x}{1099308672000} \\ x^{20}-\frac{5 x^{17}}{3}+\frac{559 x^{14}}{504}-\frac{3373 x^{11}}{9072}+\frac{20123489 x^8}{311351040}-\frac{17248789 x^5}{3269185920}+\frac{84600899 x^2}{596767564800} \\ \end{array} $$

There might be a simpler recurrence relation to generate those polynomials. This is left as an exercise for the reader. :-)