Computing slightly oscillatory series to high precision?

Question

Suppose I have the following interesting function: $$ f(x) = \sum_{k\geq1} \frac{\cos k x}{k^2(2-\cos kx)}. $$ It has some unpleasant properties, like its derivative not being continous at rational multiples of $\pi$. I suspect a closed form does not exist.

I can compute it by computing partial sums and using Richardson extrapolation, but the problem is that it is too slow to compute the function to a good number of decimal digits (100 would be nice, for example).

Is there a method that can handle this function better?

Here's a plot of $f'(\pi x)$ with some artifacts:

Answer 1

If analytic techniques are disallowed but the periodic structure is known, here is one approach. Let $$g(x) = \frac{\cos x}{2-\cos x}$$ be periodic with period $2 \pi$, so that $$g(x) = \sum_j w_j e^{ijx}$$ where $$w_j = \frac{1}{2\pi}\int_0^{2\pi} g(x) e^{-ijx} dx$$ Thus, $$\begin{aligned} f(x) &= \sum_{k \ge 1} \frac{g(kx)}{k^p} \\ &= \sum_{k \ge 1} \frac{1}{k^p} \sum_j w_j e^{ijkx} \\ &= \sum_j w_j \sum_{k \ge 1} \frac{\left(e^{ijx}\right)^k}{k^p} \\ &= \sum_j w_j \operatorname{Li}_p(e^{ijx}) \end{aligned}$$ You can either approximate the integrals $w_j$ directly or compute a bunch of $f(x)$ values and use a DFT. In either case, you can potentially apply Richardson extrapolation to the result. Since in your case $g(x)$ is analytic within a neighborhood of $\mathbb{R}$, the final series converges exponentially even without Richardson.

Answer 2

For $x = 2\pi a/b$ with $a,b$ integer, we have $$\begin{aligned} f(x) &= \sum_{k \ge 1} \frac{\cos {kx}}{k^2 (2 - \cos {kx})} \\ &= \sum_{k=1}^b \frac{\cos {kx}}{2-\cos {kx}} \sum_{n \ge 0} \frac{1}{(k+bn)^2} \\ &= \sum_{k=1}^b \frac{\cos {kx}}{2-\cos {kx}} \frac{\psi^1(k/b)}{b^2} \end{aligned}$$ where $\psi^1(z)$ is the trigamma function (http://en.wikipedia.org/wiki/Polygamma). Here are plots of the function and its derivative with the artifacts removed:

Answer 3

How about the Levin u-transform? In addition to Fortan codes, there are several versions in the GSL: `gsl_sum_levin_u*'. Matlab's MuPAD and Maple use this scheme.