Compare analytical and numerical Sine Transform

Question

According to the FFTW Website, the Fourier Sine Transform (FST) returns:

$$Y_k = 2 \sum_{j=0}^{N-1} X_i \sin [\pi (j+1)(k+1)/(N+1)]$$

WolframAlpha defines the Fourier Sine Transform as follows: $2\sqrt\frac{\lvert b \rvert}{(2\pi)^{1-a}} \int_0^\infty f(t)\sin(b\omega t)\mathrm{d}t$

Taking $a=1$ and $b=\pi$ this becomes: $F^{W}_{sin} = 2\sqrt \pi \int_0^\infty f(t)\sin(b\omega t)\mathrm{d}t$.

Comparing the two definitions one can write:

$$Y_k = 2 \sum_{j=0}^{N-1} X_i \sin [\pi (j+1)(k+1)/(N+1)] \approx \frac{1}{\sqrt\pi} F^{W}_{s}$$

Setting $f(t) = t\mathrm{e}^{-t^2}$ and performing FourierSinTransform, Wolframalpha returns:

$$FST\{f(t)\}= \frac{1}{2}\pi^2\omega\mathrm{e}^{-(1/4)\pi^2\omega^2}$$

I implemented this in my code and I was puzzled by the results: The analytical and numerical solutions look similar but I would have expected a higher precision. What is the reason for this? Am I doing some error in reasoning?

Any help is appreciated.

// RosenbluthFourier.cpp : Definiert den Einstiegspunkt für die Konsolenanwendung.
//

#include "stdafx.h"
#include <fftw3.h>



int main() {

int const Pi = 3.14159265359;
int N;
double sup;

cout << "enter N of points "; std::cin >> N;
cout << "enter sup "; std::cin >> sup;
cout << "Interval runs from 0 to " << sup << " and will be divided into " << N
    << " intervals." << endl;

double T, Df;
T = sup / (N-1);
Df = 1 / sup;
cout << "Sampling interval T = " << T << endl;
cout << "Frequency spacing df = " << Df << endl << endl; 

double *X = new double[N];
double *Y = new double[N];

for (int i = 0; i <= N; i++) {
    X[i] = T*i;
    Y[i] = X[i] * exp(-pow(X[i],2));
    cout << "X[" << i << "] = " << X[i] << "  Y[" << i << "] = " << Y[i] << endl;
}

cout << endl << "Analitically tranformed function" << endl << endl ;

double *f = new double[N];
double *Yt = new double[N];

for (int k = 0; k <= N; k++) {
    // calc pi*w
    f[k] = Pi*k*Df;
    Yt[k] = (1./2.)*Pi*f[k] * exp(-pow(f[k]/2., 2));
    cout << "f[" << k << "] = " << f[k] << "  Yt[" << k << "] = " << Yt[k] << endl;
}

cout << endl << "FFTW-tranformed function" << endl << endl;

fftw_plan p;
p = fftw_plan_r2r_1d(N, Y, Yt, FFTW_RODFT00, FFTW_ESTIMATE);
fftw_execute(p);

for (int k = 0; k <= N; k++) {
    Yt[k] = Yt[k] * T * double(sqrt(Pi)) ;
    cout << "f[" << k << "] = " << f[k] << "  Yt[" << k << "] = " << Yt[k] << endl;
}

fftw_destroy_plan(p);

return 0;
}

Answer 1

int const Pi = 3.14159265359;

The above is equivalent to int const Pi = 3; as far as C++ is concerned, which is not a good choice when you care about precision.

Try the following, instead, which also adds more decimals, since the double type has about 16 significant digits (at least if the implementation is IEEE-754 compliant):

double const Pi = 3.14159265358979324;

---

[ EDIT ]     for (int i = 0; i <= N; i++) {

The above wil run the loop $N+1$ times due to the <= inclusive comparison. But the arrays are only allocated for $N$ elements double *X = new double[N]; so the last iteration invokes undefined behavior (out-of-bounds write). Everything that happens from that point on is technically undefined by the language standards. In practice, the code is corrupting some random memory locations which may, or may not, affect the end results.

Answer 2

The FourierSineTransform was performed with 64 points in the interval [0,5]

Analytical and numerical FST