FreeFem++ users | Main >> QuantumHarmonicOscillator

Quantum Harmonic Oscillator

We consider the time-independent Schrodinger equation with a quadratic Potential:

\begin{aligned}- \frac{\hbar^2}{2 m} \Psi''(x) + \frac{1}{2} m \omega^2 x^2 \Psi(x) = E \Psi(x)\end{aligned}

This equation can be simplified to the form $$- \psi''(p) + x^2 \psi(q) = e \psi(q)$$ by the use of the dimensionless quantities \begin{aligned} \psi & = \sqrt{m \omega /\hbar}\Psi \\ p &= \sqrt{\hbar/m \omega }x \\ e &= 2 E/\hbar \omega \end{aligned} The eigensolution of this problem are of the form

$\displaystyle{ e_n = 2 n + 1, \quad \psi_n(q) = \text{const} \cdot \exp \left(- \frac{q^2}{2} \right) H_n(q), \quad H_n(q) = e^{q^2/2} (q - \partial_q)^n e^{-q^2/2} \text{ (Hermite polynomials)}, \quad n = 0,1,2,\dots }$

The following program code calculates the first 6 solutions for the quantum harmonic oscillator (download example: QMoscillator.edp):

int n = 100;
mesh Th = square(n,1,[10.*(2.*x-1.),1.*y]);	
fespace Vh(Th,P1,periodic=[[1,x],[3,x]]); // the periodic boundary condition enforces a one-dimensional space 	
Vh psi, chi; // unknown and testfunction

int nev = 6; // number of eigensolutions
real e = 1.; // energy parameter (eigenvalue)
Vh[int] EigenPsi(nev); // array to store eigenvectors
real[int] EigenVal(nev); // array to store eigenvalues									

varf Schrodinger(psi,chi) = // variational form of the schrodinger equation
	int2d(Th)(dx(psi)*dx(chi) + x^2*psi*chi) 
	- int2d(Th)(e*psi*chi) + on(2,4,psi=0);
varf RHS(psi,chi) = int2d(Th)(psi*chi);
matrix A = Schrodinger(Vh,Vh);
matrix B = RHS(Vh,Vh);

// solve the eigenvalue problem:
int num = EigenValue(A,B,sym=true,sigma=e,value=EigenVal,
vector=EigenPsi,tol=1e-10,maxit=1000,ncv=0);

for(int i = 0; i < nev; i++){
	cout << "Eigenvalue["+i+"] = " << EigenVal[i] << endl;
	plot(EigenPsi[i],dim=3,fill=1,cmm="("+i+")");
	}

Plot for the solution with n = 5:

All eigenfunction drawn in one graph with numerical eigenvalues:

return to 1d examples