4
$\begingroup$

I am relatively new to Mathematica and have been trying to use the NDEigensystem command to work with some quantum systems. I am able to get the accurate energy eigenvalues but am having problems with eigenfunctions; more specifically, I am having trouble actually calling values using the interpolating functions.

m2 = 0.5;
ℏ = 1;
w = 0.5;
\[ScriptCapitalO]2 = -ℏ^2/(2 m2) Laplacian[u[x, y], {x, y}] + 
  1/2 m2 w^2 (x^2 + y^2) u[x, y];

{vals, funs} = 
  NDEigensystem[{\[ScriptCapitalO]2, 
    DirichletCondition[u[x, y] == 0, True]}, 
   u[x, y], {x, -10, 10}, {y, -10, 10}, 28,
   Method -> {"PDEDiscretization" -> {"FiniteElement", {"MeshOptions" \
-> {"MaxCellMeasure" -> 0.5}}}}];

As we can see above, I am using a simple 2d harmonic oscillator as my Hamiltonian, and then using the NDEigensystem command I am generating eigenvalues and eigenfunctions. I am able to get the right eigenvalues for my system as we can see below Energy Eigen Values

However; the eigenfunctions aren't usable, I am trying to extract the data from the interpolating functions to no avail. From what I understand the syntax is:

ψ = funs[[1]]
ψ[2,3]

The above code should print the values of the first eigenfunction as {2,3} but it doesn't seem to be working. I'm hoping to eventually integrate these functions to calculate expectation values, I would be very grateful for any help or advice.

| improve this question | |
$\endgroup$
5
$\begingroup$

1. Use ReplaceAll to inject values for x and y:

funs[[1]] /. {x -> 2, y -> 3}

-0.0555609

2. Construct pure functions (funs2) from funs:

funs2 = Function[{x, y}, #] & /@ funs;

funs2[[1]][2, 3]

-0.0555609

3. use u instead of u[x,y] in the second argument of NDEigensystem

ClearAll[vals, funs3]
{vals, funs3} = NDEigensystem[{\[ScriptCapitalO]2, 
    DirichletCondition[u[x, y] == 0, True]}, 
   u, {x, -10, 10}, {y, -10, 10}, 28, 
   Method -> {"PDEDiscretization" -> {"FiniteElement", {"MeshOptions" -> 
    {"MaxCellMeasure" -> 0.5}}}}];

funs3[[1]][2, 3]

-0.0555609

Plot3D[funs3[[1]][x, y], {x, -10, 10}, {y, -10, 10}]

enter image description here

Plot3D[Evaluate@Through[funs3[[{1, 2, 5, 10, 15}]][x, y]], 
 {x, -10, 10}, {y, -10, 10}, 
 PlotLegends -> (Row[{"funs3[[", #, "]][x,y]"}] & /@ {1, 2, 5, 10, 15}), 
 PlotRange -> All]

enter image description here

| improve this answer | |
$\endgroup$
  • $\begingroup$ Thank you so so much for the detailed and prompt response. You're a real lifesaver, this is working perfectly. $\endgroup$ – rahul menon Jul 13 at 4:24
  • $\begingroup$ Hey, I was wondering if you would be able to help me out again. I was able to use this command perfectly to get the functions but am having some trouble integrating them. Something as basic as this command: funs2 = Function[{x, y}, #] & /@ funs; [Psi]1 = funs2[[1]]; [Psi]2 = funs2[[2]]; a = NIntegrate[[Psi]1 [Psi]2, {x,-10,10},{y,-10,10}] isn't giving any results. $\endgroup$ – rahul menon Jul 16 at 4:19
  • $\begingroup$ @rahulmenon, try a = NIntegrate[[Psi]1[x,y] [Psi]2[x,y], {x,-10,10},{y,-10,10}]? $\endgroup$ – kglr Jul 16 at 5:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.