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

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.

[Florian Rössing]

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.

Answer 1

[David Goodmanson]

Hi Daniel,

Nice animation. Looks like in your 'hamiltonian' function, the kinetic energy is missing a factor of 1/2. Changing the line to

L = (1/2)*A / (dx^2) 

gives a period right around 200.

[Daniel Niu]

You have solved a big problem for me. Thank you!