Let us rewrite that, introducing a timestep time, $\delta t$
$$\psi_{\hbox{after timestep}} \simeq \psi_{\hbox{before}} - i \delta t \left(-\frac{\partial^2 }{\partial x^2}\psi + ax^2\psi\right).$$ You can write $$\frac{\partial^2 }{\partial x^2}\psi(x) \simeq \frac{ \psi(x+1) - 2 \psi(x) + \psi(x-1) }{\delta x^2}.$$ These two approximations give us a simple timestepping function. This is not the only (or best) timestepping function that solves the problem, but perhaps it's one of the simplest. You, the reader, should try to understand the function below. Hopefully, it will make sense with the two equations above.// This returns an empty wavefunction of length n. function wavefunction(n){ var psi = []; for(var i = 0; i < n; i++){ psi[i] = {real: 0, imag: 0} } return psi; // for example [{real:0, imag:0}, {real:0, imag:0}, ...] } // This takes the starting wavefunction, and returns the wavefunction a short time later. function timestep(psi) { // This is how many units of time we're going to try to step forward. var dt = 0.002; var n = psi.length; // This is the wavefunction we're going to return eventually. var ret = wavefunction(n); // We miss off the first and last elements because it looks at the element to the // left and right of this point. for(var i = 1; i < n-1; i++) { // This is the x that is in the equation above. var x = (i-n/2) // This is the potential at this point. var V = x*x * 0.0015; // a here = 0.0015 // We start from the original wavefunction (and later add (dt * dpsi/dt) to it). ret[i].real = psi[i].real; ret[i].imag = psi[i].imag; // This is the kinetic energy applied to psi. var KPsi = { real: psi[i].real * 2 - psi[i-1].real - psi[i+1].real, imag: psi[i].imag * 2 - psi[i-1].imag - psi[i+1].imag }; // This is the potential, applied to psi var VPsi = { real: psi[i].real * V, imag: psi[i].imag * V }; // This is the whole right hand side of the schrodinger equation. var rhsReal = KPsi.real + VPsi.real; var rhsImag = KPsi.imag + VPsi.imag; // This adds it, multiplied by dt, and multiplied by i. // The multiplication by i is what swaps the real and imaginary parts. ret[i].real += rhsImag * dt; ret[i].imag -= rhsReal * dt; } return ret; } // This returns the initial state of the simulation. function init(){ var n = 128; var psi = wavefunction(n); for(var i = 0; i < n; i++){ psi[i].real = Math.exp(-(i-20)*(i-20)/(5*5))*0.75; psi[i].imag = 0; } return psi; }The simulation above is the probability density of an electron in an $x^2$ potential. This is similar to if it was in a bowl, and was rolling backwards and forwards in that bowl.
// This returns an empty wavefunction of length n. function wavefunction(n){ var psi = []; for(var i = 0; i < n; i++){ psi[i] = {real: 0, imag: 0} } return psi; // for example [{real:0, imag:0}, {real:0, imag:0}, ...] } // This takes the starting wavefunction, and returns the wavefunction a short time later. function timestep(psi) { // This is how many units of time we're going to try to step forward. var dt = 0.02; psi = timestepV(psi, dt); psi = timestepT(psi, dt); return psi; } function timestepV(psi, dt) { var n = psi.length; for(var i = 0; i < n; i++) { // This is the x that is in the equation above. var x = (i-n/2) // This is the potential at this point. var V = x*x * 0.0015; var theta = dt * V; var c = Math.cos(theta); var s = Math.sin(theta); var re = psi[i].real * c - psi[i].imag * s; var im = psi[i].imag * c + psi[i].real * s; psi[i].real = re; psi[i].imag = im; } return psi; } function timestepT(psi, dt){ psi = FFT.fft(psi); var n = psi.length; for(var i = 1; i < n/2; i++) { var k = 2 * 3.1415927 * i / n; var theta = k * k * dt; var c = Math.cos(-theta); var s = Math.sin(-theta); var re = psi[i].real * c - psi[i].imag * s; var im = psi[i].imag * c + psi[i].real * s; psi[i].real = re; psi[i].imag = im; var j = n - i; re = psi[j].real * c - psi[j].imag * s; im = psi[j].imag * c + psi[j].real * s; psi[j].real = re; psi[j].imag = im; } return FFT.fft(psi); } // This returns the initial state of the simulation. function init(){ var n = 128; var psi = wavefunction(n); for(var i = 0; i < n; i++){ psi[i].real = Math.exp(-(i-20)*(i-20)/(5*5))*0.75; psi[i].imag = 0; } return psi; }This runs much better. It's faster, mainly because the timestep was longer. But the timestep can be longer because the timestep is more accurate. Likewise, you can get away with a coarser grid, and still get the same accuracy.
BodyWorks: Kinematics ToolA page with a javascript application where you can set body positions, calculate joint angles and animate human motion. |
Physics Algorithms BookThis is a work in progress to write a book on physics algorithms. At the moment, it is about 1/3 finished though, but the latest version can be downloaded for free. |
The double slit and observersA look at the double slit experiment. The first half is meant to be a clear explanation, using simulations. The second half discusses some of the philosophy / interpretations of quantum mecahnics. |
© Hugo2015. Session @sessionNumber