Confirm that ( x + d/dx ) exp( -x^2 / 2 ) = 0

Question

I have written a small code that is supposed to verify that ( x + d/dx ) exp(-x^2 / 2 ) = 0. The idea is to use a Fourier series exp( 2 * pi j n x / L ) with sufficiently large L to represent the Gaussian and perform the operation there.

The algorithm in Matlab works as follows:

function[] = verify

    epsilon = 0.05; % step-size numerical integration

    N = 40; % number of Fourier coefficients

    L = 30; % window length numerical integration Fourier basis

    X = -L / 2:epsilon:L / 2; % grid

    xFourier = zeros(2 * N + 1); %Allocate space for Fourier coefficients of f(x)=x

    inix = zeros(2 * N + 1); % Allocate space for Fourier coefficients of f(x)=exp(-x^2/2)

    % Compute Fourier coefficients of f(x)=x
    for i1=-N:N
        A = X.*exp(-2 * pi * 1i * i1. * X / L ) / sqrt( L );
        xFourier(i1 + ( N + 1 ) ) = trapz( X, A ); 
    end

    % Compute Fourier coefficients of f(x)=exp(-x^2/2)
    for i1 = -N : N
        A = 1 / sqrt(L) * exp(-X.^2 / 2 ). * exp(-2 * pi * 1i * i1. * X / L );
        inix( i1 + N + 1 ) = trapz( X, A ); % These are the Fourier coefficients of the |x|^2*Gaussian part
    end
    TO = Hamilton( N, xFourier, L );
    norm( TO * inix' )
end

So the heart of the above algorithm is the function Hamilton that I am calling, it contains the matrix representation of the operator x d/dx, which is why norm( TO * inix' ) should return something close to zero, but it does not(?) and the function Hamilton is as follows

function [ Hamilton ] = Hamilton( N, xFourier, L)
    Hamilton = zeros( ( 2 * N + 1 ),( 2 * N + 1 ) );
    for i1 = -N : N
        for i2 = -N : N
            if i1 == i2
                Hamilton( 
                    (i1 + ( N + 1 ) ), ( i2 + ( N + 1 ) ) 
                ) = Hamilton(  
                    ( i1 + ( N + 1)),( i2 + ( N + 1 ) ) 
                ) + 1i * 2 * pi / L * i1;
            end
            if abs( i2 - i1 ) <= N 
                Hamilton( 
                    ( i1 + ( N + 1 ) ), ( i2 + ( N + 1 ) ) 
                ) = Hamilton(
                    (i1 + ( N + 1 ) ), ( i2 + ( N + 1 ) ) 
                ) + xFourier( i1 - i2  + ( N + 1 ) );
            end             
        end
    end
end

Does anybody see a mistake?

Answer 1

While not into Matlab , I somewhat miss a few terms in the code, like the factor 2 pi j k for the derivative. So here I put a Python version of what I think it should look like (sorry for the Python, but I guess it translates to Matlab quite easily):

import numpy as np

## non-normalized gaussian with sigma=1
def gauss( x ):
    return np.exp( -x**2 / 2 )

## interval on which the gaussian is evaluated
L = 10
## number of sampling points
N = 21
## sample rate
dl = L / N
## highest frequency detectable
kmax= 1 / ( 2 * dl )

## array of x values
xl = np.linspace( -L/2, L/2, N )
## array of k values
kl = np.linspace( -kmax, kmax, N )

## matrix of exponents
## the Fourier transform is defined via sum f * exp( -2 pi j k x)
## i.e. the 2 pi is in the exponent
## normalization is sqrt(N) where n is the number of sampling points
## this definition makes it forward-backward symmetric
## "outer" also exists in Matlab and basically does the same
exponent = np.outer( -1j * 2 * np.pi * kl, xl ) 
## linear operator for the standard Fourier transformation
A = np.exp( exponent ) / np.sqrt( N )

## nth derivative is given via partial integration as  ( 2 pi j k)^n f(k)
## every row needs to be multiplied by the according k
B = np.array( [ 1j * 2 * np.pi * kk * An for kk, An in zip( kl, A ) ] )

## for the part with the linear term, every column needs to be multiplied
## by the according x or--as here---every row is multiplied element 
## wise with the x-vector
C = np.array( [ xl * An for An in  A ] )

## thats the according linear operator
D = B + C

## the gaussian
yl = gauss( xl )

## the transformation with the linear operator
print(  np.dot( D, yl ).round( decimals=9 ) ) 
## ...results in a zero-vector, as expected

provides:

[ 0.+4.61e-07j  0.-3.75e-07j  0.+1.20e-08j  0.+3.09e-07j -0.-5.53e-07j
  0.+6.95e-07j -0.-7.28e-07j  0.+6.54e-07j -0.-4.91e-07j -0.+2.62e-07j
 -0.+0.00e+00j -0.-2.62e-07j -0.+4.91e-07j -0.-6.54e-07j  0.+7.28e-07j
 -0.-6.95e-07j  0.+5.53e-07j -0.-3.09e-07j  0.-1.20e-08j  0.+3.75e-07j
  0.-4.61e-07j]

This is basically zero.