Software 9.6.1 (Eigenvalues by shooting). For \(-y'' = \lambda (1 + t)y\), \(y(0) = y(1) = 0\), integrate from \(y(0) = 0\), \(y'(0) = 1\) and find the values of \(\lambda \) at which \(y(1) = 0\):

>>> import numpy as np
>>> from scipy.integrate import solve_ivp, quad
>>> from scipy.optimize import brentq
>>> w = lambda t: 1 + t                    # -y'' = lam w y, y(0) = y(1) = 0
>>> def end(lam):                          # shoot from y(0) = 0, y'(0) = 1
...     f = lambda t, z: [z[1], -lam * w(t) * z[0]]
...     return solve_ivp(f, (0, 1), [0, 1], rtol=1e-11, atol=1e-12).y[0, -1]
...
>>> grid = np.linspace(1, 260, 400); vals = [end(l) for l in grid]
>>> shoot = [brentq(end, a, b) for a, b, fa, fb
...          in zip(grid, grid[1:], vals, vals[1:]) if fa * fb < 0]
>>> np.round(shoot[:4], 4)
array([  6.5484,  26.4649,  59.6742, 106.17  ])
>>> S = quad(lambda t: np.sqrt(w(t)), 0, 1)[0]          # WKB: (n pi / S)^2
>>> np.round([(n * np.pi / S)**2 for n in range(1, 5)], 2)
array([  6.64,  26.57,  59.78, 106.28])

Reading the output. The first four eigenvalues are \(6.5484\), \(26.4649\), \(59.6742\) and \(106.17\). The last two lines compare them with the large-\(n\) approximation \(\lambda _n \approx (n\pi /S)^2\), \(S = \int _0^1\sqrt {w}\), of the Liouville–Green (WKB) method: \(6.64, 26.57, 59.78, 106.28\). The approximation is too high by about \(0.1\) for every \(n\), so its relative error falls like \(1/n^2\). Shooting turns the boundary value problem into a sequence of initial value problems, using the existence theory of the earlier sections, and the sign changes of \(y(1)\) bracket the eigenvalues.

MATLAB/Octave: bvp4c with an unknown parameter, or fzero on the shooting function as here.

Software 9.6.2 (Finite differences and the Rayleigh bound). Replace \(-y''\) by the second difference on \(399\) interior points, which turns the problem into the generalised matrix eigenvalue problem \(Ky = \lambda Wy\) with \(W = \operatorname {diag}(1 + t_i)\):

>>> from scipy.linalg import eigh
>>> N = 400; h = 1 / N; t = np.linspace(h, 1 - h, N - 1)   # interior points
>>> K = (np.diag(2 * np.ones(N - 1)) - np.diag(np.ones(N - 2), 1)
...      - np.diag(np.ones(N - 2), -1)) / h**2           # -y'', y = 0 at ends
>>> lam_fd = eigh(K, np.diag(w(t)), eigvals_only=True)[:4]
>>> np.round(lam_fd, 4)
array([  6.5484,  26.4644,  59.6713, 106.161 ])
>>> u = lambda t: t * (1 - t)              # Rayleigh trial function
>>> num = quad(lambda t: (1 - 2 * t)**2, 0, 1)[0]      # integral of u'^2
>>> R = num / quad(lambda t: w(t) * u(t)**2, 0, 1)[0]
>>> R, shoot[0]
(6.666666666666667, 6.548395305999766)

Reading the output. The matrix eigenvalues agree with the shooting values to four or five digits, the error growing with \(n\) like \(h^2\lambda _n^2\), as the second difference has error \(O(h^2y'''')\). The symmetric matrix \(K\) is the discrete form of the self-adjoint operator, and its eigenvectors are the discrete eigenfunctions drawn in Figure 13. The Rayleigh quotient of the trial function \(t(1 - t)\) is \(6.667\), above \(\lambda _1 = 6.548\) by \(1.8\%\), as Theorem 9.4.1 requires.

MATLAB/Octave: eig(K, W) for the generalised eigenvalue problem.

Software 9.6.3 (An eigenfunction expansion). The sine series of \(f(t) = t(1 - t)\) on \([0, 1]\), the expansion in the eigenfunctions of \(-y'' = \lambda y\), \(y(0) = y(1) = 0\):

>>> f = lambda t: t * (1 - t)              # sine series on [0, 1]
>>> b = [2 * quad(lambda t: f(t) * np.sin(n * np.pi * t), 0, 1)[0]
...      for n in range(1, 8)]
>>> np.round(b, 6)
array([ 0.258012,  0.      ,  0.009556, -0.      ,  0.002064,  0.      ,
        0.000752])
>>> [round(8 / (n * np.pi)**3, 6) if n % 2 else 0 for n in range(1, 8)]
[0.258012, 0, 0.009556, 0, 0.002064, 0, 0.000752]
>>> s3 = lambda t: sum(bn * np.sin((n + 1) * np.pi * t)
...                    for n, bn in enumerate(b[:3]))
>>> max(abs(s3(x) - f(x)) for x in np.linspace(0, 1, 201))
0.0029594198519884707

Reading the output. The coefficients \(b_n = 2\int _0^1f(t)\sin n\pi t\,dt\) agree with the exact values \(8/(n\pi )^3\) for odd \(n\) and \(0\) for even \(n\). They decay like \(n^{-3}\), fast, because \(f\) satisfies the boundary conditions and is smooth, and three terms already give a maximum error of \(0.003\) against a maximum value of \(f\) of \(0.25\). A function that did not vanish at the ends would have coefficients decaying only like \(1/n\), and its series would converge slowly and show the Gibbs overshoot.

MATLAB/Octave: integral(@(t) 2*t.*(1 - t).*sin(n*pi*t), 0, 1).