Numerical evaluation of zeta function $\zeta(-1)=1+2+3+\cdots= -1/12$

Question

I'm learning special functions by programming them, largely for fun. This post discussed why or how $\zeta(-1)=-1/12$, while my interest is how to numerically evaluate the zeta function at -1 and other points on $\mathbb{R}$.

Corresponding to the definition of Riemann zeta function $$ \zeta(z)= \sum_{n=1}^\infty n^{-z} $$ where $n^{a+bi}=n^a[\cos(b \ln n)+i\sin (b \ln n)]$, my codes read:

public static double[] R_powerZ(double n, double[] v){
    double a=v[0];
    double b=v[1];
    double theta= b* Math.log(n);
    double[] c={ Math.cos(theta), Math.sin(theta)};
    return M.scale(c,  Math.pow(n, a));
}

public static double[] zeta(double[] tau, int interation) {
    double[] re = new double[2];
    for (int i = 1; i <= interation; i++) {
        double[] v = R_powerZ(i, M.scale( tau ,-1));
        M._add(re, v);
    }
    return re;
}   

Let $z=-1$ and $n=50$, the numerical result is (1275,0). I can see that the value will increase when n increases, does it means that the series is divergent when $z=-1$?. I also tried $\zeta(-0.99999)=(1274.9563...,0)$ which indicates that the function is continuous around -1.

Such numerical evaluation is far from the result of $\zeta(-1)=1+2+3+\cdots=-1/12$. Further, it's crazy to get a negative value from this summation!

I read this old post. Partially convinced by the proof there, my question is: is there a constructive way (by programming, by literal summation as the definition the zeta function) to calculate the value of $\zeta(z)$? If not, why we need a definition of $\sum_{n=1}^\infty n^{-z}$ at all?

Put it differently, how to get a good approximation of the zeta function by programming from scratch? This post did not explain any experiments (e.g. physical experiments or numeric algorithms) from which the value of the zeta function can be measured.

A convenient formula from wiki page is using Bernoulli number $$ \zeta(2n)= \frac{(-1)^{n+1} B_{2n}(2\pi)^{2n}}{2 (2n)!} $$ for positive even integer, and for nonpositive integer $$ \zeta(-n)= (-1)^n\frac{B_{n+1}}{n+1} $$ both give the "right" results of $\zeta(0)$ and $\zeta(-1)$. But how can the two formula comply with the original definition $\sum_{n=1}^\infty n^{-z}$?

Then I tried a 'globally' convergent series $$ \zeta(z)=\frac{1}{1-2^{1-s}} \sum_{n=0}^\infty \frac{1}{2^{n+1}} \sum_{k=0}^n \binom{n}{k} \frac{-1^k}{(k+1)^s} $$ from the wiki page, the corresponding codes yield this image (the color denotes the argument)

Hopefully I am on the right direction.

Answer 1

Numerical evaluation of special functions is a big subject — simply plugging in numbers into the simplest things often gives poor results.

Trying to estimate $\zeta(-1)$ by the partial sums of this series should be expected to give especially poor results — the limit of its partial sums only exists when $\mathrm{Re}(z) > 1$, and $z= -1$ does not have that property.

Similarly, estimating $-1/12$ by the partial sums of the series $1 + 2 + 3 + \ldots$ is also doomed to failure; the summation used here is not "limit of partial sums", it is "zeta regularization", so you should not expect partial sums to have anything to do with the value of the sum.

Wikipedia gives a great many formulas and identities that $\zeta(z)$ satisfies. All of the infinite sums on the wikipedia page, I believe, are the ordinary "limit of partial sums" meaning, so you can use partial sums to estimate the value.