The ground state of a quantum mechanical problem can be found
by minimizing the expectation value of the energy.
The first excited state can be found the same way,
by making it orthogonal to the ground state, etc.
The following function takes three arguments,
the potential energy function, the number of states to calculate
and the spacing between the discretized values of x.
It returns a list of the energy of the states
and a list of a list of the wavefunction values.
It assumes [HBar]^2/(2m) = 1.
The second derivative in kinetic energy is approximated
using a finite difference. The calculation is done between x = 0 and x = 1.
The optimization is performed using ParametricIPOPT,
using a function defined in the attached notebook.
Only one random search (with random seed = 0)is needed for each energy level
cmpt[pefunc_, n_Integer, d_] :=
Block[{i, vars, psi, vwb, ener, enres, varres, cons1, cons2, allcons,
sln, wavefunctions},
vars = Table[psi[x], {x, 0, 1, d}];
vwb = {#, -5, 5} & /@ vars;
ener = d Sum[-psi[x] ((psi[x + d] - 2*psi[x] + psi[x - d])/d^2) +
psi[x] pefunc psi[x], {x, d, 1 - d, d}];
enres = Table[0, n];
varres = Table[0, n];
wavefunctions = Table[0, n];
cons1 = {psi[0] == 0, psi[1] == 0, vars.vars == d^-1};
Do[cons2 = Table[vars.varres[[j]] == 0, {j, i - 1}];
allcons = Join[cons1, cons2];
sln = iMin[ener, allcons, vwb, 1, 0];
enres[[i]] = sln[[1]];
wavefunctions[[i]] = sln[[2]];
varres[[i]] = vars /. sln[[2]],
{i, n}
]; {enres, wavefunctions}
]
calculate energy levels and wavefunctions for first 5 states
in an infinite square well potential with triangle potential inside.
In[6]:= tripot = cmpt[If[x < 1/2, 200*x, 200 (1 - x)], 5, 1/100];
During evaluation of In[6]:= {{Solve_Succeeded,1}}
During evaluation of In[6]:= {{Solve_Succeeded,1}}
During evaluation of In[6]:= {{Solve_Succeeded,1}}
During evaluation of In[6]:= {{Solve_Succeeded,1}}
During evaluation of In[6]:= {{Solve_Succeeded,1}}
energies
In[7]:= optenergies = tripot[[1]]
Out[7]= {73.9328, 86.7515, 144.504, 208.488, 297.988}
wavefunctions
In[8]:= Table[ListPlot[Table[{x, psi[x]}, {x, 0, 1, 1/100}] /. tripot[[2, n]],
Joined -> True], {n, 5}]
use energy levels found by optimization
to solve Schrodinger’s Equation, using ParametricNDSolve.
In[9]:= schrod = ParametricNDSolve[{-psi''[x] +
If[x < 1/2, 200*x, 200 (1 - x)]*psi[x] == e psi[x], psi[0] == 0,
psi'[0] == 1}, psi, {x, 0, 1}, {e}];
The energy levels from Schrodinger’s Equation
are close to those found using optimization.
In[10]:= schrodenergies = FindRoot[(psi[e](https://community.wolfram.com//c/portal/getImageAttachment?filename=opt_result.jpg&userId=29126) /. schrod) == 0, {e, #}] & /@ optenergies
Out[10]= {{e -> 73.9364}, {e -> 86.7689}, {e -> 144.564}, {e -> 208.699}, {e ->
298.49}}
The wavefunctions found using Schrodinger’s Equation
are close to those found by optimization.
Note that these functions are not normalized,
since the initial condition [Psi]'[0] == 1 was used.
Table[Plot[psi[e][x] /. schrod /. schrodenergies[[i]], {x, 0, 1}], {i,
5}]
ParametricIpopt Code.nb (12.5 KB)
