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

Question

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.

Thanks in advance

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];

    //read parameters (even-parity)
    cout < "psi0 = ";
    cin > psi0;
    cout < "dpsi0 = ";
    cin > dpsi0; 

    dx = x_end / number_steps;

    value_x [0] = 0;
    psi [0] = psi0;
    V [0] = 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;
    }

Answer 1

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

Answer 2

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
    Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eigensolver(H);
    std::cout<< eigensolver.eigenvectors().transpose() << std::endl;
}