5
$\begingroup$

Following my previous question, I want a graph that must look like this (the lower image is taken from the Wikipedia page for the Microcanonocal Ensemble-->> Quantum Mechanical),

enter image description here

The code for solving Schrodinger's equation:

a0 = 0.02; a1 = 0.64; L = 10;
U[x_] := a0*x^4 - a1*x^2 + a1^2/(4 a0)
{vals, funs} = 
  NDEigensystem[{-1/2 Laplacian[u[x], {x}] + U[x]*u[x], 
    DirichletCondition[u[x] == 0, True]}, u[x], {x, -L, L}, 25, 
   Method -> {"PDEDiscretization" -> {"FiniteElement", {"MeshOptions" \
-> {"MaxCellMeasure" -> 0.001}}}}];

How do we plot the graphs showing level density as shown in the above figure?

N.B.

The author used Python to plot the graph shared above. Please see the Source.

$\endgroup$

1 Answer 1

6
$\begingroup$

You can use ColorFunction to plot the wavefunction densities:

numEig = 25;
Show[Plot[U[x], {x, -2 L/3, 2 L/3}, PlotStyle -> Red, 
  AspectRatio -> 1.2, PlotHighlighting -> None], 
 Table[Plot[{vals[[i]]}, {x, -10, 10}, PlotStyle -> Thick, 
   ColorFunction -> Function[{x, y}, Opacity[2 funs[[i]][x]^2, Blue]],
    ColorFunctionScaling -> False, PlotHighlighting -> None], {i, 1, 
   numEig}]]

enter image description here

As presented in the previous answer, use SmoothHistogram for the density of states. To rotate the plot, see #145789.

SmoothHistogram[vals, 0.1]

enter image description here

$\endgroup$
5
  • $\begingroup$ Thank you @Domen for this quick response. the first graph is absolutely beautiful. However, I believe, the second graph is not correct. The adjacent graph in the Wikipedia page is a probability graph with inverted axes, i.e. the x-axis shows the state probability while y shows its position on the double well as a function of energy. $\endgroup$
    – user444
    Feb 29 at 10:53
  • $\begingroup$ Please see this Python code for the above Wikipedia graph : commons.wikimedia.org/wiki/… $\endgroup$
    – user444
    Feb 29 at 11:02
  • 2
    $\begingroup$ @user444, the plot in your question shows just the density of states, $g(E)$, and that is shown on my plot . The plot in the link you've pasted now shows a canonical ensemble, where you weight the states by the Boltzmann factor to get $\rho(E) = g(E) e^{-\beta E}$. So the two plots are different. $\endgroup$
    – Domen
    Feb 29 at 11:05
  • 2
    $\begingroup$ For the canonical ensemble, you could probably use something like β = 1/10; SmoothHistogram[WeightedData[vals, Exp[-β vals]], 0.1]. $\endgroup$
    – Domen
    Feb 29 at 11:33
  • $\begingroup$ Thank you @Domen, for you kind help. $\endgroup$
    – user444
    Feb 29 at 12:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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