确认 ( x + d/dx ) exp( -x^2 / 2 ) = 0

Question

我写了一个小代码来验证( x + d/dx ) exp(-x^2 / 2 ) = 0。这个想法是使用傅里叶级数 exp( 2 * pi j n x / L ) 具有足够大的 L 来表示高斯并在那里执行操作。

Matlab中的算法工作原理如下：

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

所以上述算法的核心是我调用的函数 Hamilton，它包含运算符 x d/dx 的矩阵表示，这就是为什么 norm( TO * inix' ) 应该 return 一些东西接近于零，但不是（？）并且函数Hamilton如下

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

有人看错了吗？

Answer 1

虽然不喜欢 Matlab ，但我有点遗漏了代码中的一些术语，例如导数的因子 2 pi j k 。所以我在这里放了一个 Python 版本的我认为它应该是什么样子的（对 Python 感到抱歉，但我想它很容易转换为 Matlab）：

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

提供：

[ 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]

这基本上是零。

确认 ( x + d/dx ) exp( -x^2 / 2 ) = 0

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

matlab

math

numerical-methods

differential-equations