A workspace for
* Trying out ipython (which is very cool, by the way).
* Explaining and coding the topics that Alex and I have been discussing.
Jim Mahoney | cs.marlboro.edu | Oct 2013 | MIT License
second order numerical differential equation
There are several approaches (including particularly the 4th order Runge-Kutta) but here I'll just do the simplest thing. If \(y\) is the dependent variable and \(x\) is the independent, put \(x\) on a grid with spacing \(dx = x[n] - x[n-1]\) where \(n\) is a discrete integer index labeling the grid points.
Then the first derivative is approximately
\[ y'(n + 1/2) = (y(n+1) - y(n))/dx \]
and likewise the second is (to \(O(dx^2)\))
\[ y''(n) = (y(n + 1/2) - y(n - 1/2))/dx \].
Plugging in gives
\[ y''(x_n) = \frac{y(x_{n+1}) - 2 y(x_n) + y(x_{n-1})}{dx^2} \]
or
\[ y_{n+1} = 2 y_n - y_{n-1} + dx^2 y''_n \]
which is implemented in diffeq2
harmonic oscillator
This 2nd order differential equation
\[ \frac{d^2 y}{dx^2} = - y \]
has exact solutions
\[ y = A sin(x) + B cos(x) \]
'Nuff said.
exact pendulum
So let's look at a problem that doesn't have an exact solution.
The equation of the exact pendulum is
\[ m \frac{d^2 L \theta}{dt^2} = - m g \sin(\theta) \]
or
\[ \frac{d^2 \theta}{dt^2} = - \frac{g}{L} \sin(\theta) \]
where \(m\) the mass, \(L\) and \(g\) are constants, while the angle \(\theta\) and time \(t\) are the dynamical variables. (The force is negative because it is in the direction of decreasing \(\theta\).)
Defining dimensionless variables \(x\) and \(y\) in terms of a time scale \(T\) gives
\[ t = x T \] \[ y = \theta \] \[ \frac{d^2 y}{dx^2} = - (T^2 \frac{g}{L}) \sin(y) \]
Choosing \[ T = \sqrt{\frac{g}{L}} \] leaves just
\[ \frac{d^2 y}{dx^2} = - \sin(y) \]
which is nearly the harmonic oscillator but with a nonlinear \(sin(y)\) term. For large amplitudes, the solution won't look like \(sin(x)\).
quantum one dimensional stationary states via "shooting"
The well-known 1D time independent Schrodinger equation is
\[ - \frac{\hbar^2}{2 m} \frac{d^2 \psi_n(x)}{dx^2} + V(x) \psi_n(x) = E_n \psi_n(x) \]
where \(V(x)\) is a binding potential energy, \(\psi_n\) and \(E_n\) are the \(n^{th}\) bound states.
Let's first do the quantum harmonic oscillator, which has the potential
\[ V(x) = \frac{1}{2} k x^2 \]
The first step before we can work this numerically is to change it to a dimensionless form. Define a distance scale \(X\), and energy scale \(\varepsilon\), and dimensionless variables \(s, C_n, y_n\) :
\[ s \equiv \frac{x}{X} \] \[ C_n \equiv \frac{E_j}{\varepsilon} \] \[ y_n \equiv \Psi_n \sqrt{X} \]
(The units on \(\Psi\) can be a bit confusing. Just remember that \(\Psi^2 dx\) is a probability, which is dimensionless.)
Defining two energy constants
\[ Kinetic = \frac{\hbar^2}{2 m X^2} \] \[ Potential = \frac{1}{2} k X^2 \]
lets us write the dimensionless equation as
\[ \frac{d^2 y_n(s)}{ds^2} = \left( \frac{Potential}{Kinetic} s^2 - \frac{\varepsilon}{Kinetic} C_n \right) y_n(s) \]
Since we have two scale factors \(X\) and \(\varepsilon\) which we can set, both of these energy ratios can be set to 1. In other words, the following two equations define the constants \(X\) and \(\varepsilon\) in terms of \(m, \hbar, k\) :
\[ \varepsilon = \frac{\hbar^2}{2 m X^2} \] \[ \frac{1}{2} k X^2 = \frac{\hbar^2}{2 m X^2} \]
and we want to solve this for \(y(s)\) :
\[ \frac{d^2 y_n(s)}{ds^2} = \left( s^2 - C_n \right) y_n(s) \]
for the quantum numbers \(n = (0,1,2,3,...)\).
Given \(N\) points \(x_n\) where \(0 \le n \le N-1\), the complex DFT is defined as
\[ X_k \equiv \sum_{n=0}^{n=N-1} x_n \omega^{- k n} \]
where \(\omega = e^{(0+1j) 2\pi / N}\) is the \(n^{th}\) complex root of 1. (I'm using Pythons notation for complex numbers here, with \(\sqrt{-1}=(0+1j)\).
The inverse transform is
\[ x_n = \frac{1}{n} \sum_{k=0}^{k=N-1} x_n \omega^{ k n} \]
TODO :
Talk about linear algebra, quantum bras and kets, and all that?
Discuss FFT?
The fit looks good by eye, even with a fairly large step. (That may well be because this diffeq Leaving the last one out, here's how the exact and approximate solutions differ.
If we didn't have an exact solution, the typical approach would be to look at how consistent the result is as dx gets smaller. Setting dx = dx/2 let's us compare every second point of the smaller grid with all the points of the bigger grid.
Not as accurate as before, but looks pretty good for the smaller \(dx\). (This runs pretty quickly, but smaller \(dx\) values start to take longer.)
quantum shooting
First, lets test the method with a problem where we know the answer.
The quantum harmonic oscillator has evenly spaced energy levels which with the dimensionless units here are at \(C_n = (1, 3, 5, 7, ...)\).
Here's what a wavefunction looks like when we guess \(C_n = 4.0\) . \(y(s)\) starts at \(y=0\), arcs and crosses the axis before \(s=2\), and then heads down towards \(-\infty\).
Here's what we get solving for the 3rd excited state via shooting. The routine keeps going until of C is consistent to 1e-6. The calculated value isn't correct to that tolerance, but it's pretty close. The wavefunction looks good until near the end, when it starts to blow up. (These solutions always do that eventually, since they're of the form \(A exp(k x) + B exp(- k x)\), and even tiny roundoff errors in \(A\) eventually dominate.)
So ... it mostly works.
And so let's try our DFT on this thing.
First thing to notice: it's slow.
Second thing: all the power is in the lowest few components.
Quick quiz:
* Which wavelengths are the non-zero parts?
* Why are only the imaginary parts nonzero,
and why this form?
* Figure out by hand the wavelengths corresponding
to the values of k which are nonzero.
* Using only the handful of nonzero Y[k]
for k near 0 or N, reconstruct an approximation
of the original function.
