Azzera filtri
Azzera filtri

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

4 visualizzazioni (ultimi 30 giorni)
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 Commento
Florian Rössing
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
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.

Più risposte (0)

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!

Translated by