% Usage: [E,psi,x,Vs] = schrodinger_galerkin(N, rho, V, m) % % Return the m lowest eigenvalues E and eigenfunctions psi of the Schrodinger % operator for a potential function V(x), on [-1,1] with periodic boundaries. % % Uses a linear-element Galerkin scheme with N points whose spacing % is proportional to the function rho(x). (See makegrid.m). % % E(i) is an array of eigenvalues, and the columns psi(:,i) are the % eigenfunctions corresponding to E(i) at the points x(:). Vs(:) % is an array giving the average potential V(x) at the x(:) elements. function [E,psi,x,Vs] = schrodinger_galerkin(N, rho, V, m) x = makegrid(N, rho); % x_1 .. x_N xm = [ x(N)-2, x(1:N-1) ]; % x_0 .. x_{N-1} xp = [x(2:N), x(1)+2]; % x_2 .. x_{N+1} dx = xp - x(1:N); dxm = x(1:N) - xm; B = diag(dxm+dx)/3 + diag(dx(1:N-1),1)/6 + diag(dxm(2:N),-1)/6; B(1,N) = dxm(1)/6; B(N,1) = dx(N)/6; % check: B should be Hermitian if max(max(abs(B - B'))) > 1e-8 error('B is not Hermitian!'); end K = diag(1./dxm + 1./dx) - diag(1./dx(1:N-1),1) - diag(1./dxm(2:N),-1); K(1,N) = -1/dxm(1); K(N,1) = -1/dx(N); Vm = zeros(N,N); for n = 1:N P = quadl(@(y) V(y + xm(n)) .* y.^2, 0, dxm(n), 1e-8) / dxm(n)^2; M = quadl(@(y) V(y + x(n)) .* (dx(n) - y).^2, 0, dx(n), 1e-8) / dx(n)^2; Vm(n, n) = P + M; end for n = 1:N-1 Vm(n,n+1) = quadl(@(y) V(y + x(n)) .* y .* (dx(n)-y), 0, dx(n), 1e-8) / dx(n)^2; Vm(n+1,n) = Vm(n,n+1); end Vm(N,1) = quadl(@(y) V(y + x(N)) .* y .* (dx(N)-y), 0, dx(N), 1e-8) /dx(N)^2; Vm(1,N) = Vm(N,1); A = K + Vm; % check: A should be Hermitian if max(max(abs(A - A'))) > 1e-8 error('A is not Hermitian!'); end [psis,Es] = eig(A,B); E = diag(Es); [E,Ei] = sort(E); E = E(1:m); psi = psis(:,Ei(1:m)); Vs = diag(Vm) ./ diag(B); % potential values on grid