Strange set of errors

[Kurt Smith]

On 8/3/07, Stephen Webb <sdavis.webb at gmail.com> 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)