Software 12.6.1 (Rayleigh and Ritz). The Rayleigh quotient of \(x(1 - x)\) on \((0, 1)\), and the Ritz method with the two trial functions \(x(1 - x)\) and \(x^2(1 - x)^2\):

>>> import sympy as sp
>>> x = sp.symbols('x')
>>> def R(u):                              # Rayleigh quotient on (0, 1)
...     top = sp.integrate(u.diff(x)**2, (x, 0, 1))
...     return top / sp.integrate(u**2, (x, 0, 1))
...
>>> R(x * (1 - x)), sp.N(sp.pi**2, 8)      # one trial function
(10, 9.8696044)
>>> phi = [x * (1 - x), x**2 * (1 - x)**2]   # Ritz: best combination of two
>>> I = lambda g: sp.integrate(g, (x, 0, 1))
>>> K = sp.Matrix(2, 2, lambda i, j: I(phi[i].diff(x) * phi[j].diff(x)))
>>> M = sp.Matrix(2, 2, lambda i, j: I(phi[i] * phi[j]))
>>> lam = sp.symbols('lam')
>>> roots = sp.solve((K - lam * M).det(), lam)
>>> sorted(sp.N(r, 10) for r in roots)     # against pi^2 and 9 pi^2
[9.869749621, 102.1302504]

Reading the output. One trial function gives \(10\); two give \(9.869750\), within \(1.5 \times 10^{-5}\) of \(\pi ^2 = 9.869604\), and always above it. Both trial functions are symmetric about \(\frac 12\), so the second Ritz value, \(102.13\), is an upper bound for the next symmetric eigenvalue, \(9\pi ^2 = 88.83\), not for \(4\pi ^2\).

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

Software 12.6.2 (Linear finite elements). The finite element method for \(-u'' = \pi ^2\sin \pi x\), whose solution is \(\sin \pi x\), assembled element by element with five-point Gauss quadrature; the \(L^2\) error, the energy-norm error, and the largest error at the nodes, for \(N\) elements:

>>> import numpy as np
>>> gx, gw = np.polynomial.legendre.leggauss(5)  # Gauss points on [-1, 1]
>>> def fem(N, f):  # -u'' = f on (0, 1), linear elements
...     h = 1 / N; nodes = np.linspace(0, 1, N + 1)
...     K = (np.diag(2*np.ones(N - 1)) - np.diag(np.ones(N - 2), 1)
...          - np.diag(np.ones(N - 2), -1)) / h
...     b = np.zeros(N - 1)
...     for e in range(N):                 # load vector, element by element
...         xq = nodes[e] + h * (gx + 1) / 2; wq = gw * h / 2
...         for loc, node in ((1 - (xq - nodes[e]) / h, e - 1),
...                           ((xq - nodes[e]) / h, e)):
...             if 0 <= node < N - 1:
...                 b[node] += np.sum(wq * f(xq) * loc)
...     return nodes, np.r_[0, np.linalg.solve(K, b), 0]
...
>>> def errors(nodes, U, u, du):  # L2 and energy-norm errors, by Gauss
...     eL2 = eH1 = 0.0
...     for e in range(len(nodes) - 1):
...         a, h = nodes[e], nodes[e + 1] - nodes[e]
...         xq = a + h * (gx + 1) / 2; wq = gw * h / 2
...         Uq = U[e] + (U[e + 1] - U[e]) * (xq - a) / h
...         eL2 += np.sum(wq * (u(xq) - Uq)**2)
...         eH1 += np.sum(wq * (du(xq) - (U[e + 1] - U[e]) / h)**2)
...     return np.sqrt(eL2), np.sqrt(eH1)
...
>>> f = lambda s: np.pi**2 * np.sin(np.pi * s)        # exact u = sin(pi x)
>>> for N in (8, 16, 32, 64):
...     nodes, U = fem(N, f)
...     e0, e1 = errors(nodes, U, lambda s: np.sin(np.pi*s),
...                     lambda s: np.pi * np.cos(np.pi*s))
...     print(N, f"{e0:.3e}", f"{e1:.3e}",
...           f"{np.abs(U - np.sin(np.pi*nodes)).max():.1e}")
...
8 9.921e-03 2.512e-01 2.2e-16
16 2.487e-03 1.258e-01 2.8e-16
32 6.220e-04 6.295e-02 2.8e-15
64 1.555e-04 3.148e-02 9.1e-15

Reading the output. Each doubling of \(N\) divides the \(L^2\) error by \(4\) and the energy error by \(2\): orders \(2\) and \(1\), as predicted. The nodal errors are at rounding level. This is special to one dimension: the Green’s function \(G(x_i, \cdot )\) of Section 9 is piecewise linear with a corner at \(x_i\), so it lies in \(V_h\), and Galerkin orthogonality then gives \(u_h(x_i) = u(x_i)\) exactly.

MATLAB/Octave: the same assembly loop; K = (2*eye(N-1) - diag(ones(N-2,1),1) - diag(ones(N-2,1),-1))/h.

Software 12.6.3 (Weierstrass’s example). The value of \(\int _{-1}^1x^2y'^2\,dx\) for \[ y_\varepsilon (x) = \frac {\arctan (x/\varepsilon )}{\arctan (1/\varepsilon )} , \] which satisfies \(y(\pm 1) = \pm 1\):

>>> from scipy.integrate import quad
>>> def weierstrass(eps):                  # E[y] = int_{-1}^{1} x^2 y'^2 dx
...     a = np.arctan(1 / eps)             # y = arctan(x/eps)/arctan(1/eps)
...     dy = lambda s: eps / (eps**2 + s**2) / a
...     g = lambda s: s**2 * dy(s)**2
...     return quad(g, -1, 1, points=[0], limit=200)[0]
...
>>> [f"{weierstrass(e):.3e}" for e in (0.1, 0.01, 0.001, 0.0001)]
['6.340e-02', '6.366e-03', '6.366e-04', '6.366e-05']

Reading the output. The values fall like \(2\varepsilon /\pi \), so the infimum is \(0\), but it is not attained: the minimising sequence steepens into a step at the origin.

MATLAB/Octave: integral(@(x) x.^2.*dy(x).^2, -1, 1).