why the simulation period is wrong about schrodinger equation in a harmonic potential

Question

Daniel Niu 2023-6-22

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1986709-why-the-simulation-period-is-wrong-about-schrodinger-equation-in-a-harmonic-potential

评论： Daniel Niu 2023-6-25

clear;  % Clear workspace variables
N = 1024;
x = linspace(-64, 64, N);
dx = x(2) - x(1);
psi = gaussian_wavepacket(x, -32.0, 3.0, 0.0);
V = (x / 32.0).^2 / 2;
H = hamiltonian(N, dx, V);
simulate = create_simulator(H, 1.0);  % Create a new simulator function each time the script runs
figure;
fill([x, fliplr(x)], [V/10.0, zeros(1, numel(V))], 'r', 'FaceAlpha', 0.1, 'EdgeColor', 'none')
hold on;
h_plot = plot(x, probability_density(psi));
hold on;
title('Time Evolution of a Gaussian Wavepacket');
xlabel('Position');
ylabel('Probability Density');
axis([-64 64 0 0.15])
h_text = text(min(x), max(probability_density(psi)), sprintf('t = %.2f', 0), 'VerticalAlignment', 'middle');
M=500;
for k = 1:M
    [psi, time] = simulate(psi);
    set(h_plot, 'YData', probability_density(psi));
    set(h_text, 'String', sprintf('t = %.2f', time));
    movieVector(M) = getframe;
    pause(0.1);
end
function H = hamiltonian(N, dx, V)
    e = ones(N, 1);
    A = spdiags([e -2 * e e], -1:1, N, N);
    L = A / (dx^2);
    H = -L + spdiags(V', 0, N, N);
    H = sparse(H);
end
function psi = gaussian_wavepacket(x, x0, sigma0, p0)
    A = (2 * pi * sigma0^2)^(-0.25);
    psi = A * exp(1i * p0 * x - ((x - x0) / (2 * sigma0)).^2);
    psi = psi.';  % Return as a column vector
end
function U = time_evolution_operator(H, dt)
    U = expm(-1i * H * dt); 
    U(abs(U) < 1E-12) = 0; 
    U = sparse(U);
end
function prob_density = probability_density(psi)
    prob_density = abs(psi).^2;  
end
function simulate = create_simulator(H, dt)
    U = time_evolution_operator(H, dt);
    t = 0;
    simulate = @simulate_func;  % Return handle to nested function
    function [psi, time] = simulate_func(psi)
        time = t * dt;
        psi = U * psi;
        t = t + 1;
    end
end
 %my potential as V = (1/2) * (x/32)^2. I create a harmonic potential of the form (1/2) * m * ω^2 * x^2, 
 %and assuming m = 1 for simplicity, then  ω^2 = 1/32^2, and  ω = 1/32. 
 %This means my classical period T should be T = 2π * 32 = ~200.6
 %but in the simulation, the period is about 142.
 %I don't know why?
 %Your help would be highly appreciated.

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Florian Rössing 2023-6-22

I have studied physics, so I have seen the problem already. But: This is a Matlab Forum. Could you please provide the Math that yields your expected result so we can compare that against the code.

请先登录，再进行评论。

请先登录，再回答此问题。