Numerical integration - Fourier transform or brute force?

[coccoinomane]

Hi everybody!

I kindly request your help in optimizing the numerical integration of the following expression:

\xi (r)=\frac{1}{2\pi ^2}\int_{-\infty}^{\infty}f(k)\cdot \sin(k\cdot r)\cdot dk

f(k) vanishes outside the boundaries k=0 and k=2; I have got k and f(k) as float arrays, so we are talking about a numerical problem.

I obtained \xi (r) by using the trapezoidal integration method. I would have loved to use something more sophisticated such as polynomial or spline integration, but I could not find anything written in C++ to do that.

The problem with the brute force approach is that sin(kr) make the integrand function oscillate wildly at high r's. Therefore, I need to interpolate the f(k) in order to have enough points (~millions) over which to integrate. This is very time consuming and pretty inelegant, as brute force always is :)

I tried a more "intelligent" approach by expanding

\sin (k\cdot r)=\frac{\exp(ikr)-\exp(-ikr)}{2i}

and then computing the resulting expression

\xi (r)=\frac{1}{4\pi^2 i}\cdot\left(\int f(k)\cdot \exp(ikr)\cdot dk - \int f(k)\cdot \exp(-ikr)\cdot dk\right)

making use of the Discrete Fourier Transform algorithm in Octave and Matlab. This way, however, I do not have any idea on how to relate the output array (\xi(r)) with the r's. Furthermore, I compared the results I obtained with the brute force solution and it turned out that the DFT algorithm solution was wrong.

Here are my questions:
1) is it possible to solve the above integral by making use of some Fourier transform algorithm?
2) if this is not the case, can somebody suggest me a numerical integration method more effective than the trapezoidal one?

Thank you very much for reading this long post :)

Cheers

P.S. this is not homework.

 

[coccoinomane]

Hi,

did I omit some important information? Or maybe I posted in the wrong forum?

By the way, I found a routine in the GNU Scientific Library called "QAWF adaptive integration for Fourier integrals" that could be what I am looking for. I only need to get acquainted with GSL, and it is going to take some time.

Any advice is still welcome :)

Ciao,

Guido

 

[Dr.D]

Are you trying to do a numerical integration over doubly infinite limits? How do you propose to pull that off?

 

[coccoinomane]

Are you trying to do a numerical integration over doubly infinite limits? How do you propose to pull that off?

Hi Dr.D, thank you for the answer.

The function f(k) vanishes for k<0 and k>2, so the infinite integration limits should not be a problem.

 

[Dr.D]

Sorry I overlooked that part of the problem statement. All I saw was that infinite integral.

You might want to look at an old book from about 1970 titled Numerical Methods that Work (the author's name slips away from me right now).

 

[uart]

coccoinomane, why don't you just make a linear interpolation of f between sample points and integrate (not numerically but mathematically in closed form) for each of these trapezoids.

That is, between f(k) and f(k+1) use :

f(x) = f_k + \frac{(f_{k+1}-f_k) \, (x-x_k)}{(x_{k+1} - x_k)} = ax + b

Now you can integrate (ax + b) sin(rx) dx easily enough for each of these linear regions.

BTW.
\int x \sin(rx) \, dx = \frac{\sin(rx)}{r^2} - \frac{x \cos(rx)}{r}

 

[coccoinomane]

uart, thank you very much for the idea! This way I can evaluate the integral almost analytically.
I am trying a slightly modified version of your approach by using a Splines interpolation of f(k) instead of a linear one; I will let you know the outcome. By the way, do you have any idea on how to quantify the error this way?

Thank you again, you saved my day :)

Cheers,

Guido

coccoinomane, why don't you just make a linear interpolation of f between sample points and integrate (not numerically but mathematically in closed form) for each of these trapezoids.

That is, between f(k) and f(k+1) use :

f(x) = f_k + \frac{(f_{k+1}-f_k) \, (x-x_k)}{(x_{k+1} - x_k)} = ax + b

Now you can integrate (ax + b) sin(rx) dx easily enough for each of these linear regions.

BTW.
\int x \sin(rx) \, dx = \frac{\sin(rx)}{r^2} - \frac{x \cos(rx)}{r}

 

[coccoinomane]

Hi Everybody!

I am raising this post from the dead just to list a few algorithm I used with success to integrate the above function:

1) Andrew Hamilton's FFTLog algorithm from http://casa.colorado.edu/~ajsh/FFTLog/. It permits you to integrate not only f(x)*sin(x) but any integrand in the form f(x)*J_u(x) where J_u(x) is the u-th order regular Bessel function. The only downsize it's that it is written in Fortran 77.

2) QAWO routine from QUADPACK. A C version is contained in the GNU Scientific Library (http://www.gnu.org/software/gsl/) while a Fortran 77 version can be downloaded from Netlib (http://www.netlib.org/quadpack). It's a rock solid package that has been tested over the years.

3) Numerical Recipes C++ routine called "dftint" (Discrete Fourier Transform INTegration). It can be found on page 692 of Numerical Recipes 3rd edition; alas it is not free software.

If you need any clarification or help with the above algorithms please feel free to send a Personal Message to me or to e-mail me at coccoinomane gmail dot com.

Cheers,

Guido

 

[maxplankj]

There is a method: "Gauss-Hermite Integration" to aproximate integrals of the form

<br /> \int_{-\infty}^{\infty} \exp(-x^{2}) f(x) dx = \sum_{i=1} ^{N} w_{N,i} f(x_{i})<br />

The N grid point x_{i} can be obtained as the zeros of the N -point Hermite polynomial

find the w_{N,i} in tables

Search the book: Applied Numerical Methods Using MATLAB - Yang Cao Chung and Morris

 

[Polyamorph]

This looks like the sine Fourier transform of a reciprocal space total scattering function into real space.

I'm interested in knowing how the results from using a trapezoidal integration compare to that of the other methods you found? The trapezoidal integration method is an approximation but is it a good approximation? And why is using FFT better, is it just more accurate?

cheers.