计算科学 - 薛定谔方程的数值解--多井 - 吾爱随笔录

薛定谔方程的数值解--多井

计算科学有限差分特征值本征系统量子力学

2021-11-24 21:51:14

我正在尝试求解一维四次势阱的允许波函数和能量。

为此，我使用了修补方法（https://engineering.dartmouth.edu/microeng/otherweb/henning/papers/jap-sub95.pdf）。我使用 Numerov 方法从 x = 0 到 x = patch 以及从 x = patch 到 x = end 积分。然后我找到函数的零点

ψ_{l e f t} (p a t c h) * \frac{d ψ_{r i g h t} (p a t c h)}{d x} - ψ_{r i g h t} (p a t c h) * \frac{d ψ_{l e f t} (p a t c h)}{d x}

$\psi_{left}(patch) *\frac{\mathrm{d} \psi_{right}(patch) }{\mathrm{d} x} - \psi_{right}(patch) *\frac{\mathrm{d} \psi_{left}(patch) }{\mathrm{d} x}$

其中左表示从 x = 0 到 x = patch 的积分，右表示从 x = end 到 x = patch 的积分。

我能够为小深度的四次井找到正确的波函数。潜力：

以及发现的前四个波函数：

我遇到的问题是，当我增加井深时，我无法找到正确的波函数。对于以下更深的井：

对于这个更深的井，第一个发现的波函数（下图）显然不是第一个允许的能量，因为它有太多的过零：

当我检查找到根的最小化函数时，我发现它在低能量时不会过零。

我怀疑这个问题是在更深的井中，两个潜在的井可以分开考虑，因为波函数重叠很少。然而，我希望我的求解器足够强大，能够找到孤立的波函数，或者至少能够识别何时可以分别考虑两个势阱。我还希望我的求解器能够求解不对称井，我不确定它现在是否可以做到。有没有办法可以调整我的代码来解决这两个问题？

请在附件中找到我的代码：https ://drive.google.com/drive/folders/0B1OuUJ2C7IHAcHRkcm1JLVpLMWM

此外，这是下面复制的 MATLAB 代码。主要功能是第一个代码块，是脚本。

%% Outputs
% E_allowed; % Allowed eigenvalues
% wavefunctions_allowed; % Allowed eigenfunctions

%% Inputs
potential_function; % potential function, energy in eV
z_axis; % the 1D problem
h; % z_step, z(2) - z(1), distance in meters
mass_particle;
E_guess_max; E_guess_min;
E_guess_step;
E_disp;

%% Set up
% Find eigenfunctions and eigenenergies for an arbitrary potential function
if(~any(strcmp(who,'MAX_SHRO_LOOP')))
    MAX_SHRO_LOOP = 200; % Maximum runs for the bisection loop
end

if(~any(strcmp(who,'TWO_PASS_PATCHING')))
    TWO_PASS_PATCHING = 1; % Maximum runs for the bisection loop
end

if(~any(strcmp(who,'E_guess_step')))
    E_guess_step = 1e-4; % MUST BE SET ACCORDINGLY
    %MUST BE SMALLER THAN THE SPACING BETWEEN ADJACENET ENERGY LEVELS
    %MUST BE LARGE SUCH THAT WE INCREMENT E AFTER EIGEN FOUND WE SEE CHANGE
end

if(~any(strcmp(who,'z_patch_vec')))
    % the x patches that can be tried
    z_patch_vec = 3 : numel(z_axis) - 2;
end

if(~any(strcmp(who,'z_patch_index')))
    z_patch_index = round(numel(z_patch_vec) / 2); % two pass z patching
end

if(~any(strcmp(who,'init_small')))
    init_small = 1e-100; % have to adjust to get good results
    % if determinant is too small then increase this number--if the eigenvalues
    % are found too quickly
    % if NaNs in wavefunctions then decrease
end


if(~any(strcmp(who,'num_to_find')))
    % adaptive eigenenergy finding...
    num_to_find = [30 60 100 500 Inf];
end

if(~any(strcmp(who,'accuracy_to_find')))
    accuracy_to_find = [1e-16 1e-8 1e-8 1e-8 1e-8]; 
    % after each num_to_find has been reached, change the accuracy
end

% E in eV
if(~any(strcmp(who,'E_guess_min')))
    E_guess_min = min(potential_function);   
end


if(~any(strcmp(who,'E_guess_max')))
    E_guess_max = max([potential_function(1) potential_function(end)]);    
end

if(~any(strcmp(who,'E_disp')))
    E_disp = 0;
end




%%
mass_particle = mass_particle / q;

Es_guessed = [];
crossfires_guessed = [];

% the calculated stuff...
E_allowed = [];
wavefunctions_allowed = [];
left_allowed_wavefunctions = [];
right_allowed_wavefunctions = [];

E_guess = E_guess_min;
E_bracket_low = E_guess_min;
E_bracket_high = E_guess_min;

% for displaying stuff
E_guess_range = abs(E_guess_max - E_guess_min);
curr_percent = 10;
fprintf('\nShrodinger: ');

z_patches = []; % the x_patches tried
if_loop_counter = 1; % the counter for the adaptive x pacthing

while E_guess < E_guess_max

    % adapative eigenenergy error finding
    EIGEN_ERR = accuracy_to_find(sum(numel(E_allowed) > num_to_find) + 1);  
    % adaptive x patching
    z_patch = z_patch_vec(z_patch_index);
    % print progress bar
    if( 100 * abs(min(potential_function) - E_guess) / E_guess_range >= curr_percent )
        fprintf(1,'%2d%%... ',(curr_percent));
        curr_percent = curr_percent + 10;
    end

    % Find a good x_patch value

    % Low energy bracket
    T = h^2 * 2 * mass_particle * (potential_function - E_bracket_low) / (12 * h_bar^2); % T function
    guess_left_bracket_low  = Numerov_Left(0, init_small, T);
    guess_right_bracket_low = Numerov_Right(0, init_small, T);
    [crossfire_mat_bracket_low val_low] = make_crossfire_mat(guess_left_bracket_low, guess_right_bracket_low, z_patch, h);

    % high energy bracket
    E_bracket_high = E_bracket_low + E_guess_step;
    T = h^2 * 2 * mass_particle * (potential_function - E_bracket_high) / (12 * h_bar^2); % T function
    guess_left_bracket_high  = Numerov_Left(0, init_small, T);
    guess_right_bracket_high = Numerov_Right(0, init_small, T);
    crossfire_mat_bracket_high = make_crossfire_mat(guess_left_bracket_high, guess_right_bracket_high, z_patch, h);
    %--this code can be deleted...

    % poor patching crossfire returns NaN, try another patching value
    if(isnan(crossfire_mat_bracket_low) ||  isnan(crossfire_mat_bracket_high))

        if_loop_counter = 1 + if_loop_counter;
        z_patch_index = mod(z_patch_index + 1, numel(z_patch_vec)) + 1;
        z_patch = z_patch_vec(z_patch_index);

        if(if_loop_counter > numel(z_patch_vec))
            disp('patching'); % all the x patchings were tried, none worked...
            break;
        end

    end
    %--this code can be deleted...

    %--this if statement (but not the code in the block) can be deleted
    % only continue if the x patching is okay and neither crossfire is not NaN
    if(~isnan(crossfire_mat_bracket_low) && ~isnan(crossfire_mat_bracket_high))

        if_loop_counter = 1; % reset the x patch counter so we can do it again later if needed
        z_patches = [z_patches; z_patch]; % the x_patches we have tried

        crossfires_guessed = [crossfires_guessed; val_low];
        Es_guessed = [Es_guessed; E_bracket_low];  

        %-------------------------- bracket guesses --------------------------%
        % Want to bound the eigenenerngy between the low and high brackets
        % Now having found a good x_patch, let us continue to bracket the eigenvalue guess

        % low energy bracket is either the min energy (if 1st try) or the
        % previous eigenenergy + guess_step

        % high energy bracket search
        % determinant has to flip signs
        while(crossfire_mat_bracket_low * crossfire_mat_bracket_high > 0)

            E_bracket_high = E_bracket_high + E_guess_step;

            T = h^2 * 2 * mass_particle * (potential_function - E_bracket_high) / (12 * h_bar^2); % T function

            guess_left_bracket_high  = Numerov_Left(0, init_small, T);
            guess_right_bracket_high = Numerov_Right(0, init_small, T);

            [crossfire_mat_bracket_high val] = make_crossfire_mat(guess_left_bracket_high, guess_right_bracket_high, z_patch, h);
            crossfires_guessed = [crossfires_guessed; val];
            Es_guessed = [Es_guessed; E_bracket_high];  

        end

        counter = 1; % counter for the binary search
        found_flag = 0; % 1 = found an eigenenergy

        % set the first guess to be the low bracket
        E_guess_previous = E_bracket_low;
        crossfire_mat_previous = crossfire_mat_bracket_low;

        % guess
        E_guess = (E_bracket_high + E_bracket_low) / 2;

        % binary search
        % search until MATLAB resolution reached or iteration counter satisfied
        while(abs(E_guess_previous - E_guess) > EIGEN_ERR && counter < MAX_SHRO_LOOP && ~found_flag)

            if(E_guess > E_guess_max)
                found_flag = 1;
                break
            end

            % for the bisection guess
            T = h^2 * 2 * mass_particle * (potential_function - E_guess) / (12 * h_bar^2); % T function
            guess_left  = Numerov_Left(0, init_small, T);
            guess_right = Numerov_Right(0, init_small, T);        
            [crossfire_mat val] = make_crossfire_mat(guess_left, guess_right, z_patch, h);

            crossfires_guessed = [crossfires_guessed; val];
            Es_guessed = [Es_guessed; E_guess];  

            % E guess correct
            if(crossfire_mat == 0)

                E_allowed = [E_allowed; E_guess];
                right_allowed_wavefunctions = [right_allowed_wavefunctions, guess_right];
                left_allowed_wavefunctions = [left_allowed_wavefunctions, guess_left];
                found_flag = 1; % will not add to the allowed stuff at the end of the bisection loop

            % guessed too low, no sign change
            elseif(crossfire_mat * crossfire_mat_previous > 0)

                E_bracket_low = E_guess;

            else % guessed too high and the sign changed

                E_bracket_high = E_guess;

            end

            % diagnostic display
            if(E_disp)
                disp(E_guess);
            end

            % go onto the next bisection loop iteration
            E_guess_previous = E_guess;
            crossfire_mat_previous = crossfire_mat; 
            E_guess = 0.5 * (E_bracket_low + E_bracket_high); % New guess

            counter = counter + 1;            

        end

        % if the binary search stopped
        % stopped because the difference between two eigenenergies was so
        % small or because the max iterations was reached...not because
        % determinant was zero
        if(~found_flag)

            E_allowed = [E_allowed; E_guess];
            right_allowed_wavefunctions = [right_allowed_wavefunctions, guess_right];    
            left_allowed_wavefunctions = [left_allowed_wavefunctions, guess_left];

        end

        % increment the brackets for the next run
        E_bracket_low = E_allowed(end) + E_guess_step;
        E_bracket_high = E_bracket_low + E_guess_step;

    end

end


%%
% Clean up wavefunctions
% wavefunction correction
wavefunctions_allowed = zeros(numel(z_axis), numel(E_allowed));
patching = [];
for i = 1 : numel(E_allowed)

     [wavefunctions_allowed(:,i) patch_out] = wavefunction_patch(left_allowed_wavefunctions(:,i), right_allowed_wavefunctions(:,i), z_axis, z_patches(i), E_allowed(i), potential_function);
     patching = [patching; patch_out];

end

% flag_zerocross_eigenlevel = 1; % so far ok
% % num nodes ~= the energy level
% for i = 1 : numel(E_allowed)
%    
%     % zero crossings not commensurate with eigenvalue level
%     if(numel(find_roots(wavefunctions_allowed(:,i))) ~= i - 1)
%         flag_zerocross_eigenlevel = 0;
%         break
%     end
%     
% end
% 
% if(~flag_zerocross_eigenlevel)
%     
%     %disp('ERROR');
%     % the wells are too seperate, consider them as isolated!!
%     
% end


%% Two pass patching, now we do it again...having found the optimal patching where the agreement is closest. Necessary?
% the calculated stuff...
if(TWO_PASS_PATCHING)

    E_allowed = [];
    wavefunctions_allowed = [];
    left_allowed_wavefunctions = [];
    right_allowed_wavefunctions = [];

    % E in eV
    E_guess_min = min(potential_function);
    E_guess_max = max([potential_function(1) potential_function(end)]);

    E_guess = E_guess_min;
    E_bracket_low = E_guess_min;
    E_bracket_high = E_guess_min;

    % for displaying stuff
    E_guess_range = abs(E_guess_max - E_guess_min);
    curr_percent = 10;
    fprintf('\nShrodinger: ');

    while E_guess < E_guess_max && numel(E_allowed) < numel(patching)

        % adapative eigenenergy error finding
        EIGEN_ERR = accuracy_to_find(sum(numel(E_allowed) > num_to_find) + 1);  

        % second pass patching
        z_patch = patching(numel(E_allowed) + 1);

        % print progress bar
        if( 100 * abs(min(potential_function) - E_guess) / E_guess_range >= curr_percent )
            fprintf(1,'%2d%%... ',(curr_percent));
            curr_percent = curr_percent + 10;
        end

        % Low energy bracket
        T = h^2 * 2 * mass_particle * (potential_function - E_bracket_low) / (12 * h_bar^2); % T function
        guess_left_bracket_low  = Numerov_Left(0, init_small, T);
        guess_right_bracket_low = Numerov_Right(0, init_small, T);
        [crossfire_mat_bracket_low val_low] = make_crossfire_mat(guess_left_bracket_low, guess_right_bracket_low, z_patch, h);

        % high energy bracket
        E_bracket_high = E_bracket_low + E_guess_step;
        T = h^2 * 2 * mass_particle * (potential_function - E_bracket_high) / (12 * h_bar^2); % T function
        guess_left_bracket_high  = Numerov_Left(0, init_small, T);
        guess_right_bracket_high = Numerov_Right(0, init_small, T);
        crossfire_mat_bracket_high = make_crossfire_mat(guess_left_bracket_high, guess_right_bracket_high, z_patch, h);

        crossfires_guessed = [crossfires_guessed; val_low];
        Es_guessed = [Es_guessed; E_bracket_low];  

        %-------------------------- bracket guesses --------------------------%
        % Want to bound the eigenenerngy between the low and high brackets
        % Now having found a good x_patch, let us continue to bracket the eigenvalue guess

        % low energy bracket is either the min energy (if 1st try) or the
        % previous eigenenergy + guess_step

        % high energy bracket search
        % determinant has to flip signs
        while(crossfire_mat_bracket_low * crossfire_mat_bracket_high > 0)

            E_bracket_high = E_bracket_high + E_guess_step;

            T = h^2 * 2 * mass_particle * (potential_function - E_bracket_high) / (12 * h_bar^2); % T function

            guess_left_bracket_high  = Numerov_Left(0, init_small, T);
            guess_right_bracket_high = Numerov_Right(0, init_small, T);

            [crossfire_mat_bracket_high val] = make_crossfire_mat(guess_left_bracket_high, guess_right_bracket_high, z_patch, h);
            crossfires_guessed = [crossfires_guessed; val];
            Es_guessed = [Es_guessed; E_bracket_high];  

        end

        counter = 1; % counter for the binary search
        found_flag = 0; % 1 = found an eigenenergy

        % set the first guess to be the low bracket
        E_guess_previous = E_bracket_low;
        crossfire_mat_previous = crossfire_mat_bracket_low;

        % guess
        E_guess = (E_bracket_high + E_bracket_low) / 2;

        % binary search
        % search until MATLAB resolution reached or iteration counter satisfied
        while(abs(E_guess_previous - E_guess) > EIGEN_ERR && counter < MAX_SHRO_LOOP && ~found_flag)

            if(E_guess > E_guess_max)
                found_flag = 1;
                break
            end

            % for the bisection guess
            T = h^2 * 2 * mass_particle * (potential_function - E_guess) / (12 * h_bar^2); % T function
            guess_left  = Numerov_Left(0, init_small, T);
            guess_right = Numerov_Right(0, init_small, T);        
            [crossfire_mat val] = make_crossfire_mat(guess_left, guess_right, z_patch, h);

            crossfires_guessed = [crossfires_guessed; val];
            Es_guessed = [Es_guessed; E_guess];  

                % E guess correct
            if(crossfire_mat == 0)

                E_allowed = [E_allowed; E_guess];
                right_allowed_wavefunctions = [right_allowed_wavefunctions, guess_right];
                left_allowed_wavefunctions = [left_allowed_wavefunctions, guess_left];
                found_flag = 1; % will not add to the allowed stuff at the end of the bisection loop

                % guessed too low, no sign change
            elseif(crossfire_mat * crossfire_mat_previous > 0)

                E_bracket_low = E_guess;

            else % guessed too high and the sign changed

                E_bracket_high = E_guess;

            end

            % diagnostic display
            if(E_disp)
                disp(E_guess);
            end

            % go onto the next bisection loop iteration
            E_guess_previous = E_guess;
            crossfire_mat_previous = crossfire_mat; 
            E_guess = 0.5 * (E_bracket_low + E_bracket_high); % New guess

            counter = counter + 1;            

        end

        % if the binary search stopped
        % stopped because the difference between two eigenenergies was so
        % small or because the max iterations was reached...not because
        % determinant was zero
        if(~found_flag)

            E_allowed = [E_allowed; E_guess];
            right_allowed_wavefunctions = [right_allowed_wavefunctions, guess_right];    
            left_allowed_wavefunctions = [left_allowed_wavefunctions, guess_left];

        end

        % increment the brackets for the next run
        E_bracket_low = E_allowed(end) + E_guess_step;
        E_bracket_high = E_bracket_low + E_guess_step;

    end

    %%
    % Clean up wavefunctions
    % wavefunction correction
    wavefunctions_allowed = zeros(numel(z_axis), numel(E_allowed));
    patching_second_pass = [];
    for i = 1 : numel(E_allowed)

         [wavefunctions_allowed(:,i) patch_out] = wavefunction_patch(left_allowed_wavefunctions(:,i), right_allowed_wavefunctions(:,i), z_axis, patching(i), E_allowed(i), potential_function);
         patching_second_pass = [patching_second_pass; patch_out];

    end

end

clearvars MAX_SHRO_LOOP TWO_PASS_PATCHING E_guess_step z_patch_vec z_patch_index init_small...
    accuracy_to_find num_to_find E_guess_min E_guess_max

function u = Numerov_Left(init1, init2, T)
    % init1 = x(1) wavefunction value, must be zero for bound particle
    % init2 = x(2) wavefunction value, can be anything b/c normalization
    % T function for Numerov

    u = zeros(numel(T), 1);
        u(1) = init1; % boundary condition
        u(2) = init2; % boundary condition

    for i = 2 : numel(T) - 1
        % Solve for u(x + h)
        u(i + 1) = (( (2 + 10 * T(i)) * u(i) - (1 - T(i - 1)) * u(i - 1) ) / (1 - T(i + 1)));

    end

end

function [sign_val, val] = make_crossfire_mat(left, right, patch_index, h)

    crossfire = [left(patch_index) right(patch_index);...
        (left(patch_index) - left(patch_index - 1)) / h ...
        (right(patch_index) - right(patch_index -1)) / h];
    val = det(crossfire);

    sign_val = sign(val);

end

function u = Numerov_Right(init1, init2, T)
    % init1 = x(1) wavefunction value, must be zero for bound particle
    % init2 = x(2) wavefunction value, can be anything b/c normalization
    % T function for Numerov

    u = zeros(numel(T), 1);
        u(end) = init1; % boundary condition
        u(end - 1) = init2; % boundary condition

    for i = numel(T) - 1 : -1 : 2
        % Solve for u(x - h)
        u(i - 1) = (( (2 + 10 * T(i)) * u(i) - (1 - T(i + 1)) * u(i + 1) ) / (1 - T(i - 1)));

    end

end

function [out patch] = wavefunction_patch(left, right, z_axis, z_patch, eigenenergy, band_function)
% patch the left and right wavefunctions together

    % for symmemtric potentials
    if(sum(band_function - fliplr(flipud(band_function))) == 0)
        patch = round(numel(band_function) / 2);
    else

        % patch at location that is closest and nonzero
        % the division will make it nonzero because zero will go to inf!!
        left = left * right(z_patch) / left(z_patch);
        [val patch] = min((abs(left.^2 - right.^2)) ./ (left.^2));

        % patch at a peak
        temp = find_roots(gradient(left));
        patch = temp(1); % return the patching point so we can go back and rerun the crossfire with this new patching point
    end 

    out = [left(1 : patch); right(patch + 1 : end)];
    out = out / sqrt(trapz(z_axis, out.^2));

end

2个回答

虽然我无法为您的具体实施提供帮助，但我想指出一种替代方法（如对 phil 的回答的评论中已指出的那样）：Marston 的“傅立叶网格哈密顿量”（FGH）。

FGH 是一种伪谱方法，它为构造有界系统的离散哈密顿矩阵提供了一个简单的方法。像往常一样，离散哈密顿量的特征向量和特征值表示离散波函数和相应的能量。对于 FGH，哈密顿量 $H = T+V$ 被分裂成动能部分 $T$ 和潜在的部分 $V$ . 在位置空间中的矩阵元素 $\left< x\vert V\vert x'\right>$ 只是 $V(x)\delta(x- x')$ ，所以势能只贡献哈密顿矩阵的对角线。对于动能算子的矩阵元素，带限（离散 $x_i$ ！）使用平面波基础和矩阵元素 $\left< x_i\vert T\vert x_j\right>$ 是解析计算的。

您会在 David Tannor 的“量子力学简介：与时间相关的观点”一书中找到对此方法和各种相关傅立叶方法的精彩讨论。那里也给出了原始参考资料。

您会在此处找到一个带有简单 Python/Cython 实现的 Jupyter 笔记本，该实现解决了四次势。要更改您使用的电位，请找到与的线pot = lambda x: x**4 - 20*x**2并将其更改为您感兴趣的任何其他电位。确保仅计算约束状态 $\left<\psi\vert\psi\right>$ 在你的边界附近衰减到零 $x$ 领域。

寻找薛定谔方程特征值的标准方法称为“虚时间传播”。您更改坐标 t=-i\tau，并在 \tau 方向上积分。任何随机初始条件都会收敛到最低能量本征态。得到的方程通过分裂方法求解：首先使用快速傅里叶变换传播动能，然后传播势能。这是可能的，因为解快速（指数）衰减为 |x| 增长，您可以通过周期性条件有效地替换边界条件。

这是一些二阶方法的matlab代码

% in order to set up the problem, you need a spatial grid and a momentum grid, 
% since you will be using FFTs, the number of grid points should be a power of 2
Dx  = 0.01; 
x   = [-1:Dx:1-Dx];

% Now you define the eigenvalues in Fourier space
k           = zeros(1,Len).';
k(2:Len/2)  = 1:Len/2-1;                %  (x, 1, ..., Len/2-1, x, x, ..., x)
k((Len/2+1):Len) = -k((Len/2+1):-1:2);  %  (x, ...,x, -Len/2+1, ..., -1)
k(Len/2+1) = 0;

% from this you derive the momentum operator in Fourier space
p = -1i k.^2;

% This allows you to define the full problem, kinetic energy and potential, adjust as you please
E_kin = p.^2 / 2/mass;
V     = x.^4 - x.^2 ;

% you want to solve e^(h(T+V))Psi
% you approximate by the Strang splitting e^(h/2 V)e^(h T)e^(h/2 V)
for j=1:1:n
        Psi = ifft( Psi);
        Psi =  exp(- 1i*h* E_Kin ) .* Psi;
        Psi =  fft( Psi);
    % Apply potential V
        Psi = exp(  -1i*h/2* V)
    % Normalize the result in each step
         Psi = Psi / ( norm(Psi,2)^2*Dx);
   end
end

可以在这里找到更复杂（更高阶）的方法，以及对该问题的许多参考：http: //personales.upv.es/serblaza/2013JChemPhys.pdf

由于您想要更多的本征态，因此有几个过程： 1：迭代。计算第一个本征态。然后，您重新启动该过程，但在每一步中，您都将投影减去此状态。例如，对于第二个状态 psi_2，在每一步中你写 psi_2 = psi_2 - *psi_1。

作为一个完整的替代方案，您可以使用有限差分（或再次使用 FFT）离散拉普拉斯算子，然后使用逆幂方法。给定您的哈密顿量 H，迭代 n： (H-\lambda^(-1)I)^(-1)\psi_(n+1) = \psi_n 对特征值 \lambda 进行一些猜测。您可以在维基百科上找到详细信息。

其它你可能感兴趣的问题

上一篇约束模拟退火下一篇初学者的弹性动力学/波动方程的速度-应力公式