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.