% ================================================================ % John Kerl % kerl at math dot arizona dot edu % 2008-02-08 % ================================================================ % ---------------------------------------------------------------- % Euler: % psi = psi + i * dt * (c Dpsi - V psi) % t = t + dt % ---------------------------------------------------------------- % Fourth-order Runge-Kutta ("RK4") approximation to % dy/dt = f(t,y) % with % y(t_0) = y_0 % is % y_{n+1} = y_n + h/6(k1 + 2 k2 + 2 k3 + k4) % t_{n+1} = t_n + h % where % k1 = f(t_n, y_n) % k2 = f(t_n + h/2, y_n + k1 h/2) % k3 = f(t_n + h/2, y_n + k2 h/2) % k4 = f(t_n + h, y_n + k3 h) % ---------------------------------------------------------------- % SE: % i hbar dpsi/dt = -hbar^2 / 2m Dpsi + V psi % i dpsi/dt = -hbar / 2m Dpsi + V/hbar psi % dpsi/dt = i (hbar / 2m Dpsi - V/hbar psi) % dpsi/dt = f(psi) -- run RK4's at each x. % % y += y + dt/6(k1 + 2 k2 + 2 k3 + k4) % t += t + dt % where % k1 = f(psi) % k2 = f(psi + k1 h/2) % k3 = f(psi + k2 h/2) % k4 = f(psi + k3 h) % i.e. % k1 = i (hbar / 2m D(psi) - V/hbar (psi)) % k2 = i (hbar / 2m D(psi + k1*h/2) - V/hbar (psi + k1*h/2)) % k3 = i (hbar / 2m D(psi + k2*h/2) - V/hbar (psi + k2*h/2)) % k4 = i (hbar / 2m D(psi + k3*h) - V/hbar (psi + k3*h)) % With hbar = 1 and mass --> c: % k1 = i (c D(psi) - V (psi)) % k2 = i (c D(psi + k1*h/2) - V (psi + k1*h/2)) % k3 = i (c D(psi + k2*h/2) - V (psi + k2*h/2)) % k4 = i (c D(psi + k3*h) - V (psi + k3*h)) xlo = -10; xhi = 10; dx = 0.5; ylo = -10; yhi = 10; dy = 0.5; dt = 0.01; h = 0.01; c = 1.0; niter = 1000; make_eps = 0 [x, y] = meshgrid(xlo:dx:xhi, ylo:dy:yhi); psi0 = exp(-x.^2 + -y.^2); V = 0.0 * x * y; psi0 = psi0 / sqrt(sum(sum(conj(psi0).*psi0))); % normalize psie = psi0; psir = psi0; t = 0; for iter = 1:niter % Euler: psie = psie + i * dt * (c * del2(psie) - V.*psie); % RK4: k1 = i * (c * del2(psir) - V .* psir); k2 = i * (c * del2(psir + k1*h/2) - V .* (psir + k1*h/2)); k3 = i * (c * del2(psir + k2*h/2) - V .* (psir + k2*h/2)); k4 = i * (c * del2(psir + k3*h) - V .* (psir + k3*h)); psir = psir + dt/6*(k1 + 2*k2 + 2*k3 + k4); %disp(abs(psie)); %disp(abs(psir)); %disp(sqrt(sum(sum(conj(psie).*psie)))); %disp(sqrt(sum(sum(conj(psir).*psie)))); psie = psie / sqrt(sum(sum(conj(psie).*psie))); % normalize psir = psir / sqrt(sum(sum(conj(psir).*psir))); % normalize t = t + dt; end figure; surf(x,y,V); shading flat; if make_eps set(gcf, 'PaperPositionMode', 'Manual'); set(gcf, 'PaperPosition', [0 0 5.0 2.0]) print -depsc se2dv end figure; surf(x,y,(abs(psi0))); shading flat; if make_eps set(gcf, 'PaperPositionMode', 'Manual'); set(gcf, 'PaperPosition', [0 0 5.0 2.0]) print -depsc se2d0 end figure; surf(x,y,(abs(psie))); shading flat; if make_eps set(gcf, 'PaperPositionMode', 'Manual'); set(gcf, 'PaperPosition', [0 0 5.0 2.0]) print -depsc se2de end figure; surf(x,y,(abs(psir))); shading flat; if make_eps set(gcf, 'PaperPositionMode', 'Manual'); set(gcf, 'PaperPosition', [0 0 5.0 2.0]) print -depsc se2dr end