确认(x + d/dx)exp(-x ^ 2/2)= 0 [英] Confirm that ( x + d/dx ) exp( -x^2 / 2 ) = 0

问题描述

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

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.

Matlab中的算法如下:

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

所以上述算法的核心是我要调用的汉密尔顿函数，它包含运算符 xd/dx 的矩阵表示，这就是为什么 norm(TO * inix')应该返回接近零的值，但不会返回(?)并且汉密尔顿函数如下所示

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

有人看到错误吗?

推荐答案

虽然不在Matlab中，但我还是有点想念代码中的一些术语，例如导数的因子 2 pi j k .所以在这里，我放置了一个我认为应该是什么样的Python版本(对不起Python，但是我想它很容易转换为Matlab):

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

提供:

[ 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的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！