Handbook of Mathematical Functions

For simulations of many physical systems, one might encounter problems that require the use of special functions and mathematical formula. I found the book by Abramowitz and Stegun "Handbook of Mathematical Functions" is very useful. An electronic copy of the tenth printing of this book can be found in the following link:
http://www.math.sfu.ca/~cbm/aands/

Solving Schrodinger Equation and Computing Density Matrix

I have written two papers about the finite difference time domain (FDTD) method for solving the Schroedinger equation and for computing the single particle density matrix. I have put these papers in the arXiv.org.

[1] arXiv:0712.4317
Title: The Finite Difference Time Domain Method for Computing Single-Particle Density Matrix
Authors: I. Wayan Sudiarta, D. J. Wallace Geldart
Comments: 19 pages, 8 figures

[2] arXiv:0712.4201
Title: Solving the Schrodinger Equation for a Charged Particle in a Magnetic Field using the Finite Difference Time Domain Method
Authors: I. Wayan Sudiarta, D. J. Wallace Geldart
Comments: 8 pages, 4 figures

Solving Schrodinger Equation - Linear Algebra with LAPACK

Today I solve the Time-independent Schrodinger equation (TISE) using the finite difference method that I explained in previous blog and the Matrix method. Being able to solve the TISE numerically is important since only small idealized system can be solved analytically. I consider here only one dimensional case for a particle in a potential V(x). The one dimensional TISE is

-(1/2) [d^2 psi(x) /dx^2] + V(x) psi(x) = E psi(x)

Note I use natural unit m = hbar = 1.

The first step we need to do is to transform differential equation into matrix equation. To use finite difference, we evaluate the wavefunction, psi(x) i.e. the solution of the TISE, at N+1 grid points with spacing dx, {psi(i)}. Similarly the potential is evaluated at i*dx i.e. V(i). Note i=0, 1, ..., N. The computational length is L = N*dx.

Using central finite difference,

d^2 psi(x)/dx^2 = (psi(x+dx) - 2psi(x) + psi(x-dx))/dx^2
= ( psi(i+1) - 2psi(i) + psi(i-1) )/dx^2

This transforms the 1D TISE into
-( psi(i+1) - 2psi(i) + psi(i-1) )/(2dx^2) + V(i)psi(i) = E psi(i)

after simple manipulation, we have

-psi(i+1) + (2 + (2dx^2)V(i))psi(i) - psi(i-1) = (2dx^2)E psi(i)

For computational convenience we use the boundary conditions where psi(0)=psi(L)=0;

The above discretized 1D TISE and using boundary conditions, in matrix form, is given by





This is an eigen matrix problem. To solve it I use LAPACK function dstev_() to compute eigenvalues for symmetric matrix. As an example I consider a finite square well potential given by V(x) = 1 for |x-5|>1 and V(x) = 0 for |x-5|<1 and grid spacing dx=0.1.


--- eigen.c ---

#include <stdio.h>
#include <math.h>

#define N 100

int main(void)
{
int i;
double dx = 0.1;
double dxsq = dx*dx;
double hdiag[(N-1)],hoff[(N-2)]; //hold diagonal and off-diagonal
double v[(N-1)];
char jobType = 'V';
long int n = (N-1);
long int info,ldz = N-1;
double work[2*N-4];
double z[(N-1)][(N-1)];

/* initialize V(x) */
for(i=0;i<(N-1);i++){
if(fabs((i+1)*dx - 5)<1){ v[i]= 0;}
else v[i] = 1;
}

/* initialize hoff */
for(i=0;i<(N-2);i++){
hoff[i] = -1;
}
/* initialize hdiag */
for(i=0;i<(N-1);i++){
hdiag[i] = 2.0 + 2.0*dxsq*v[i];
}

dstev_(&jobType,&n,hdiag,hoff,z,&ldz,work,&info);

/* Print results for four eigen states */
printf("info = %ld\n",info);
printf("The Lowest Four Eigen Values \n");
for(i=0;i<4;i++){
printf("%ld %lf\n",i,hdiag[i]/(2.0*dxsq));
}

printf("The Lowest four Eigen Vectors\n");
printf("x groundstate 1st 2nd 3rd \n");
for(i=0;i<(N-1);i++){
printf("%lf %lf %lf %lf %lf\n",i*dx,z[0][i],z[1][i],z[2][i],z[3][i]);
}

return 0;
}



To compile:

gcc eigen.c -lm -llapack -o try