6
$\begingroup$

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?

$\endgroup$
7
$\begingroup$

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. :-)

$\endgroup$
10
$\begingroup$

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

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.