# Fast way to compute integral of type $\int dx f(x) \cos(n \pi x)$ in SciPy

I have an integral of the form $$I(n) = \int_0^1 dx f(x) \cos(n \pi x) ,$$ where $n$ is an integer. In other words, I calculate the cosine Fourier coefficients of function $f$, which is real and continuous on the interval. I need to calculate this for a large number of large values of $n$. Currently I'm just doing this:

iint = array([ integrate.quad(lambda x: f(x)*cos(n*pi*x), 0, 1, limit=1000)[0] for n in range(0,nlim) ])


where $nlim$ is typically of the order of $\sim 10^3$. I have the feeling that this is not the most optimal way of doing this. How do I make my calculation faster?

• what is f(x) ....? – Nasser Apr 11 '14 at 11:31
• Is there a reason you don't want to use one of the Fourier Transform libraries? Most of them can compute Discrete Cosine Transforms, too. – Bill Barth Apr 11 '14 at 11:53
• $f(x)$ can be one of several well-behaved (smooth, no fast oscillations, real valued, positive) functions. The FFT library in scipy requires me to input the function at discrete space points. If I want to get the output for large values of $n$, I'd need to first make a very dense grid, calculate this function on this grid and then give it to the function that calculates FFT. I doubt this can be much faster than what I'm doing now, but correct me if I'm wrong. – Echows Apr 11 '14 at 13:02
• It depends, do you need this for only certain values of $n$ or for all $n$ in a range $[1,N]$ where $N$ is big? – Bill Barth Apr 11 '14 at 14:10
• i'd build my own quadrature – Aksakal almost surely binary Apr 11 '14 at 19:02

This is part of the (complex-valued) Fourier transform for which there is (provably) no more efficient way than the Fast Fourier Transform (FFT) if you want to compute the integrals for at least a significant fraction of all frequencies $n$ between zero and the largest frequency you care about.

In other words, while your intuition may tell you that evaluating the function at a large number of points is probably in efficient, the intuition is wrong in this case.

• To amplify Wolfgang: At the very least, you will need a large number of function evaluations to capture the oscillatory nature of $cos(n\pi x)$ if $N$ is large ($O(N)$ of them) and an equivalent number of arithmetic operations on the results. The only thing that using quadrature buys you is the ability not to store the $f(x)$ values at the integration points. By switching to FFTs, you get the advantage of getting all the cosine transform values between 1 and $N$. – Bill Barth Apr 12 '14 at 14:23
• Thats a really good way to think about this. – meawoppl Apr 17 '14 at 16:15
• @BillBarth: Your comment about large numbers of function evaluations isn't necessarily true, e.g. if the integrals are evaluated using Filon transforms. In that case, the accuracy actually increases as the frequency increases. – Lysistrata Mar 13 '16 at 12:47
• @Lysistrata, maybe you could link to an example? – Bill Barth Mar 13 '16 at 16:30
• @BillBarth: I think the best place for examples is in some of the publications and reports by Arieh Iserles and his group at the University of Cambridge. damtp.cam.ac.uk/user/ai/Arieh_Iserles/Publications.html In particular, see: Iserles, A., "On the numerical quadrature of highly-oscillating'integrals I: Fourier transforms", DAMTP Tech. Report NA2003/05, University of Cambridge, UK, June 6, 2003, pp. 26. damtp.cam.ac.uk/user/na/NA_papers/NA2003_05.ps.gz – Lysistrata Mar 13 '16 at 18:16

Isn't this the same as the real part of the Fourier transform of f(x) \begin{align} I(n) &= \int^1_0dxf(x)\cos(n\pi x) \\ &=\Re{{\int^1_0dx f(x) e^{-j n \pi x}}} \\ &= \Re\int^1_0dx f(x) (cos(n\pi x) - j\sin(n \pi x)) \\ &= \int^1_0dxf(x)\cos(n\pi x) \end{align}

If $f$ is very nice, then you can approximate $f$ with a sequence of piecewise polynomials and integrate over the resulting intervals exactly. This may be much, much cheaper than using an FFT. This is also true for approximations of $f$ by any functions where you can easily know or determine the exact antiderivative in each interval. If $f$ is a black box, then you'll have more trouble with this method.