p08: Eigenvalues of harmonic oscillator

-u'' + x^2 u = \lambda u, \qquad x \in \mathbb{R}

(1)

The exact eigenvalues are $1,2,5,7,\ldots$ .

#%config InlineBackend.figure_format='svg'
from numpy import pi,arange,linspace,sin,zeros,diag,sort,real
from scipy.linalg import toeplitz
from numpy.linalg import eig

L = 8.0
for N in range(6,37,6):
    h = 2.0*pi/N; x = h*linspace(1,N,N); x = L*(x-pi)/pi
    col = zeros(N)
    col[0] = -pi**2/(3.0*h**2) - 1.0/6.0
    col[1:] = -0.5*(-1.0)**arange(1,N)/sin(0.5*h*arange(1,N))**2
    D2 = (pi/L)**2 * toeplitz(col)
    evals,evecs = eig(-D2 + diag(x**2))
    eigenvalues = sort(real(evals))
    print("-------- N = %d --------" % N)
    for e in eigenvalues[0:4]:
        print("%20.14f" % e)

-------- N = 6 --------
    0.46147291699548
    7.49413462105052
    7.72091605300656
   28.83248377834010
-------- N = 12 --------
    0.97813728129860
    3.17160532064718
    4.45593529116678
    8.92452905811993
-------- N = 18 --------
    0.99997000149930
    3.00064406679582
    4.99259532440771
    7.03957189798150
-------- N = 24 --------
    0.99999999762904
    3.00000009841087
    4.99999796527329
    7.00002499815654
-------- N = 30 --------
    0.99999999999999
    3.00000000000074
    4.99999999997561
    7.00000000050862
-------- N = 36 --------
    0.99999999999998
    3.00000000000000
    4.99999999999999
    7.00000000000001