# Computing slightly oscillatory series to high precision?

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:

• Maybe you can use the fact that $cos(kx) = T_k(x)$, where $T_k(x)$ is a Chebyshev polynomial. Then the summation starts to look like a series of rational polynomials. Then if you can turn the series in to a rational polynomial in a Chebyshev basis, it would allow a very efficient way to sum it up. If you aren't familiar with Chebyshev polynomials and basis, the Numerical Recipes in C has a good primer, as well as this: www2.maths.ox.ac.uk/chebfun/ATAP/ATAPfirst6chapters.pdf – Jay Lemmon Jan 2 '14 at 23:56
• er, that should say $cos(kx) = T_k(cos(x))$ – Jay Lemmon Jan 3 '14 at 0:03
• @JayLemmon Thanks you for that link. I'll have a look and see if it helps. – Kirill Jan 5 '14 at 3:51
• I'm joining this party a bit late, but have you tried using Padé approximants, i.e. the $\varepsilon$-Algorithm instead of Richardson extrapolation? – Pedro Feb 4 '14 at 19:57
• By analogy with the case of highly oscillatory integrals, I don't think you'll be able to do a good job without some knowledge of the separation between oscillatory and nonoscillatory parts. If you have such a separation, the Fourier series answer gives you easy exponential convergence. – Geoffrey Irving Feb 5 '14 at 2:31

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.
• I presume you mean $g(x) = \cos(x) / (2 - \cos(x))$? – Geoff Oxberry Feb 4 '14 at 0:51
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: