Eigenvalue problem and plotting its eigenfunctions

Question

How many different ways can one solve an eigenvalue problem and plot its corresponding eigenfunctions in Mathematica? For example for Harmonic Oscillator? Which one is the most accurate one?

Thanks in advance.

Answer 1

Here is how you can solve the simple harmonic oscillator — i.e. quadratic potential — eigenvalue problem using Mathematica. For simplicity, I set all of the constants to unity.

Define the differential (Schrödinger) equation.

deqn = -(1/2) y''[x] + 1/2 x^2 y[x] == e y[x];

Solve the differential equation.

sol = DSolve[deqn, y, x][[1]]

(*
  {y -> Function[{x}, 
    C[2] ParabolicCylinderD[1/2 (-1 - 2 e), I Sqrt[2] x] + 
    C[1] ParabolicCylinderD[1/2 (-1 + 2 e), Sqrt[2] x]]}
*)

This is the general solution, parameterised by the eigenvalue e and two constants of integration C[1] and C[2].

We have to ensure that the solution is square integrable, so it had better behave itself as x goes to +Infinity and as x goes to -Infinity.

Do a series expansion about x = +Infinity.

y[x] /. sol // Series[#, {x, Infinity, 1}] & // Expand

(* E^(x^2/2) C[2] (expression1) + E^(-(x^2/2)) C[1] (expression2) *)

So, rather than sol, you must now use sol /. C[2] -> 0, so that C[2] = 0 suppresses the divergent E^(x^2/2) term.

Do a series expansion about x = -Infinity.

y[x] /. sol /. C[2] -> 0 // Series[#, {x, -Infinity, 1}] & // Expand

(* E^(x^2/2) (expression1a) + E^(-(x^2/2)) (expression2a) *)

Somehow you must arrange things so that expression1a suppresses the divergent E^(x^2/2) term. In this case expression1a contains the factor

suppress = -((I 2^(1/4 (1 - 2 e)) E^(-I e \[Pi]) Sqrt[\[Pi]])/ Gamma[1/2 (1 - 2 e)]);

If you plot suppress versus e you find that it has zeros at half-integer values of e — i.e. e = 1/2 + m where m = 0, 1, 2, ... — which gives you the familiar ladder of simple harmonic oscillator eigenvalues.

Plot[Abs[suppress], {e, 0, 5}]

(* graphics *)

So the (yet to be normalised) solution of the differential equation is finally

y[x] /. sol /. C[2] -> 0 /. e -> 1/2 + m // Simplify

(* C[1] ParabolicCylinderD[m, Sqrt[2] x] *)

Verify that this is indeed a solution of the differential equation.

deqn /. sol /. C[2] -> 0 /. e -> 1/2 + m // FullSimplify[#, Assumptions -> m \[Element] Integers] &

(* True *)

Variants of this approach can be used solve for the eigenfunctions and eigenvalues of other (sufficiently simple) potentials.

Answer 2

You can solve the characteristic polynomial of your matrix (equation(s)):

CharacteristicPolynomial[ m, \[Lambda]]
Solve[0 == %, \[Lambda]]

or just use this function

Eigenvalues[m]