Perturbation Theory in Mathematica


Hello,
I've wanted to implement perturbation theory in mathematica for some time now. Here, I've done so for a non-degenerate time independent case. Specifically for the one dimensional harmonic oscillator (not quantum because symbolic constants just make the computation take longer). I tested it on some problems from Shankar's quantum mechanics book and it works.
The things you'll need are:
1. Eigenvalues and eigenfunctions for a solved case (in our case the harmonic oscillator).
$$H0=\frac{\text{$\#$1} x^2}{2}-\frac{\nabla _{\{x\}}^{}\text{$\#$1}}{2}\&$$
\[ScriptCapitalH]0 = -(1/2) Laplacian[#1,{x}] + 1/2 (x^2) #1 &;
\[ScriptCapitalH]1 = x^2 #1 &;
\[ScriptCapitalH] = \[ScriptCapitalH]0[#1] + \[ScriptCapitalH]1[#1] &;
{eigList, wfnList} = DEigensystem[\[ScriptCapitalH]0[u[x]],u[x],{x}\[Element] FullRegion[1], nmax, Method -> "Normalize"];
where nmax is the maximum number of eigenvalues/eigenfunctions to compute. The above code stores the eigenvalues in eigList and the eigenfunctions in wfnList.
2. A perturbation
$$H1=\text{$\#$1} x^4\&;$$
Now let's get cooking! To make things easy we'll just use the notation found here Let's define a function, V, that computes the energy of an unperturbed state with the perturbation's hamiltonian: V[n_, m_] := FullSimplify[Integrate[ExpandAll[FullSimplify[Conjugate[wfnList[[n]]] \[ScriptCapitalH1[wfnList[[m]]]]], {x, -\[Infinity], \[Infinity]}]];

And now a function Ei that computes the difference in energy between two unperturbed states.
Ei[n_, m_] := \[ScriptCapitalE]0[n] - \[ScriptCapitalE]0[m];

Now the energy corrections up to fourth order are given by:

\[ScriptCapitalE]0[n_] := FullSimplify[Integrate[ExpandAll[FullSimplify[Conjugate[wfnList[[n]]]\[ScriptCapitalH]0[wfnList[[n]]]]], {x, -\[Infinity], \[Infinity]}]]

\[ScriptCapitalE]1[n_] := Total[Table[V[n, k2], {k2, DeleteCases[Range[nmax], n]}]]

\[ScriptCapitalE]2[n_] := Total[Table[Abs[V[n, k2]]^2/Ei[n, k2], {k2, DeleteCases[Range[nmax], n]}]]

\[ScriptCapitalE]3[n_] := Total[Flatten[ Table[((V[n, k3]*V[k3, k2]*V[k2, n])/(Ei[n, k2]* Ei[n, k3])) - ((V[n, n] (Abs[V[n, k3]])^2)/(Ei[n, k3]^2)), {k3, DeleteCases[DeleteCases[Range[nmax], n], k2]}, {k2, DeleteCases[Range[nmax], n]}]]]

\[ScriptCapitalE]4[n_] := Total[Flatten[ Table[(V[n, k4]*V[k4, k3]*V[k3, k2]*V[k2, n])/(Ei[n, k2]*Ei[n, k3]* Ei[n, k4]) - ((Abs[V[n, k4]]^2/(Ei[n, k4]^2)) (Abs[ V[n, k2]]^2/(Ei[n, k2]))) - (V[n, n] ((V[n, k4]*V[k4, k3]* V[k3, n])/((Ei[n, k3]^2) Ei[n, k4]))) - (V[n, n] ((V[n, k4]*V[k4, k2]* V[k2, n])/((Ei[n, k4]^2) Ei[n, k2]))) + ((V[n, n]^2) ((V[n, k4]^2)/((Ei[n, k4]^3)))), {k4, DeleteCases[DeleteCases[DeleteCases[Range[nmax], n], k2], k4]}, {k3, DeleteCases[DeleteCases[Range[nmax], n], k2]}, {k2, DeleteCases[Range[nmax], n]}]]]

\[ScriptCapitalE][ n_] := \[ScriptCapitalE]0[n] + \[ScriptCapitalE]1[ n] + \[ScriptCapitalE]2[n] + \[ScriptCapitalE]3[ n] + \[ScriptCapitalE]4[n]

Our ground stateto fourth order is then 0.828125. Pretty close.

In the interest of full disclosure I'm kind of over perturbation theory at this point because I figured out how to implement creation and annihilation operators in Mathematica. That may or may not be next. Tune in for more next time.

Thanks,
Alex Krotz

Comments

Post a Comment

Popular posts from this blog

Using Airplane Transponder Data to Characterize Antenna Performance

Image
Hi All, Recently I obtained some low-noise amplifiers in order to better receive NOAA GOES satellites. A low-noise amplifier however is no good without a decent antenna. Most often people use these 2.4 GHz wifi antennas modified to operate at the GOES downlink frequency. This typically gives very good results but it's pretty challenging to DIY a satellite dish. An alternative to satellite dish type antennas is the highly directional yagi-uda antenna. Yagi-uda antennas are pretty (seemingly) simple to build so one night I made a first pass at throwing one together. Alas the satellite signal still did not appear (probably due to my poor assembly). Instead of scrapping the inadequate antenna I decided to try and figure out if adding the passive elements actually did anything to the directionality of the yagi-uda antenna versus a regular dipole antenna. To do this without some sort of signal source is typically impossible and so people resort to modelling programs that calculate the ga
Read more

Microstate table for a d3 electronic configuration.

This was just a test to see if I could put LaTeX formatted tables into this webpage MathJax example $$\begin{array}{|c | c | c | c | c|} \hline \begin{array}{c}\end{array} & \frac{3}{2} & \frac{1}{2} & -\frac{1}{2} & -\frac{3}{2} \\ \hline 5 & \begin{array}{c}\end{array} & ^{2}H\begin{array}{c}^{2}H(2^+,2^-,1^+)\end{array} & \begin{array}{c}(1^-,2^-,2^+)\end{array} & \begin{array}{c}\end{array} \\ \hline 4 & \begin{array}{c}\end{array} & ^{2}G\begin{array}{c}^{2}H(2^+,2^-,0^+)\\^{2}G(2^+,1^+,1^-)\end{array} & \begin{array}{c}(0^-,2^-,2^+)\\(1^-,2^-,1^+)\end{array} & \begin{array}{c}\end{array} \\ \hline 3 & ^{4}F\begin{array}{c}^{4}F(2^+,1^+,0^+)\end{array} & ^{2}F\begin{array}{c}^{2}H(2^+,2^-,-1^+)\\^{2}G(2^+,1^+,0^-)\\^{4}F(2^+,1^-,0^+)\\^{2}F(2^-,1^+,0^+)\end{array} & \begin{array}{c}(1^-,2^-,0^+)\\(0^-,2^-,1^+)\\(0^-,1^-,2^+)\\(-1^-,2^-,2^+)\end{array} & \begin{array}{c}(0^-
Read more