13

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:

Function's derivative, $f'(\pi x)$

Kirill
  • 11,438
  • 2
  • 27
  • 51
  • 1
    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: http://www2.maths.ox.ac.uk/chebfun/ATAP/ATAPfirst6chapters.pdf – Jay Lemmon Jan 02 '14 at 23:56
  • 1
    er, that should say $cos(kx) = T_k(cos(x))$ – Jay Lemmon Jan 03 '14 at 00:03
  • @JayLemmon Thanks you for that link. I'll have a look and see if it helps. – Kirill Jan 05 '14 at 03: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 04 '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 05 '14 at 02:31

3 Answers3

7

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.

Geoffrey Irving
  • 3,969
  • 18
  • 41
3

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: Values and derivatives for the series

Geoffrey Irving
  • 3,969
  • 18
  • 41
  • Thank you. The problem is that I selected this specific function as a model for another more complicated function that I actually wanted to evaluate, having similar features, but not actually the same. I'm aware of the closed form from this question on MSE. I meant this as a question about summing an infinite series numerically without closed form. – Kirill Feb 03 '14 at 20:04
  • Maybe my other answer is better then? – Geoffrey Irving Feb 04 '14 at 00:29
0

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.

Community
  • 1
horchler
  • 763
  • 1
  • 5
  • 17