Evaluación numérica de la función zeta de Hurwitz

Question

¿Existe una manera de evaluar numéricamente la función zeta de Hurwitz

$$\zeta(s,a) = \sum_{n=0}^\infty \frac{1}{(n+a)^{s}}$$

que sea más eficaz (es decir, rápida y precisa) que la simple suma explícita de los términos uno por uno?

Answer 1

Tenga en cuenta que los nombres de sus variables no son comunes, voy a utilizar los nombres estándar de DLMF $$\zeta(s,a) = \sum\limits_{n=0}^{\infty}\frac{1}{(n+a)^s}$$ Un buen algoritmo numérico utiliza el Suma de Euler-Maclaurin ( DLMF o Wolfram ) en la forma $$\zeta(s,a) = \sum\limits_{k=0}^{n}\frac{1}{(a+k)^s} + \frac{(a+n)^{1-s}}{s-1} + \sum\limits_{k=0}^{\infty} \frac{B_{k+1}}{(k+1)!}\: \frac{(s+k-1)_k}{(a+n)^{s+k}}$$ donde $(\cdot)_k$ es el símbolo de Pochhammer. Utilizando las propiedades de los números de Bernoulli $B_1=-1/2\;$ y $B_{2k+1}=0\;$ esto se puede reescribir como $$\zeta(s,a) = \sum\limits_{k=0}^{n}\frac{1}{(a+k)^s} + \frac{(a+n)^{1-s}}{s-1} - \frac{1}{2(a+n)^{s}} + \sum\limits_{k=1}^{\infty} \frac{B_{2k}}{(2k)!}\: \frac{(s+2k-2)_{2k-1}}{(a+n)^{s+2k-1}} \cdot $$ El límite superior de la primera suma debe seleccionarse en función de su aritmética (por ejemplo $n=8$ o $n=9$ para el doble, etc.). La segunda suma (infinita) puede calcularse sumando los términos hasta que no haya ningún cambio (es decir, el término actual es de menor magnitud que la suma actual).

Este algoritmo es eficiente para el rango primario $s>0, a>0:\;$ Para la aritmética de punto flotante de 80 bits obtengo $\zeta(2.5,0.75)\approx 2.49154238551193522\;$ con $10$ términos de la segunda suma, y $\zeta(10,0.25)\approx 1048576.10768311475\;$ sólo necesita dos términos.

Para los negativos $s$ puede utilizar Polinomios de Bernoulli (para $s=-1,-2\dots$ ) o el fórmula de reflexión de DLMF .

---

Editar : El número de términos depende principalmente de $s$ y $n$ . Aquí algunos valores intermedios de mi implementación para un doble IEEE de 53 bits. Para los valores $(s,a)\;$ pares $n=9.\;$ s1 es la primera suma de $k=0\dots n,\;$ s2 = s1 + los dos siguientes términos simples, y s2 la suma total después de terminar la última suma.

ZetaH(10, 0.25)
s1 = 1048576.10768311471
s2 = 1048576.10768311494
s3 = 1048576.10768311494       1 term

ZetaH(5.0, 0.25)
s1 = 1024.34894710160870
s2 = 1024.34897386673833
s3 = 1024.34897452658106       6 terms

ZetaH(2.5, 0.75)
s1 = 2.47125718157242114
s2 = 2.49147059743736321
s3 = 2.49154238551193474       8 terms

ZetaH(0.75, 0.25)
s1 =  6.13122575577889695
s2 = -0.938864666031481665
s3 = -0.937591964187151961     9 terms

En las implementaciones reales, se puede ajustar el $n$ y, por supuesto, hay otros métodos de aceleración: Terminar la primera suma, si no hay cambios, precalcular un cierto número de los complicados coeficientes de la segunda suma, etc.

Encontrará una implementación básica en C en el archivo zeta.c de la conocida biblioteca Cephes http://www.moshier.net/double.zip

---

Edición 2 : Aquí los datos de las dos sumas. La primera columna muestra n (el límite superior de la primera suma), j es el índice de la segunda suma, para la que se alcanza el criterio de convergencia. La segunda La segunda suma también termina, si los términos son no decrecientes, lo que es el caso de n=4,5 y los tres últimos pares s,a, en la tabla se muestra como ---.

        s,a=     s,a=     s,a=       s,a=
       10,0.25  5,0.25  2.5,0.75   0.75,0.25
n=4      j=7      ---      ---       ---
n=5      j=3      ---      ---       ---
n=6      j=1      j=8      j=11      j=12
n=7      j=1      j=6      j=9       j=10
n=8      j=0      j=5      j=8       j=8
n=9      j=0      j=5      j=7       j=8
n=10              j=4      j=7       j=7
n=11              j=4      j=6       j=7
n=12              j=4      j=6       j=6
n=13              j=3      j=6       j=6
n=14              j=3      j=5       j=6
n=15              j=3      j=5       j=6
n=16              j=3      j=5       j=5
n=17              j=3      j=5       j=5
n=18              j=3      j=5       j=5
n=19              j=2      j=4       j=5
n=20              j=2      j=4       j=5
n=21              j=2      j=4       j=5
n=22              j=2      j=4       j=5
n=23              j=2      j=4       j=5
n=24              j=2      j=4       j=5
n=25              j=2      j=4       j=4
n=26              j=2      j=4       j=4

Answer 2

La primera fórmula dada por @gammatester concuerda con la dada en el sitio de Wolfram, pero la segunda forma, reescrita, parece estar en error.

El primer y simple error es que el exponente de la $(a+n)$ el término en el denominador de la serie de Bernoulli está equivocado por 2 - debería ser $s+2k-1$ pero un error más complicado es la interpretación del símbolo de Pochhammer al sumar sobre los índices pares, $2k$ . Tal y como se ha dado, en la fórmula reescrita equivale a $(s)2k+1$ . Esto es incorrecto: considere primero la forma Pochhammer en la forma original como $(s + k-1)k$ . Esto significa que el primer factor del símbolo de Pochhammer avanza en uno cada vez, y el último en 2. Por lo tanto, en la suma par debe ser $(s + 2k-2)2k-1$ . El primer término es $(s)1 =s$ el segundo es $(s+2)(s+3)(s+4)$ etc. El último factor es $s+4k-4$ .

Numéricamente la forma más rápida de evaluar en el bucle es incrementar $k$ por uno, dos veces. $(s+k-1)k$ se avanza un paso dividiendo por el primer factor y añadiendo dos más.

He implementado el algoritmo correcto en Fortran y Mathematica (con una precisión superior a la doble), y estoy de acuerdo con @allard en que el término de error después de la primera suma explícita -es decir, con la serie que incluye números Bernoulli- no disminuye como se esperaba. El error se satura mucho antes de que los términos de la serie empiecen a divergir, y no es un problema de precisión.