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

3 visualizzazioni (ultimi 30 giorni)
Daniel Niu il 22 Giu 2023
Commentato: Daniel Niu il 25 Giu 2023
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 CommentoMostra NessunoNascondi Nessuno
Florian Rössing il 22 Giu 2023
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.

Accedi per commentare.

### Risposta accettata

David Goodmanson il 23 Giu 2023
Modificato: David Goodmanson il 23 Giu 2023
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.
##### 1 CommentoMostra NessunoNascondi Nessuno
Daniel Niu il 25 Giu 2023
You have solved a big problem for me. Thank you!

Accedi per commentare.

### Categorie

Scopri di più su Signal Processing in Help Center e File Exchange

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!