% ================================================================ % John Kerl % kerl at math dot arizona dot edu % 2008-02-08 % ================================================================ % ---------------------------------------------------------------- % 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.2; dt = 0.01; h = 0.01; c = 1.0; niter = 10000; x = [xlo:dx:xhi]; V = 0.0 * x; %% ---------------------------------------------------------------- mn = size(x); n = mn(2); H = disclap(n) + diag(V); [eigenvectors,D] = eig(H); %psi0 = exp(-x.^2); psi0 = eigenvectors(:,1); % First eigenvector %psi0 = eigenvectors(:,2); % Second eigenvector %psi0 = 0.4*eigenvectors(:,1) + 0.6*eigenvectors(:,2); psi0 = psi0 / sqrt(sum(conj(psi0).*psi0)); % normalize %% ---------------------------------------------------------------- plot(x, abs(psi0),'b') psie = psi0; psir = psi0; t = 0; for iter = 1:niter % Euler: % psi = psi + i * dt * (del2(psi,dx) - V.*psi) % t = t + dt psie = psie + i * dt * (c * del2(psie,dx) - V'.*psie); % RK4 SE: % k1 = i * (c * del2(psi,dx) - V * psi) % k2 = i * (c * del2(psi + k1*h/2,dx) - V * (psi + k1*h/2)) % k3 = i * (c * del2(psi + k2*h/2,dx) - V * (psi + k2*h/2)) % k4 = i * (c * del2(psi + k3*h,dx) - V * (psi + k3*h)) % psi += dt/6(k1 + 2 k2 + 2 k3 + k4) k1 = i * (c * del2(psir,dx) - V' .* psir); k2 = i * (c * del2(psir + k1*h/2,dx) - V' .* (psir + k1*h/2)); k3 = i * (c * del2(psir + k2*h/2,dx) - V' .* (psir + k2*h/2)); k4 = i * (c * del2(psir + k3*h,dx) - V' .* (psir + k3*h)); psir = psir + dt/6*(k1 + 2*k2 + 2*k3 + k4); %disp(abs(psie)); %disp(abs(psir)); %disp(sqrt(sum(conj(psie).*psie))); %disp(sqrt(sum(conj(psir).*psie))); psie = psie / sqrt(sum(conj(psie).*psie)); % normalize psir = psir / sqrt(sum(conj(psir).*psir)); % normalize t = t + dt; end %disp(abs(psie)); %disp(abs(psir)); hold on;plot(x, abs(psie),'r') hold on;plot(x, abs(psir),'g')