p08: Eigenvalues of harmonic oscillator¶

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

In [1]:

%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

In [2]:

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("%24.15e" % e)

N = 6
   4.614729169954720e-01
   7.494134621050521e+00
   7.720916053006558e+00
   2.883248377834012e+01
N = 12
   9.781372812986103e-01
   3.171605320647185e+00
   4.455935291166790e+00
   8.924529058119932e+00
N = 18
   9.999700014993121e-01
   3.000644066795828e+00
   4.992595324407724e+00
   7.039571897981499e+00
N = 24
   9.999999976290406e-01
   3.000000098410861e+00
   4.999997965273280e+00
   7.000024998156542e+00
N = 30
   9.999999999999772e-01
   3.000000000000737e+00
   4.999999999975596e+00
   7.000000000508624e+00
N = 36
   9.999999999999927e-01
   2.999999999999998e+00
   5.000000000000000e+00
   6.999999999999988e+00