# 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?

-

## marked as duplicate by Jens, Öskå, rasher, m_goldberg, bobthechemistJun 13 at 22:53

You might look here: link –  eldo Jun 13 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 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 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 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 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][[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.

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