Overview
I think the big take-away for me on this project was realizing how much I miss C++.
It was an absolute pleasure to program, and so fast!
Crank Nicolson
For this exercise I used the Crank Nicolson method, which is an unconditionally stable finite difference method for solving partial differential equations.
To get the Crank Nicolson equation you basically take the standard FTCS equation (Forward Time Central Space) both forward and backwards and combine the equations.
$$ \psi(x, t+h) - h\frac{i \hbar}{4 m a^2}[\psi(x+a, t+h) + \psi(x-a, t+h) - 2 \psi(x, t+h)] \\ = \psi(x,t) + h \frac{i \hbar}{4 m a^2}[\psi(x+a, t) + \psi(x-a, t) - 2\psi(x, t)] $$
If you massage the coefficients you can construct a matrix equation similar to:
$$ \bm{A} \psi(t+h) = \bm{B} \psi(t) $$
And, so for each step we just want to solve for \( \psi(t+h) \). But what’s neat is that the matrixes are Tridiagonal.
It’s easy to calculate \( \bm{B} \psi(t) \), and then we have an equation of the form of \( \bm{A} x = v \). And since it’s tridiagonal we can use Thomas’s algorithm (just gaussian elimination with backpropogation). (Source Code) \( O(n) \) instead of \( O(n^3) \).
Results
An electron bouncing in an infinite well. Not that exciting, but look that it’s 5000 grid points, and simulating thousands of steps takes less than a second. Pretty neat.
What went well
- At first I hand rolled my own matrix math, and then switched to the Eigen package, which was wayyyyy faster.
- The SparseMatrix implementation is also amazingly great since the \( A, B \) matrixes would have been 5000x5000
- Hand rolling the Thomas algorithm was surprisingly smooth (although Eigen does the back propogation for me).
What didn’t go well
- Plotting in c++ is terrible, I was about to use ROOT but eventually just decided to skip that mess and use gnuplot. So plotting is done outside the program. But gnuplot makes gifs which is nice.
- I haven’t done much c++ programming on MacOS, but it’s a bummer that the profilers are so hard to setup. Eventually I found out MacOS has a tool call “Instrument” which worked alright.
- There weren’t any leaks, but it was still nice to check
Future
Two ideas:
- I have a hydrogen orbital simulation which samples the analytically solved solutions. Why not let crank-nicolson do the heavy lifting and solve the problem just from a potential? Is that even possible?
- And then, why not take my c++ code, and compile it to web assembly and see how fast it is to invoke from javascript?
I also saw that there’s a way to numerically solve the Schrodinger equation using some spectral method that I’ll probably look at.
Either way I’m excited! These last few simulations have been my favorite.