# Solving the 1D Particle-in-a-Box using C++

I've just finished learning the physics behind the problem and would like to write a program in C++ than can solve the problem. I'm actually stuck at the start. I've quite a bit of research, the problem is there is not too many examples of code used to solve the problem.

I'm going to solve the problem using finite-difference form.

$$\dfrac{d^2\psi}{dx^2}\approx \dfrac{\psi_{n+1}+\psi_{n-1}-2\psi_{n}}{(\Delta x)^2}$$

Which allows us to rearrange in the form

$$\psi_{n+1}=2\psi_{n} - \psi_{n-1}-2(\Delta x)^2(E-V_n)\psi_n$$

Using the even-parity solution, we have

$$\psi(0)=1 \quad \quad \psi'(0)=0$$

at $n=0 \implies x=0$

Letting $m=1$ and $\hbar = 1$ I know that at the ground state $E_g =\frac{\pi^2}{8}$

I'm struggling to write code that can help me find $\psi_n$ for values of $n$.

Can I assume that $V_0 = 0$ since we want the wave function to be inside the well?

I'm not asking for someone to write the code for me, I'd just like tips on what I need to define and which method I should use.

EDIT: This is what I have so far

#include <iostream>
#include <string>
#include <stdio.h>
#include <unistd.h>
#include <math.h>
#include <stdlib.h>
#include <stdarg.h>
#include <assert.h>

using namespace std;

const double PI = 3.14159265358979323846264338327950;
int main() {

cout < "1D Particle-in-a-Box"\n;

double psi0, dpsi0, N, dx, x_end, E;
int number_steps, nr;

double value_x [number_steps];
double psi [number_steps];
double V [number_steps];

cout < "psi0 = ";
cin > psi0;
cout < "dpsi0 = ";
cin > dpsi0;

dx = x_end / number_steps;

value_x  = 0;
psi  = psi0;
V  = 0.0
//finite-difference form

while (nr < number_steps)
{
value_x [nr+1] = dx*(nr+1);
psi[nr+1] = 2*psi[nr] - psi[nr-1] - (2*(dt*dt))*(E-V[nr])*(psi[nr]);
nr = nr + 1;
}

• What method do you want to use? Check this answer, for example Nov 21 '17 at 21:19
• Shooting method and runge kutta since that is the easiest apparently. Nov 22 '17 at 18:25
• Finite difference is the easiest, in my opinion. Nov 23 '17 at 2:31
• I don't understand your boundary conditions. Is the box bounded by infinite barriers? If so, assuming the box is between $-L$ and $L$ , you must have on the one hand $\psi(-L)=\psi(L)=0$. On the other hand, the depth of the box being arbitrary, you can put $V(x)=0\ ;\ x \in \{-L;L\}$ no matter what. The solution being even-parity tells you that it is symmetrical with regards to the center of the box, not the value at any $x$. This depends on the potential. Plus, in the ground state your solution is a wave. Saying $n=0 \rightarrow x=0$ makes no sense in this problem. Dec 29 '17 at 11:56
• @Patrick Moloney I think you should focus more on the physics of the basic problem before moving to try to program it. As pointed out by G. Clavier, It seems you have some misunderstandings of the underlying problem. I would also recommend looking over this paper as one way of solving generic 1d potential problems:google.com/url?sa=t&source=web&rct=j&url=http://… Feb 26 '18 at 3:20

I'm not too familiar with QM, but I would use FDM. You could treat it as a steady-diffusion equation with a source term. You would have

$$-\frac{h^{2}}{2 m} \frac{\psi_{i-1} - 2 \psi_{i} + \psi_{i+1}}{\Delta x^{2}} + V \psi_{i} = E \psi_{i}$$

You can rearrange this to be

$$-\frac{h^{2}}{2 m} \frac{\psi_{i-1} - 2 \psi_{i} + \psi_{i+1}}{\Delta x^{2}} = (E-V)\psi_{i}$$

You can use TDMA (Tridiagonal Matrix Algorithm) to solve this ODE as the left-hand side is a tridiagonal matrix and the right-hand side is a forcing function. Therefore for the tridiagonal matix you have

$$a_{i} = -\frac{h^{2}}{2m \Delta x^{2}}\\ b_{i} = \frac{h^{2}}{m \Delta x^{2}}\\ c_{i} = -\frac{h^{2}}{2m \Delta x^{2}}\\ d_{i} = (E - V) \psi_{i}$$

For the boundary condition $\psi(0) = 0$ you make the adjust to the matrix (assuming 1 is your starting index)

$$a_{1} = 0\\ b_{1} = 1\\ d_{1} = 0$$

For the Neumann BC $\psi'(0) = 0$ you will need to evaluate $\psi''$ at $i=1$ and then use a first order approximate $\psi'_{\partial D}$ and then plug in your first order approximate to $\psi''$ evaluate at the boundary.

EDIT: I noticed you said you wanted to use RK4, but I'm confused as you would do that for a Time Dependent Schrodinger Equation TDSE as opposed to TISE. If you are solving TDSE then you could use RK4 for $\frac{\partial \psi}{\partial t},$ however since type of PDE is parabolic I would suggest using implicit time stepping, such as Crank-Nicholson. In the transient case you would have

$$-\frac{h^{2}}{2 m} \frac{\psi_{i-1} - 2 \psi_{i} + \psi_{i+1}}{\Delta x^{2}} + V_{i} \psi_{i} = i h \frac{\psi^{n+1}_{i} - \psi^{n}}{\Delta t}$$

This would give you a $O(\Delta x^{2}, \Delta t)$ approximation.

In explicit form this would give you

• The equation is similar to diffusion considering imaginary time Nov 25 '17 at 4:41
• I agree, to me it seems like it is a parabolic PDE with a forcing function. Nov 27 '17 at 19:05

Easiest solution: set up a grid $$x_0, \ldots, x_N$$ (with spacing $$\Delta x$$) large enough that your wavefunction (approximately) vanishes at the boundary. Then set up the (symmetric tridiagonal) Hamilton matrix

$$H_{ij}= \begin{cases} V(x_i) + \frac{1}{\Delta x^2} && i=j\,,\\ -\frac{1}{2\Delta x^2} && |i-j| = 1\,, \\ 0 && \text{otherwise,} \end{cases}$$

where $$V(x)$$ is the potential and $$\left\{-\frac{1}{2\Delta x^2}, \frac{1}{\Delta x^2}, -\frac{1}{2\Delta x^2}\right\}$$ is the usual three-point finite-difference approximation to the kinetic energy operator $$-\frac12 \frac {\partial^2}{\partial x^2}$$.

Diagonalize this matrix (e.g. using Eigen) to obtain $$M = U E U^T$$. The $$i$$-th column of $$U$$ contains the eigenstate of you system, and the eigenvalue $$E_i$$ is the corresponding eigenenergy.

That's as simple as it can get but at the same time very accurate for most one-dimensional problems.

From your question it seems you want to use the shooting method or similar things. Don't do that. There is no reason in applying the shooting method to a linear one-dimensional problem: it's inefficient and not very stable.

Here is some code to try yourself:

#include <iostream>
#include <Eigen/Dense>

int main()
{
//parameters
int N = 200;
double x0 = -10.0;
double xN = 10.0;
auto potential = [] (double x) {return 0.5*x*x;};

//set up the Hamiltonian
Eigen::MatrixXd H = Eigen::MatrixXd::Zero(N,N);
double dx = (xN - x0) / (N-1);

for(int i=0; i<N; ++i)
{
double x= x0 + i*dx;

if(i>0) { H(i-1,i) = -0.5/(dx*dx);  }
H(i,i) = potential(x) + 1.0/(dx*dx);
if(i<(N-1)) { H(i+1,i) = -0.5/(dx*dx);  }
}

//... and diagonalize it