Eigenvalue problem and plotting its eigenfunctions [duplicate]

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?

• You might look here: link
– eldo
Jun 13 '14 at 9:48
• Actually, I am not looking for Harmonic Oscillator. I am intersted how can one in general solve Eigenvalue Problem numerically by Mathematica.
– user14782
Jun 13 '14 at 9:52
• What sort of eigen problem are you trying to solve? Is it a partial differential equation or a set of ordinary differential equations? Some code from you would help. Have you looked at Eigensystems and its variants in Help?
– Hugh
Jun 13 '14 at 10:20
• it is an ordinary 1 dimensional second order differential equation. Imagine potential is 1/2x^2 and I want to obtain eigenvalues and plot eigenfunctions. I would like to learn the general procedure of solving such problems in Mathematica.
– user14782
Jun 13 '14 at 10:23
• @user14782 the method will depend on what Schroedinger equation you want to solve. Single particle or many-body? In the continuum or on a lattice? Is the Hamiltonian time-dependent or not? Even if you specify single-particle, in a continuum and time-independent, you can 1) discretise space and solve, 2) solve exactly part of the problem (eg the potential part), expand the rest in terms of that, and trunacate (1 is a special case of this) etc. Basically you need to be more specific, otherwise the answer to your question would be the bulk of modern computational physics.
– acl
Jun 13 '14 at 10:39

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][]

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

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

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 (expression1) + E^(-(x^2/2)) C (expression2) *)

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

Do a series expansion about x = -Infinity.

y[x] /. sol /. C -> 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 -> 0 /. e -> 1/2 + m // Simplify

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

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

deqn /. sol /. C -> 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.

• It's Great! Thanks!
– user14782
Jun 13 '14 at 17:12
• @user14782 - you should click the up- and down-buttons of Stephen Luttrel's answer.
– eldo
Jun 13 '14 at 18:00
• You are right! I was not familiar with that.
– user14782
Jun 13 '14 at 18:29
• @user14782 ideally you only click the up button if you find the answer helpful, of course :)
– acl
Jun 13 '14 at 19:42

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

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

or just use this function

Eigenvalues[m]
• I would like to solve Schrodinger differential equation and plot the eigenfunctions.
– user14782
Jun 13 '14 at 10:15