- 19

- 0

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

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

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

I obtained [tex]\xi (r)[/tex] 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 [tex]f(k)[/tex] 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

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

and then computing the resulting expression

[tex]\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)[/tex]

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 ([tex]\xi(r)[/tex]) 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.