# How to implement a numerically efficient Airy Zeta Function

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? ## 2 Answers Here is a version that produces exact results and memoizes the derivatives and polynomials for efficiency: airyZetadAA[n_ /; IntegerQ[n] && n >= 0] := airyZetadAA[n] = If[n == 0, AiryAi[airyZetaz]^2, D[airyZetadAA[n - 1], airyZetaz]] airyZetadAB[n_ /; IntegerQ[n] && n >= 0] := airyZetadAB[n] = If[n == 0, AiryAi[airyZetaz] AiryBi[airyZetaz], D[airyZetadAB[n - 1], airyZetaz]] airyZetadBoA[n_ /; IntegerQ[n] && n >= 0] := airyZetadBoA[n] = If[n == 0, AiryBi[airyZetaz]/AiryAi[airyZetaz], D[airyZetadBoA[n - 1], airyZetaz]] airyZetapoly[n_ /; IntegerQ[n] && n >= 2] := airyZetapoly[ n] = π/ Gamma[n] (airyZetadBoA[n - 1]/(3^(2/3) Gamma[1/3]^2) + airyZetadAB[n - 2] - Sum[Binomial[n - 1, j] airyZetadBoA[n - 1 - j] airyZetadAA[ j - 1], {j, 1, n - 1}]) /. airyZetaz -> 0 /. Gamma[2/3] -> (2 π)/(Sqrt[3] Gamma[1/3]) /. Gamma[1/3] -> Sqrt[(2 π)/(3^(1/6) airyZetax)] // Expand airyZeta[n_ /; IntegerQ[n] && n >= 2] := airyZeta[n] = airyZetapoly[n] /. airyZetax -> (2 π)/(3^(1/6) Gamma[1/3]^2) airyZetaN[n_ /; IntegerQ[n] && n >= 2] := airyZetaN[n] = airyZetapoly[n] /. airyZetax -> 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. :-) • You might be interested in this. – J. M. is away Oct 26 '15 at 6:24 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.)

• Your second version is not giving the correct answers. – Mark Adler Oct 26 '15 at 16:23
• @Mark, prolly due to the inadvertently inserted line break. Try again. – J. M. is away Oct 26 '15 at 16:29
• Yep, that fixed it. – Mark Adler Oct 26 '15 at 17:00