% Usage: [E,psi,x,Vs] = schrodinger_galerkin_sparse(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). This % function differs from schrodinger_galerkin.m in that it uses % sparse-matrix iterative solvers in order to require O(Nm) storage % and O(Nm^2) time, instead of O(N^2) and O(N^3) for dense matrix methods. % % 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_sparse(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 = spdiags([dx'/6 (dxm+dx)'/3 dxm'/6], -1:1, N, N); 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 % Force B to be exactly Hermitian (in case of numerical errors in makegrid): B = (B + B')/2; K = spdiags([-1./dx' (1./dxm+1./dx)' -1./dxm'], -1:1, N, N); K(1,N) = -1/dxm(1); K(N,1) = -1/dx(N); V0 = zeros(N,1); 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; V0(n) = P + M; end V1 = zeros(N,1); for n = 1:N-1 V1(n) = quadl(@(y) V(y + x(n)) .* y .* (dx(n)-y), 0, dx(n), 1e-8) /dx(n)^2; end Vm = spdiags([V1 V0 [0;V1(1:end-1)]], -1:1, N, N); 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); % hack: we'll shift eigenvalues by +40 to make them all positive % (since we know smallest eig is around -35), and then we can % use 'sm' (smallest-magnitude) Arnoldi, which seems better behaved % than 'sa' for some reason I don't understand. A = K + Vm + 40*B; % check: A should be Hermitian if max(max(abs(A - A'))) > 1e-8 error('A is not Hermitian!'); end opts.disp = 0; opts.maxit = 1000; [psis,Es] = eigs(A,B, m, 'sm', opts); E = diag(Es) - 40; [E,Ei] = sort(E); E = E(1:m); psi = psis(:,Ei(1:m)); Vs = V0 ./ diag(B); % potential values on grid