p08: Eigenvalues of harmonic oscillator

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

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

%config InlineBackend.figure_format='svg'
from numpy import pi,arange,linspace,sin,zeros,diag,sort
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(evals)
    print("-------- N = %d --------" % N)
    for e in eigenvalues[0:4]:
        print("%20.14f" % e)

-------- N = 6 --------
    0.46147291699547
    7.49413462105052
    7.72091605300656
   28.83248377834012
-------- N = 12 --------
    0.97813728129861
    3.17160532064719
    4.45593529116679
    8.92452905811993
-------- N = 18 --------
    0.99997000149931
    3.00064406679583
    4.99259532440772
    7.03957189798150
-------- N = 24 --------
    0.99999999762904
    3.00000009841086
    4.99999796527328
    7.00002499815654
-------- N = 30 --------
    0.99999999999998
    3.00000000000074
    4.99999999997560
    7.00000000050862
-------- N = 36 --------
    0.99999999999999
    3.00000000000000
    5.00000000000000
    6.99999999999999