Software 8.4.1 (Harmonic functions). The fundamental solutions in two and three dimensions, the Poisson kernel in polar coordinates, and the total mass of the kernel at \(r = \frac 12\):

>>> import sympy as sp
>>> x, y, z = sp.symbols('x y z', real=True)
>>> r, th, ph = sp.symbols('r theta phi', positive=True)
>>> lap2 = lambda u: u.diff(x, 2) + u.diff(y, 2)
>>> lap3 = lambda u: u.diff(x, 2) + u.diff(y, 2) + u.diff(z, 2)
>>> sp.simplify(lap2(sp.log(x**2 + y**2)))
0
>>> sp.simplify(lap3(1 / sp.sqrt(x**2 + y**2 + z**2)))
0
>>> P = (1 - r**2) / (1 - 2*r*sp.cos(th - ph) + r**2) / (2*sp.pi)  # Poisson
>>> lap_polar = lambda u: u.diff(r, 2) + u.diff(r) / r + u.diff(th, 2) / r**2
>>> sp.simplify(lap_polar(P))
0
>>> sp.integrate(P.subs({r: sp.Rational(1, 2), ph: 0}), (th, 0, 2*sp.pi))
1

Reading the output. \(\log (x^2 + y^2)\) and \(1/|x|\) are harmonic away from the origin, the Poisson kernel is harmonic in \((r, \theta )\), and it integrates to \(1\), as the proof of Theorem 8.3.1 requires.

MATLAB/Octave: with the Symbolic Math Toolbox, simplify(laplacian(log(x^2+y^2), [x y])).

Software 8.4.2 (Poisson’s integral). The solution of Example 8.3.2 from Poisson’s formula by quadrature, against the closed form, at the centre and three other points (radius, angle):

>>> import numpy as np
>>> from scipy.integrate import quad
>>> def disc(g, rr, t):                    # unit disc, boundary values g(phi)
...     k = lambda p: ((1 - rr**2) / (2*np.pi)
...                    / (1 - 2*rr*np.cos(t - p) + rr**2))
...     return quad(lambda p: k(p) * g(p), 0, 2*np.pi, limit=200,
...                 points=[np.pi])[0]
...
>>> g = lambda p: 1.0 if p < np.pi else 0.0  # hot upper half, cold lower
>>> exact = lambda rr, t: 0.5 + np.arctan(2*rr*np.sin(t) / (1-rr*rr)) / np.pi
>>> for rr, t in ((0, 0), (0.5, np.pi/2), (0.9, np.pi/2), (0.5, -np.pi/4)):
...     print(rr, round(t, 4), round(disc(g, rr, t), 8),
...           round(exact(rr, t), 8))
...
0 0 0.5 0.5
0.5 1.5708 0.79516724 0.79516724
0.9 1.5708 0.96652458 0.96652458
0.5 -0.7854 0.25936746 0.25936746

Reading the output. The two agree to eight digits. The centre value is \(\frac 12\), the mean of the boundary values. At \((0.9, \frac \pi 2)\), close to the hot arc, \(u = 0.967\), and at \((0.5, -\frac \pi 4)\), in the cold half, \(u = 0.259\).

MATLAB/Octave: integral(@(p) (1-r^2)./(1-2*r*cos(t-p)+r^2)/(2*pi), 0, pi).

Software 8.4.3 (The five-point scheme). The square with one hot side, discretised with \(h = 1/200\): the five-point Laplacian on \(199^2 = 39\,601\) unknowns, assembled from Kronecker products and solved by sparse elimination. The centre value, then the series against the grid at \((0.5, 0.9)\), then the extreme values:

>>> import scipy.sparse as sps, scipy.sparse.linalg as spl
>>> n = 199  # interior points per side, h = 1/200
>>> h = 1 / (n + 1)
>>> T = sps.diags([-1, 2, -1], [-1, 0, 1], shape=(n, n))
>>> A = (sps.kron(sps.eye(n), T) + sps.kron(T, sps.eye(n))) / h**2   # -Delta
>>> b = np.zeros((n, n)); b[-1, :] = 1 / h**2      # u = 1 on top, 0 elsewhere
>>> U = spl.spsolve(A.tocsc(), b.ravel()).reshape(n, n)  # U[j, i]: (x_i, y_j)
>>> c = n // 2; U[c, c]                    # centre: exactly 1/4 by symmetry
0.24999999999992328
>>> k = np.arange(1, 4000, 2)  # series in sin(k pi x) sinh(k pi y)
>>> def series(X, Y):  # sinh ratio written to avoid overflow
...     ratio = (np.exp(k*np.pi*(Y - 1)) * (1 - np.exp(-2*k*np.pi*Y))
...              / (1 - np.exp(-2*k*np.pi)))
...     return np.sum(4 / (k*np.pi) * np.sin(k*np.pi*X) * ratio)
...
>>> series(0.5, 0.5), series(0.5, 0.9), U[179, c]  # y = 0.9 is row 179
(0.25000000000000006, 0.8016894653419548, 0.8016823459709765)
>>> U.max(), U.min()                       # maximum principle: inside (0, 1)
(0.9899251532058908, 2.735635261044178e-05)

Reading the output. The grid gives \(\frac 14\) at the centre to \(13\) digits, because the symmetry argument holds exactly for the discrete problem too. Near the hot side the grid value \(0.801682\) differs from the series value \(0.801689\) by \(7 \times 10^{-6}\), the \(O(h^2)\) error of the scheme. Inside the square the discrete solution lies strictly between \(0\) and \(1\): the discrete maximum principle, which follows from the five-point mean value rule just as Theorem 8.2.2 follows from Theorem 8.2.1.

MATLAB/Octave: A = kron(speye(n),T) + kron(T,speye(n)); U = A∖b.