On 8/3/07, Stephen Webb <[EMAIL PROTECTED]> wrote: > Greetings all, > > I've recently begun using Python to do scientific computation, and I wrote > the following script to find approximate eigenvalues for a semi-infinite > matrix: > > > from pylab import * > from numpy import * > from scipy import * > > def bandstructure(N,s): > > b = s/4.0 > > jmax = 10 + N**2 > > spectrum1 = [0]*2*N > spectrum2 = [0]*2*N > spectrum3 = [0]*2*N > > > for k in arange(1, 2*N+1, 1): > > A = zeros( (jmax,jmax) ) > > i = 0 > while i <= jmax-1: > if i <= jmax-2: > A[i,i+1] = b > A[i+1,i] = b > A[i,i] = ((k + 2.0*i*N)/N)**2 > i = i+1 > else: > A[i,i] = ((k + 2.0*i*N)/N)**2 > i = i+1 > > #This portion of the code builds a matrix twice as large to check > against > > B = zeros( (2*jmax,2*jmax) ) > > i = 0 > while i <= 2*jmax-1: > if i <= 2*jmax-2: > B[i,i+1] = b > B[i+1,i] = b > B[i,i] = ((k + 2.0*i*N)/N)**2 > i = i+1 > else: > B[i,i] = ((k + 2.0*i*N)/N)**2 > i = i+1 > > x = linalg.eigvals(A) > y = linalg.eigvals(B) > > j = 1 > while j<=3: > if abs(y[j]-x[j]) <= 10.0**(-5): > j = j + 1 > else: > print 'jmax not large enough to obtain accurate results' > > spectrum1[k-1] = x[0] + 0.5*s > spectrum2[k-1] = x[1] + 0.5*s > spectrum3[k-1] = x[2] + 0.5*s > > plot (k, spectrum1, k, spectrum2, k, spectrum3) > > xlabel('k (energy level)') > ylabel('E/E_r') > title('Finite Size Band Structure, N = %d, s = %f' % (N, s)) > grid(true) > > > When I run this script, I get the following message, which I can't figure > out: > > Traceback (most recent call last): > File "<stdin>", line 1, in <module> > File "bandstruc.py", line 61, in bandstructure > plot (k, spectrum1, k, spectrum2, k, spectrum3) > File > "/Library/Frameworks/Python.framework/Versions/2.5/lib/python2.5/site-packages/matplotlib/pylab.py", > line 2028, in plot > ret = gca().plot(*args, **kwargs) > File > "/Library/Frameworks/Python.framework/Versions/2.5/lib/python2.5/site-packages/matplotlib/axes.py", > line 2535, in plot > for line in self._get_lines(*args, **kwargs): > File > "/Library/Frameworks/Python.framework/Versions/2.5/lib/python2.5/site-packages/matplotlib/axes.py", > line 437, in _grab_next_args > for seg in self._plot_2_args(remaining[:2], **kwargs): > File > "/Library/Frameworks/Python.framework/Versions/2.5/lib/python2.5/site-packages/matplotlib/axes.py", > line 337, in _plot_2_args > x, y, multicol = self._xy_from_xy(x, y) > File > "/Library/Frameworks/Python.framework/Versions/2.5/lib/python2.5/site-packages/matplotlib/axes.py", > line 266, in _xy_from_xy > nrx, ncx = x.shape > ValueError: need more than 0 values to unpack >
Basically, you're trying to plot spectrum[123] vs. k: the spectra are lists, but k is an integer. I made some modifications to get things working more smoothly. Note the ion() function call at the beginning, and the use of "indx = arange(1,2*N+1)" instead of plotting vs. k. There's probably some other things that could be cleaned up, but I'm short on time right now; this will get your code working at least. Hope this helps, (From a fellow physicist) Kurt ################################################################## from pylab import * from numpy import * ion() def bandstructure(N,s): b = s/4.0 jmax = 10 + N**2 spectrum1 = [0]*2*N spectrum2 = [0]*2*N spectrum3 = [0]*2*N for k in arange(1, 2*N+1, 1): A = zeros( (jmax,jmax) ) i = 0 while i <= jmax-1: if i <= jmax-2: A[i,i+1] = b A[i+1,i] = b A[i,i] = ((k + 2.0*i*N)/N)**2 i = i+1 else: A[i,i] = ((k + 2.0*i*N)/N)**2 i = i+1 #This portion of the code builds a matrix twice as large to check against B = zeros( (2*jmax,2*jmax) ) i = 0 while i <= 2*jmax-1: if i <= 2*jmax-2: B[i,i+1] = b B[i+1,i] = b B[i,i] = ((k + 2.0*i*N)/N)**2 i = i+1 else: B[i,i] = ((k + 2.0*i*N)/N)**2 i = i+1 x = linalg.eigvals(A) y = linalg.eigvals(B) j = 1 while j<=3: if abs(y[j]-x[j]) <= 10.0**(-5): j = j + 1 else: print 'jmax not large enough to obtain accurate results' spectrum1[k-1] = x[0] + 0.5*s spectrum2[k-1] = x[1] + 0.5*s spectrum3[k-1] = x[2] + 0.5*s indx = arange(1,2*N+1, 1) plot (indx, spectrum1, indx, spectrum2, indx, spectrum3) xlabel('k (energy level)') ylabel('E/E_r') title('Finite Size Band Structure, N = %d, s = %f' % (N, s)) grid(True) -- http://mail.python.org/mailman/listinfo/python-list