9.4 Green’s functions in software

Software 9.4.1 (The string under load). The one-dimensional Green’s function applied to the loads \(f = 1\) and \(f = \sin \pi x\), with checks of the equation and the boundary values:

>>> import sympy as sp
>>> x, y = sp.symbols('x y', positive=True)
>>> def solve(f):          # -u'' = f on (0, 1), u = 0 at the ends: int G f dy
...     return sp.simplify(sp.integrate(y * (1 - x) * f(y), (y, 0, x))
...                        + sp.integrate(x * (1 - y) * f(y), (y, x, 1)))
...
>>> u1 = solve(lambda s: 1); u2 = solve(lambda s: sp.sin(sp.pi * s))
>>> sp.factor(u1), u2
(-x*(x - 1)/2, sin(pi*x)/pi**2)
>>> sp.simplify(-u2.diff(x, 2)), u2.subs(x, 0), u2.subs(x, 1)
(sin(pi*x), 0, 0)

Reading the output. A uniform load gives the parabola \(x(1 - x)/2\), and the load \(\sin \pi x\) gives \(\sin (\pi x)/\pi ^2\), since \(\sin \pi x\) is an eigenfunction with eigenvalue \(\pi ^2\); then \(-u'' = \sin \pi x\) and \(u\) vanishes at both ends.

MATLAB/Octave: integral(@(y) G(x,y).*f(y), 0, 1) with G as a function handle.

Software 9.4.2 (The half-plane and the induced charge). The half-plane Poisson integral with \(g = 1\) on \([-1, 1]\) and \(0\) elsewhere, against the closed form \(\frac 1\pi (\arctan \frac {x + 1}y - \arctan \frac {x - 1}y)\), the angle subtended by the segment divided by \(\pi \); then the total charge induced on a grounded plane by a unit charge at height \(1\):

>>> import numpy as np
>>> from scipy.integrate import quad
>>> def halfplane(g, X, Y, lo, hi):  # (1/pi) int Y g(t) / ((X - t)^2 + Y^2)
...     k = lambda s: Y / ((X - s)**2 + Y**2) / np.pi
...     return quad(lambda s: k(s) * g(s), lo, hi)[0]
...
>>> exact = lambda X, Y: (np.arctan((X+1) / Y) - np.arctan((X-1) / Y)) / np.pi
>>> for X, Y in ((0, 0.1), (0, 1), (2, 1), (0, 10)):   # g = 1 on [-1, 1]
...     print(X, Y, round(halfplane(lambda s: 1.0, X, Y, -1, 1), 10),
...           round(exact(X, Y), 10))
...
0 0.1 0.9365489651 0.9365489651
0 1 0.5 0.5
2 1 0.1475836177 0.1475836177
0 10 0.0634510349 0.0634510349
>>> h = 1.0                  # unit charge at height h over a grounded plane
>>> sigma = lambda r: -h / (2 * np.pi * (r**2 + h**2)**1.5)   # induced
>>> quad(lambda r: 2 * np.pi * r * sigma(r), 0, np.inf)[0]   # total charge
-1.0

Reading the output. The quadrature and the closed form agree to ten digits. Close to the segment, at \((0, 0.1)\), \(u = 0.937\); directly above it at height \(1\) the segment subtends a right angle and \(u = \frac 12\); far away, at height \(10\), \(u = 0.0635 \approx 2/10\pi \). The induced charge totals exactly \(-1\): every field line leaving the charge ends on the plane.

MATLAB/Octave: integral(@(r) -r./(r.^2+1).^1.5, 0, Inf).

Software 9.4.3 (The disc). The Green’s function of the unit disc with source at \(x = (0.5, 0)\): its values on the boundary, the solution of \(-\Delta u = 1\) at \(x\) from \(\int G(x, y)\,dy\), and the normal derivative of \(G\) on the boundary against the Poisson kernel:

>>> def G_disc(x, y):                      # unit disc, -Delta_y G = delta_x
...     xs = x / np.dot(x, x)              # image point x* = x / |x|^2
...     d1 = np.linalg.norm(y - x)
...     d2 = np.linalg.norm(x) * np.linalg.norm(y - xs)
...     return -np.log(d1 / d2) / (2 * np.pi)
...
>>> x0 = np.array([0.5, 0.0])
>>> [round(G_disc(x0, np.array([np.cos(a), np.sin(a)])), 12) + 0.0
...  for a in (0, 1, 2, 3)]           # G = 0 on the boundary circle
[0.0, 0.0, 0.0, 0.0]
>>> from scipy.integrate import dblquad    # -Delta u = 1: u = (1 - r^2) / 4
>>> point = lambda rho, a: rho * np.array([np.cos(a), np.sin(a)])
>>> val = dblquad(lambda rho, a: G_disc(x0, point(rho, a)) * rho,
...               0, 2 * np.pi, 0, 1, epsabs=1e-10)[0]   # u = int G dy
>>> round(val, 8), (1 - 0.25) / 4
(0.1875, 0.1875)
>>> eps = 1e-6                    # -dG/dnu on the boundary: Poisson kernel
>>> for a in (0.0, np.pi / 2, np.pi):
...     b = np.array([np.cos(a), np.sin(a)])
...     dn = -(G_disc(x0, b) - G_disc(x0, (1 - eps) * b)) / eps
...     P = (1 - 0.25) / (1 - 2 * 0.5 * np.cos(a) + 0.25) / (2 * np.pi)
...     print(round(a, 4), round(dn, 6), round(P, 6))
...
0.0 0.477465 0.477465
1.5708 0.095493 0.095493
3.1416 0.053052 0.053052

Reading the output. \(G\) vanishes on the circle, as the image construction guarantees. The integral of \(G\) gives \(0.1875\), which is \((1 - |x|^2)/4\), the solution of \(-\Delta u = 1\) with \(u = 0\) on the circle. The numerical normal derivative reproduces the Poisson kernel to six digits: \(3/2\pi = 0.477\) at the nearest boundary point, \(0.053\) at the farthest.

MATLAB/Octave: integral2 in polar coordinates with the same G.

Questions on this section

Stuck on something here? Ask below and it stays attached to this topic.