Recently, @zenterix asked a question related to solving the quantum mechanics of the finite-wall-height particle-in-a-box problem. In the last two decades or so I've been teaching a quantum mechanics course every second year, in which I sketch the method of solution of this problem for the class, I but do not get into the mathematical details. A few times I've started to work on it in Maple, but abandoned it because of lack of time. So I decided to persist this time.
This post is more about some of the challenges and tradeoffs in making it work, hence the title. The title is also a nod to the fact that this problem appears in Engel's textbook [1] in a chapter entitled "Applying the particle in a box model to real-world topics". You don't need to know much about quantum mechanics to understand it. From the mathematical point of view you are solving three ordinary differential equations (in three regions) and stitching the solutions together so that the solution, a wavefunction, is continuous with continuous derivative and tends to zero at plus infinity and minus infinity.
The three regions and the wavefunctions for three eigenvalues (energies) are shown here in the figure, which is the final output of the worksheet. Next is the worksheet and I'll comment further on it below.
(It's a quirk of quantum mechanics that we plot two things with different units (wavefunctions and energy) on the same axis, and worse, we offset one by the other.)
Download finite_box_non-dim2.mw
Comments
1. The first general comment is about the issue of how much Maple can do. In principle one should be able to give dsolve the three ODEs and the boundary conditions and it should output the result. We are still far away from that. Even for one region, the boundary conditions at infinity are not handled by dsolve. The other extreme is to solve the ODEs for their general solutions, and follow the algebraic manipulations for implementing the boundary conditions as one might do it on paper. Engel for example has a comment "At this point we notice that by dividing the equations in each pair, the coefficients can be eliminated to give [...]". Maple will not easily do that trick or the manipulations that led to the special form on which the trick is applied. Levine's textbook [2] gives a different non-intuitive set of steps.
But Maple can solve complicated equations, so I wanted a more general strategy that avoids as much as possible the requirement to manipulate into special forms. On the other hand, it is tempting (as @zenterix tried) to just give the three general solutions with 6 arbitrary constants and the boundary conditions to Maple's solve, and solve for the six constants. One finds that all six are zero. This us a consequence of the fact that the ODEs are linear, and misses the crucial point that it is an eigenvalue problem. So there must be some sort of step-by-step guidance for Maple and there is a general but non-obvious strategy to solving such problems.
2. The second general comment is that there is a trade-off between presenting the solution steps without any arcane Maple code, and getting to the solution efficiently and programmatically, without too many manual manipulations. One of my suggestions to @zenterix involved manually dividing both sides of an equations by a fairly complicated expression, to which @acer expressed a preference for programmatic solutions. At the time, I mainly did that to keep close the to existing solution, and was thinking that I usually set things up so that is not necessary. But I found here that I do it a lot, and it is a natural consequence of the interactive way I use Maple. Here's two examples:
In the second line of label (4), I originally got a more complicated form of de_outside, then manually figured out the factor to put inside simplify to get the form you see. I do this frequently with Maple, especially when I use dchange.
In the last line of label (8), I manually divided by B to remove it. Had B been in the denominator, I could have removed it programmatically with numer (an operation I use frequently, though it somewhat recklessly ignores denominator zeroes).
Programmatically doing things can be relatively unreadable and disrupt the readability for the non-Maple reader. For example in label (6) there is "extraneous" code to find which coefficient to set to zero so that the solution goes to zero at +/- infinity. This was forced on me since after generating the nice figure I re-ran the worksheet only to have multiple errors. I originally used eval(outside, {c__1 = B__L, c__2 = A__L}); to change Maple's dsolve constants to read the way I wanted. But I realized that dsolve does not reproducibly assign c__1 and c__2 to the same term each time, raising the general issue of:
3. Session to session reproducibility. If I am only going to use a worksheet myself, I don't usually worry too much about this, since I can usually debug this fairly easily. On the other hand if the worksheet is for others, it needs to be foolproof. I have taken to always sorting the output of solve, e.g., the assignment to bval in the execution group after label (15). The code added for the dsolve constants works, but is non-trivial Maple (suffixed type testing). (Solving this issue for Eigenvectors is yet more difficult.)
4. Some efforts to make the worksheet more readable were hard to do.
- For some reason, Engel chose to call the two constants in the left region B and B'. Entering B*exp(-k*x) + B'*exp(k*x) in 2D leads to B(x)*exp(-k*x) + diff(B(x), x)*exp(k*x). No doubt this can be fixed, but I didn't pursue this.
- I would normally use upper case Psi for the non-dimensionalized psi, but Psi is the name of a special function in Maple. But wait, now we can use "local Phi;" to solve this problem. However diff(Psi(x),x,x) displays as Psi(2,x), so I had to settle for Phi. (SCR submitted.)
- I build up the left wavefunctions incrementally as expressions assigned to left1, left2, etc, which is rather ugly. Why not use arrow procedures so we can use psi(x), D(psi)(x) etc. Aside from some awkwardness in updating such procedures by redefining them as needed, I ran into the following problem:
limit(A*exp(-k*x) + B*exp(k*x), x = infinity) assuming k > 0; gives as expected
signum(B)*infinity;
but
psi := x -> A*exp(-k*x) + B*exp(k*x);
limit(psi(x), x = infinity) assuming k > 0;
returns unevaluated. What happened to full evaluation?
- I'm used to entering hbar as hbar in 1D; it displays as h with the bar through. In 2-D entering hbar^2 is displayed nicely, but internally it is `ℏ`. We can have the following puzzle:
- I originally wanted a vertical axis label Psi, E/V__0, but labels = [X, Psi, E/V__0] obviously won't work. It took some time to figure out the answer is to make "Psi, E/V__0" an atomic variable. I finally settled on the more honest, if cryptic, label shown.
5. solve gives an "invalid" solution. Originally I worked with solutions left2 etc (see label (6)) after simplify(eval(left2, eqns)). Every second wavefunction was about 10^9 times higher than the others, which was difficult to debug. It turns out the denominator of those wavefunctions was exactly zero but numerically 10^(-10), meaning they plotted as large-scaled versions of the right shape, without any division-by-zero error. The fix is to remove the denominators by scaling (see label (14)).
This is just a consequence of Maple giving generic solutions, and the only way around it is to be aware of it.
Summary
To solve such problems with Maple requires one to play around in an interactive way. After hacking to a solution, it can take some effort to make it presentable and usable to others. But the final result is pleasing, and is certainly much easier than working through the problem on paper. The ability to manipulate complicated expressions means circumventing tricks that might be needed in manual solutions.
References
[1] T. Engel, Quantum Chemistry and spectroscopy. Pearson 2019, 4th ed. Sec. 5.1, Prob. 5.3.
[2] I. N. Levine, Quantum Chemistry. Pearson 2009, 6th ed. Sec. 2.4.
