Message Boards Message Boards


Quantum Mechanics Without Differential Equations

Posted 1 month ago
9 Replies
1 Total Likes

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}}


In[7]:= optenergies = tripot[[1]]

Out[7]= {73.9328, 86.7515, 144.504, 208.488, 297.988}


In[8]:= Table[ListPlot[Table[{x, psi[x]}, {x, 0, 1, 1/100}] /. tripot[[2, n]], 
  Joined -> True], {n, 5}]

enter image description here

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][2] /. schrod) == 0, {e, #}] & /@ optenergies

Out[10]= {{e -> 73.9364}, {e -> 86.7689}, {e -> 144.564}, {e -> 208.699}, {e -> 

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,

enter image description here

9 Replies

Dr. Frank: Neat application of IPOPT! Just curious how the optimization approach compares to solution via collocation method. I surmise the collocation might not (gracefully) handle the eigenfunction orthogonality constraint?

Hello Frank,

looks impressive - cool.

I tried to run your code ( copy = ctrl c and paste = ctrl v , because I could not download your notebook: error message " Ooops, we couldn't ...." and I got lots of error messages.

Table::itform: Argument 5 at position 2 does not have the correct form for an iterator. >>

ReplaceAll::reps: {Join[cons1,{}]} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing. >>

Having a closer look I think the varable i in

Do[cons2 = Table[vars.varres[[j]] == 0, {j, i - 1}];

is not defined, nor do I find


in the code.

Regards, Hans

Hans, I can e-mail you the code I used for ParametricIPOPTMinimize, if you like.

Hello Frank,

I don't know what has happened earlier, but in the meantime I could download and open the notebook. Alas, to no avail, because I have only Mathematica 7 and no access to ParametricIPOPTMinimize.

Beste Grüße Hans


Some years ago I wrote a Mathematica interface to Ipopt, using MathLInk. After that I bugged Wolfram Research to do something like that in Mathematica and they did a better job than I had. Using ParametricIpoptMInimize, I can set up the problem once and then run it multiple times with different starting values for the variables, with no extra overhead. I find that it can be much more powerful than NMinimize.


Hello Frank,

is there any possibility that I can access IPOPTMinimize under Mma Version 7?


I don't know how that could be done.

I didn't try collocation because I felt it would be easier to do optimization with a complicated objective function and a small number of constraints, as opposed to solving a lot of coupled equations from collocation. .

I can approximate a simple harmonic oscillator as follows:

In[25]:= AbsoluteTiming[sho = cmpt[10000 (x - 1/2)^2, 5, 1/1000];]

During evaluation of In[25]:= {{Solve_Succeeded,1}}

During evaluation of In[25]:= {{Solve_Succeeded,1}}

During evaluation of In[25]:= {{Solve_Succeeded,1}}

During evaluation of In[25]:= {{Solve_Succeeded,1}}

During evaluation of In[25]:= {{Solve_Succeeded,1}}

Out[25]= {7.70231, Null}

 energy levels

In[26]:= sho[[1]]

Out[26]= {99.9994, 299.997, 499.992, 699.985, 899.977}

 test spacing

In[27]:= Table[sho[[1, n]]/(n - 1/2), {n, 5}]

Out[27]= {199.999, 199.998, 199.997, 199.996, 199.995}

In[28]:= Table[ListPlot[Table[{x, psi[x]}, {x, 0, 1, 1/100}] /. sho[[2, n]], 
  Joined -> True], {n, 5}]

enter image description here

Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
or Discard

Group Abstract Group Abstract