2
$\begingroup$

I have this code:

V[x_] = If[x <= -1, ∞, 1]*If[-1 < x < -a, 0, 1]*If[-a <= x <= a, V1, 1]*If[a < x < 1, 0, 1]*If[1 <= x, ∞, 1] /. {a -> 1/4}
sg = -psi''[x]/2 + V[x]*psi[x] == EE *psi[x] ;
sol = ParametricNDSolve[{sg, psi[-1]==0, psi'[-1]==1},psi, {x, -1, 1}, {EE, V1}]
Evaluate[psi[5, 30][1] /. sol]
Plot[Evaluate[psi[5, 30][x] /. sol], {x, -1, 1}]

I can't not explain the error:

ParametricNDSolve::ndnum: Encountered non-numerical value for a derivative at x == -1.. >>

and

Infinity::indet: Indeterminate expression 0. \[Infinity] encountered. >>

$\endgroup$
3
  • $\begingroup$ I don't think you can use Infinity in your potential function V. You should instead specify the boundary conditions at x = -1 and x == 1 that the function should be zero. $\endgroup$
    – march
    Dec 8, 2015 at 21:18
  • $\begingroup$ I tried to change V and set the boundary conditions: sol = ParametricNDSolve[{sg, psi[-1] == 0, psi[1] == 0, psi'[-1] == 1}, psi, {x, -1, 1}, {EE, V1}] But this can't be solved: "NDSolve is not currently able to solve boundary value problems with \ discrete variables. " $\endgroup$ Dec 8, 2015 at 21:32
  • $\begingroup$ Thank you for the accept! Usually it's a good idea to wait a day or so before accepting an answer for the purposes of attracting other answers (an accepted answer gets fewer looks). When you have enough reputation, please upvote the answers that you get that you like! Welcome to mathematica.SE. (I changed the title a bit; if you really feel like it, feel free to change it back, but it's a good idea for the titles to be as specific as possible). $\endgroup$
    – march
    Dec 9, 2015 at 0:14

1 Answer 1

3
$\begingroup$

Original Post

The problem is that you have introduced a discontinuity in the derivative of your function by specifying that the potential energy for x <= -1 is infinite. If you remove the factors of If[x <= -1, ∞, 1] and If[x >= 1, ∞, 1], your code runs fine, because it now has no problem applying the boundary conditions at x == -1.

Here is your code, a little cleaned up. You can use one or the other of the definitions of v below:

v[x_, a_?Positive] := Piecewise[{{v1, -a <= x <= a}}]
v[x_, a_?Positive] := v1 UnitBox[x/(2 a)]
Plot[v[x, 1/4] /. v1 -> 10, {x, -1, 1}, Exclusions -> False]

enter image description here

Then, the Schrodinger equation is given by

sg[a_?Positive] := -psi''[x]/2 + v[x, a]*psi[x] == EE*psi[x];

and we solve this via

sol = ParametricNDSolve[{sg[1/4], psi[-1] == 0, psi'[-1] == 1}, psi, {x, -1, 1}, {EE, v1}];

Finally, we use FindRoot to see where the function goes to zero one left boundary and plot. For instance,

pot = 10;
energy = FindRoot[Evaluate[psi[e, pot][1] == 0 /. sol], {e, 5}]
(* { e -> 5.64844} *)

and

Plot[Evaluate[psi[e /. energy, pot][x] /. sol], {x, -1, 1}]

enter image description here

Update

Let's see how well ParametricNDSolve does with this problem. The energy eigenvalues for the v1 -> 1 potential (i.e. just an infinite square well of width 2 units) should be

((π^2)/8)*Range[9]^2

enter image description here

My favorite way of quickly extracting numerical roots is as follows. Plot the function whose roots give us the energies with the roots found, i.e.

p1 = Plot[psi[e, 0][1] /. sol
  , {e, 0, 100}, PlotRange -> {-0.5, 1}
  , MeshFunctions -> Function[#2], Mesh -> {{0}}, MeshStyle -> {PointSize[0.015], Red}]

enter image description here

Extract the zeroes and refine them with FindRoot (this second step isn't entirely necessary, but it gets us better precision) via

rawEnergies = Sort@Cases[Normal@p1, Point[{a_, _}] :> a, Infinity]
refinedEnergies = e /. FindRoot[Evaluate[psi[e, 0][1] == 0 /. sol], {e, #}] & /@ rawEnergies

enter image description here

Finally, how close are they to the real values? Multiply by 8/π^2:

8/π^2 # & /@ refinedEnergies
(* {1., 4., 9., 16., 25., 36., 49., 64., 81.} *)

Excellent!

Finally:

Plot[Evaluate[psi[#, 0][x] & /@ refinedEnergies /. sol], {x, -1, 1}]

enter image description here

Those look like very nice sinusoidal functions!

Energies as a function of the height of the bump

Let's automate the procedure!

energies[v1_, maxE_] := Module[{p1, rawEnergies, refinedEnergies}
  , p1 = Plot[psi[e, v1][1] /. sol, {e, 0, maxE}, PlotRange -> {-0.5, 1}, MeshFunctions -> Function[#2], Mesh -> {{0}}]
  ; rawEnergies = Sort@Cases[Normal@p1, Point[{a_, _}] :> a, Infinity]
  ; refinedEnergies = {v1, e /. FindRoot[Evaluate[psi[e, v1][1] == 0 /. sol], {e, #}]} & /@rawEnergies
 ]

This is somewhat inefficient, but it works well enough. We do this root-finding for a number v1s via

egys = Table[energies[v1, 100], {v1, 0, 30, 0.5}];

This spits out some errors because sometimes FindRoot has some problems, but here it doesn't matter. Finally, we plot this using

Show[
  ListLinePlot[Flatten[egys, {2}]
   , PlotRange -> {0, 40}
   , Frame -> True, FrameLabel -> {"V", Style["E", Italic]}
   , BaseStyle -> {FontFamily -> "CMU Serif", FontSize -> 15}
  ]
  , Plot[Evaluate[\[Pi]^2/8 Range[9]^2], {x, 0, 30}, PlotStyle -> Dashed]
 ]

resulting in

enter image description here

The dashed curves are the particle-in-a-box energies. (By the way, the list of energies that we make might be ragged, because the energies get pushed up by the potential bump. For this reason, we have to transpose a ragged array, and we do this using the magic that is the second argument to Flatten.)

$\endgroup$
2
  • $\begingroup$ Thanks a lot! It works just fine. $\endgroup$ Dec 8, 2015 at 23:20
  • $\begingroup$ But wait, there's more! $\endgroup$
    – march
    Dec 8, 2015 at 23:21

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.