Sunday, February 07, 2010

Python vs. MATLab

Whenever I need to put together quickly some numerical code, I immediately turn to MATLab. The problem with programs like MATLab and Mathematica is that they are very "packaged", i.e. it is not very easy to integrate those codes in to a regular SHELL session. If I want to further process the calculation result using other programs (e.g. gnuplot to make very pretty LaTeX plot), then I am almost always forced to write the result to some output file using fprintf or similar, then process that file. In principle, if the output is well formatted, one can run the MATLab code like this:

$ matlab -nodisplay -nodesktop -nojvm -r "run_my_code(args);exit;"
However, MATLab would also print out a long header as well as its own command prompts. Additional work is needed to get rid of those undesirables and make use of the print-out.

Python, specifically SciPy, comes to the rescue. Take a look at this piece of code I put together in five minutes, which finds the eigenvalues---a non-trivial numerical feat---of a simple 4x4 hermitian matrix.

#!/usr/bin/python

## USAGE: GBLLandau.py n Vb g1 w2B i, gives the n-th LL of the i-th band

from numpy import sqrt,mat # square root and matrix
from scipy.linalg import eigh # for hermitian eigenfunction
import sys

# obtain parameters from input arguments
argc = len(sys.argv)
if not (argc == 5 or argc == 6):
        sys.stderr.write("USAGE: GBLLandau.py n Vb g1 w_0^2 {i}\n");
        sys.exit(1)

n   = int(sys.argv[1]);
Vb  = float(sys.argv[2]);
g1  = float(sys.argv[3]);
w2B = float(sys.argv[4]);

# compute some often used items
w0  = sqrt(w2B)
s1 = sqrt(n) * w0;
s2 = sqrt(n+1) * w0;

# construct the Hamiltonian
Hp  = mat([[-Vb/2, s1, 0, 0],
[s1, -Vb/2, g1, 0],
[0, g1, Vb/2, s2],
[0, 0, s2, Vb/2]])

# print result
if argc == 6:
        i   = int(sys.argv[5]);
        print eigh(Hp)[0][(i-1)%4]
elif argc == 5:
        print eigh(Hp)[0]
The desired result is now simply printed to the terminal and can be further processed by other commands. I am particular struck by how similar the relevant function names are to those of MATLab.

No comments: