• I assume this would be totally impractical with problems of greater than 1-dimension, though of course for separable problems, you can solve for each separable part individually.
• There's no reason, necessarily, to do this all in the x basis. We could instead solve for the fourier or laguerre series of the eigenfunction, which might converge faster and/or have a nicer form for the derivative.

v ← a random vector
for some number of iterations,
   v ← H v
   v ← v / (v dot v)   # make v a unit vector again
end

A v = λ v             definition of eigenvector/eigenvalue
A-1 A v = λ A-1 v     multiply both sides by A inverse
A-1 v = (1/λ) v       re-arrange

for jj=1:N_eigenvectors_to_find
    % start with a random vector
    v = rand(size(x));
    for ii=1:N_iterations
        % Apply the map v <-- A v
        v = Hinv * v;
        % Subtract the projections of this vector onto
        %  all previously found eigenvectors
        for kk=1:(jj-1)
            v = v - vs(:,kk) * vs(:,kk)' * v;
        end
        % Make this a unit vector by dividing by its length
        v = v/sqrt(v'*v);
    end
    % Store the newly found eigenvector in our list
    vs(:,jj) = v;
end

psi[n_, l_, m_, r_, θ_, φ_] := 
   Factorial[ n + l] / 
     (Factorial[ 2 l + 1] Sqrt[(2 κ)^3 (2 n Factorial[n - l - 1])])
   Exp[-κ r] (2 κ r)^l 
   Hypergeometric1F1[-n + l + 1, 2 l + 2, 2 κ r] 
   SphericalHarmonicY[l, m, θ, φ]

hbar = 1.05457148e-34;  % Planck's constant in m^2 kg / s

while (1),
    % Advance the state of the system    
    psidot = (-i/hbar) * hamiltonian(psi);    
    psi = psi + psidot * dt;    
end

• Functional derivative
• Functional integration
• Path integral formulation
• Wiener measure
• Bernoulli polynomial
• Perron–Frobenius theorem

In lecture today Prof. Rajeev applied the variational principle to the expected energy 〈ψ|H|ψ〉 [where H=p²/2m + V(x)], using the method of Lagrange multipliers to minimize E(ψ)=〈ψ|H|ψ〉 subject to the constraint G(ψ)=〈ψ|ψ〉=1. The result is that the condition Hψ=Eψ pops out, e.g. the lowest energy is the energy of an eigenstate of H. But didn't we already know this about the ground state energy, i.e. that we can decompose any state into eigenstates of H, and one of those eigenstates must have lowest energy?

Any time we have a commutation relation of the form [N,X] = cX for a scalar c, we will have an operator X that acts as a raising (or lowering, if c < 0) operator for the eigenvalues of N.

NX |n〉 = NX |n〉
NX |n〉 = (XN + [N, X]) |n〉
NX |n〉 = (XN + c X) |n〉
NX |n〉 = XN |n> + cX |n〉
NX |n〉 = Xn |n> + cX |n〉
NX |n〉 = (n+c) X|n〉
 X |n〉 = a |n+c〉 for some a

Submitted to Wikipedia as Raising operator.

Okay. Clebsch-Gordon coefficients. Say we have a spin j₁ particle which lives in space V₁ and a spin j₂ particle which lives in space V₂. Each of these particles can have the z-axis projection of its spin range from -j...j, so each particle has 2j+1 possible states, labelled as |j m〉. This means that for the system of two particles, there are (2j₁+1)(2j₂+1) possible states in the "product space". For instance, if we have two spin ½ particles, with the m=-½ state labelled as |↓〉 and the m=+½ state labelled as |↑〉, then a basis of product states for the tensor product space V₁⊗V₂ is {|↑〉,|↓〉}⊗{|↑〉,|↓〉}={|↑〉|↑〉, |↑〉|↓〉, |↓〉|↑〉, |↓〉|↓〉}. In this basis we've labelled states with |m₁ m₂〉 (with j₁ and j₂, which are fixed, implied). What we want to do is form a change of basis from the |m₁ m₂〉 basis to a "total spin" basis |J M〉. We accomplish this principally using two facts: (1) the "extremal states" in the two bases are identical. Namely, |J=(j₁+j₂) M=J〉=|j₁ j₁〉|j₂ j₂〉. (2) On this identity we can apply the lowering operator, J_. On the "total spin" basis, this operator just decrements M by one and spits out some coefficient (which we can ignore). On the product space basis, we must remember that J_ = J_⁽¹⁾ + J_⁽²⁾, i.e. it is the sum of two operators, one which acts on each space. More properly, this is written J_ = J_⁽¹⁾⊗I⁽²⁾ + I⁽¹⁾⊗J_⁽²⁾. Essentially we can fill in another row of the change of basis matrix each time we apply the lowering operator J_. This will get us all of the J=j₁+j₂) states. However, to get a complete total-spin basis, you also get lower values of J, which I haven't quite figured out yet.

Anyways. Here's a sample numerical application of this in Matlab. I want to compute the good-old ½⊗½, i.e. a system of two spin-½ particles. First we get the familiar pauli spin matrices as our Jx, Jy, and Jz:

j = 1/2;
[Jx,Jy,Jz,Jminus,Jplus] = Joperators(j);

Now let's prepare a "maximal" state. Since our basis is ordered from m=-j to m=+j, the maximal state is the last state in the basis:

v = [0 0 0 1]';

Now we can form our total-J lowering operator as J_ = J_⁽¹⁾⊗I⁽²⁾ + I⁽¹⁾⊗J_⁽²⁾. Note that, since j₁=j₂ for this example, the matrix representations of J_⁽¹⁾ and J_⁽²⁾ are the same.

JAYminus = kron(eye(2),Jminus) + kron(Jminus, eye(2));

This is a 4x4 matrix which operates on the product space. Now watch this magic:

>> (JAYminus^0 * v)'
ans =  0     0     0     1
>> (JAYminus^1 * v)'
ans =  0     1     1     0
>> (JAYminus^2 * v)'
ans =  2     0     0     0

Ignore the numerical values for the moment, since I've thrown away some square-roots and hbars in the J operators. What this is telling us is that, assuming that |J=1 M=1〉 is identically equal to |↑↑〉, we can lower that once to find out that |J=1 M=0〉 is made out of equal parts of |↑↓〉 and |↓↑〉, and that |J=1 M=-1〉 is identically |↓↓〉. We just have to normalize each of these states to get our new set of orthonormal eigenstates spanning the total-J space—well, the subspace with J=j₁+j₂. How to systematically get the J<j₁+j₂ states is the next step.

I never feel like I actually understand something until I can write a computer program to compute or simulate it, since writing a computer program necessarily requires a full understanding of the algorithm behind some process. Something that's long seemed rather mysterious to me is Clebsch-Gordon coefficients, which give you the decomposition of a direct product space of angular momenta into a direct sum space of angular momenta. Or something like that. Anyway, as a first start to writing a program to compute Clebsch-Gordon coefficients, I've written this little program that computes the Jx, Jy, and Jz matrices for a z-axis spin eigenstate of total spin j. This program omits a square root from Jx and Jy, and a factor of Planck's constant from everything. This is so that I can work with integers and simple fractions instead of decimal approximations.

function [Jx,Jy,Jz] = Joperators(j)

% The basis ranges from <j -j| to <j j| in increments of 1
mvalues = -j:1:j;
cardinality = length(mvalues);

% Jminus is zero except for directly above the diagonal
% There is an implicit \hbar and square root

Jminus = zeros(cardinality);
for i=2:cardinality,
    m = mvalues(i);
    Jminus(i-1, i) = j*(j+1) - m*(m-1);
end

Jplus = zeros(cardinality);
for i=1:(cardinality-1),
    m = mvalues(i);
    Jplus(i+1, i) = j*(j+1) - m*(m+1);
end

Jx = (1/2) * (Jplus + Jminus);
Jy = (-i/2) * (Jplus - Jminus);
Jz = diag(mvalues);