Green’s Functions Illustrated

There’s no shame in admitting this: Green’s functions can be a tad bit confusing. This post aims to provide a visual and descriptive overview of what a Green’s function is and how it can help solve certain types of differential equations, namely nonhomogeneous boundary value problems in one dimension. Green’s functions can be used for many other types of problems too, but it’s often helpful to start with a specific problem before looking at the general theory.

Review of Terminology

A differential equation is said to be linear if the unknown function and its derivatives appear only to the first power, multiplied by coefficients that depend only on the independent variable. For example, $x^2u^{”}(x) + u(x) = 0$ is linear since the coefficients do not depend on $u$, as opposed to $2u u^{”}(x) + u(x)= 0$ which is not (that first coefficient includes a $u$ term). In this case, the differential equation is said to be nonlinear.

Linear differential equations are usually represented using linear operators denoted $L$. For example, the differential equation above can be written using the operator $L = x^2 {d^2 \over dx^2} + 1$, for which $Lu = x^2u^{”}(x) + u(x) = 0$. In words, “$L$ operates on $u$.”

A linear differential equation is said to be homogeneous if every term involves the unknown function or its derivatives, e.g., $u^{”}(x) + u(x) = 0$, as opposed to $u^{”}(x) + u(x) = 1$ or $u^{”}(x) + u(x) = \sin(x)$. In operator form, homogeneity is expressed as $Lu=0$. An equation that’s not homogeneous is called nonhomogeneous. The nonhomogeneous term in a differential equation is called the forcing function and often represents something that acts on the system like an external force or heat source. In operator form, this is written $Lu=f$, where $f$ is the forcing function.

Linear equations obey superposition: if $u$ solves one equation and $v$ solves another with the same operator, then any linear combination of $u$ and $v$ solves the corresponding linear combination of the equations. In operator form:

If $L(u) = f$ and $L(v) = g$  then  $ L(\alpha u + \beta v) = \alpha L(u) + \beta L(v) = \alpha f + \beta g$

In order to solve a differential equation we’ll need boundary conditions (or initial conditions if time is the independent variable). The number of conditions required is the same as the order of the differential equation. For example, a second-order differential equation will require two boundary conditions. These boundary conditions describe how the function of interest, $u(x)$, behaves at its boundaries. They can also be homogeneous or nonhomogeneous depending on the value of $u$ along them. If both the differential equation and the boundary conditions are homogeneous, then the trivial solution, $u=0$, always works. Often it’s the only solution, but not always: $u^{”} + \pi^2 u = 0$ with $u(0) = u(1) = 0$ is also solved by $u = \sin(\pi x)$. As it turns out, the homogeneous problem having only the trivial solution is exactly the condition for a Green’s function to exist.

The domain is the region $\Omega$ on which $u$ is defined, together with its boundary, denoted $\partial \Omega$.

Motivation and Overview

Referring to the diagram above, Green’s functions aid in solving nonhomogeneous linear differential equations with nonhomogeneous boundary conditions (the top gray box).

Because the equation of interest is linear, the problem can be split into two parts: a nonhomogeneous linear differential equation with homogeneous boundary conditions and a homogeneous linear differential equation with nonhomogeneous boundary conditions. See the diagram above.

These are then summed together or superimposed to form the overall solution.

The Green’s function does the heavy lifting on the first of these problems: nonhomogeneous linear differential equation with homogeneous boundary conditions. When the boundary values are prescribed (Dirichlet conditions), the second problem is the well-known Dirichlet Problem. As we’ll see at the end, the Green’s function can help with that part too. If there are no nonhomogeneous boundary conditions, the Green’s function provides the entire solution.

Before we get into the solution technique, it’s useful to understand the history. Green’s functions were first introduced in 1828 by George Green, a self-taught mathematician and miller from Nottingham, in a privately published essay on the mathematics of electricity and magnetism. In modern language, the idea at the heart of his method is this: what if the potential produced by a continuous charge distribution could be represented as an infinite sum of point sources, each of whose effects already respects the boundary? Let’s start there.

The Point Source

As the name suggests, a point source is one that exists only at a single point yet imparts a measurable effect across the whole domain. In electrostatics, this would be a point charge, a fictitious object that possesses a finite charge and no size. Other examples are the point force in mechanics and the point heat source in heat transfer.

Because these point sources have zero size, the concentration (charge or force per unit length, area or volume) they produce at that point is infinite. Note that the total amount of the source is finite, and its effect on the system is not limited to the point itself. Its effects linger as you move away from the source (see figure). In one dimension, that effect stays finite everywhere; in two and three dimensions it blows up at the source point itself, just as the potential of a point charge does, but remains finite everywhere else.

This is the idea behind a Dirac Delta function, denoted $\delta(x – \xi)$, where $\xi$ is the location of the point source. You should note its appearance in the diagram above too, written there as shorthand $LG=\delta(\xi)$. In this expression, the concentration of the point source is acting as a forcing function. So depending on the differential equation under consideration, this concentration can be interpreted to be a load (mechanics, see below), a heat source (steady state heat equation), a charge density (Poisson’s equation) and so on.

Next, let’s impose homogeneous boundary conditions on this point source and require its effects to taper to zero at the boundaries. The combined expression in operator form is:

$$\begin{eqnarray}
\begin{cases}
\phantom{-}
LG(x, \xi) = \delta(x – \xi) && x \in \Omega
\\\phantom{-}
G(x, \xi) = 0 && x \in \partial \Omega
\end{cases}
\end{eqnarray}$$

The solution to this pair of equations is called the Green’s function $G(x, \xi)$.

Conceptually a Green’s function measures the effect that this tapering point source induces at any point in the domain, analogous to how the actual solution $u$ in that domain responds to the forcing function of interest. A Green’s function is sometimes called an influence function for this reason.

A Green’s function measures the effect of a point source at each point in the domain

It’s important to observe that the Green’s function is independent of any forcing function and nonhomogeneous boundary conditions on this domain. It only cares about the differential equation (per the operator $L$), the boundary geometry, and the type of boundary condition (fixed values versus fixed slopes, for example). As will be shown, this provides a baseline solution for applying arbitrary forcing functions and boundary conditions in this domain. Thus a Green’s function provides the kernel for solving arbitrary problems with the same DE and shape.

There are two important caveats. First, since this technique requires the principle of superposition, it’s only applicable when the corresponding differential equation is linear, as with Poisson’s equation, the wave equation or the heat equation. For nonlinear differential equations, other techniques will be required.

Second caveat. Most problems don’t give up a Green’s function easily. And when they do, it’s often in the land of toy problems, regions whose boundary geometry is highly symmetric, like a sphere or a rectangle. But in two dimensions, for Laplace’s and Poisson’s equations, those toy geometries can be mapped to more complex regions in a way that preserves the solution using conformal mappings. More on that in another post perhaps. For now, let’s jump into a classic example to better understand what all these words mean.

The following canonical one-dimensional Green’s function problem was adapted from [SG]§1.1 and [GBG]§I.4. Some of the formulas will be provided, not derived, to keep things moving along. The mathematical detail will be provided at the end.

A cable of length $l$ is stretched such that it has an axial tension $T$, i.e., it doesn’t sag. We assume for the sake of simplicity that the cable is massless and thus does not experience any gravitational forces (i.e., toy problem). To make it even clearer, we’ll set $l=T=1$.

Next, subject this cable to a continuous transverse load denoted $P(x)$ as shown below (the red arrows). This is measured in units of force per unit length, thus the terms “load” and “pressure” are used interchangeably. The load causes displacement in the cable, denoted $u(x)$. Our domain is the bounded region $[0, 1]$ and the boundaries are the two endpoints. The forcing function is the load.

It can be shown [MF] that, for small displacements, the relationship between load and displacement is described by a second-order linear differential equation with homogeneous boundary conditions—the cable is fixed on each end. Specifically, it’s a one-dimensional Poisson’s equation.

$$\begin{eqnarray}-T{d^2 u\over dx^2} = P(x) \,\, \,\, x \in (0,1)\end{eqnarray}$$ $$\begin{eqnarray}u(0) = u(1) = 0\end{eqnarray}$$

Relating this back to the first problem diagram above, our linear differential equation has homogeneous boundary conditions and so there’s no Dirichlet problem to solve. This means the Green’s function will provide the solution to the entire problem, assuming we can find it.

To find the Green’s function, we modify the problem. Instead of a continuous load $P(x)$, imagine there’s a point load applied at a single point $\xi$, denoted $\delta(x – \xi)$. Recall that load is force per length and that the length of a point is zero. And also recall that a point load is one that’s infinite at this point and zero everywhere else. The resulting force created by this load (remember, load = force/length) is found by integrating across the interval, $F= \int_0^1 \delta(s – \xi)ds$, which by definition is one (“unit load”). Also mentioned earlier, these properties correspond to the delta distribution, which is loosely defined as

$$\begin{eqnarray}
\delta(\xi) = \begin{cases}
\phantom{-} 0 & \text{if } \xi \ne 0
\\\phantom{-} +\infty & \text{if } \xi=0
\end{cases}
\end{eqnarray}$$

$$\begin{eqnarray}
\int^{\infty}_{-\infty} \delta(\xi) d\xi = 1
\end{eqnarray}$$

The delta represents the load of a point source: It produces a total force of one unit, concentrated at a single point of the cable. The corresponding displacement at a point $x$ from this source point at $\xi$ is measured by a Green’s function, denoted $G(x, \xi)$.

This Green’s function displacement adheres to the same differential equation as the continuous load, except here the point load is characterized by a delta distribution. And like the continuous load, the deflection is zero at the endpoints.

$$\begin{eqnarray}-T{d^2 \over dx^2} G(x, \xi)= \delta(x – \xi) \,\, \,\,\,\, x \in (0,1)\end{eqnarray}$$ $$\begin{eqnarray}G(0, \xi) = G(1, \xi) =0\end{eqnarray}$$

On either side of $\xi$, any linear function, $Ax + B$, will satisfy this differential equation (with different constants on each side). We also observe that the derivative won’t exist at $\xi$–the delta distribution represents a singularity at that point. But we want our solution to be continuous too (no break in the cable), despite the kink at that one point. In fact, integrating the equation across $\xi$ shows that the slope must drop by exactly $1/T$ there; this jump, together with continuity and the boundary conditions, pins down the solution. Finally, since the endpoints are zero, the resulting displacement should increase from zero when $x < \xi$ and decrease toward it when $x > \xi$. In other words, it’s an inverted V-shaped function. A graph of such a function is shown below.

Observe that it looks like a thin wire being stretched by a screwdriver, as expected. The exact Green’s function for this problem is described by the following piecewise function, which is derived in [SJ] and others. Note that this generalizes the problem to a cable of length $l$ and tension $T$:

$$\begin{eqnarray}G(x, \xi) = {x (l – \xi ) \over Tl} & x \le \xi\end{eqnarray}$$ $$\begin{eqnarray} \space\space\space\space\space\space\space\space\space\space= {\xi (l – x) \over Tl} & \xi \le x \end{eqnarray}$$

If you rearrange things a bit, an interesting property is revealed [SJ]:

$$\begin{eqnarray}G(x, \xi) = \frac{1}{T}(1 – \frac{\xi}{l})x \space\space\space\space\space\space x \le \xi\end{eqnarray}$$ $$\begin{eqnarray} \space\space\space\space\space\space\space\space\space\space = \frac{1}{T}(1 – \frac{x}{l})\xi \space\space\space\space\space\space \xi \le x \end{eqnarray}$$ $$\begin{eqnarray}
=G(\xi, x) \end{eqnarray} $$

Turns out the Green’s function for this problem is symmetric: $G(x, \xi) = G(\xi, x)$. This symmetry is a consequence of the operator being self-adjoint (with these boundary conditions), and it will be a useful feature once we start solving equations with a Green’s function.

Important sidebar on units
The units matter here. Per the expression above, the numerator is the product of two distance values and the denominator is the product of force and distance. In total that’s distance per unit force. This means that the displacement at $x_0$ from an arbitrary force, $F_{\xi}$, applied at the point $\xi$, is just $F_{\xi} * G(x_0, \xi)$. Why? Because the Green’s function is normalized for force. To repeat,

The displacement at $x_0$ induced by a force at the point $\xi$ is the Green’s function $G(x_0, \xi)$ multiplied by the magnitude of that force.

Now let’s see if we can think about the whole picture in terms of this weird Green’s function.

Because the governing equation is linear (Poisson’s equation), we can take advantage of the principle of superposition, which says that these effects add linearly: the effect on a point is the sum of all of the forces acting on it.

Start by dividing the cable into $n$ equal increments of length $l/n$, the center of which is $\xi_i$. The intervals themselves will be denoted $\Delta \xi_i$. Now within each interval, define a Green’s function whose singularity is located at $\xi_i$. Pay attention to this point: we are now modeling the system using a collection of Green’s functions, not just one. This is shown below.

Next, we approximate the force applied across this interval as the product of the pressure and the interval length: force = pressure * length.

$$\begin{eqnarray} F_{\xi_i} = P(\xi_i) {l \over n} = P(\xi_i) \Delta \xi_i \end{eqnarray}$$

With this, we can now estimate the displacement at point $x_0$ caused by a segment as $G(x_0, \xi_i) P(\xi_i) \Delta \xi_i$.

Finally, by superposition we can estimate the total displacement at $x_0$ caused by the load function $P$ by summing  across the complete set of $n$ intervals.

$$\begin{eqnarray} u(x_0) \approx \sum_{i=1}^{n} G(x_0, \xi_i) P(\xi_i) \Delta \xi_i \end{eqnarray}$$

Taking the limit as those intervals shrink to zero we get an integral expression for displacement:

$$\begin{eqnarray}u(x_0) = \int_{\Omega}G(x_0, \xi) P(\xi)  d\xi\end{eqnarray}$$

The displacement is thus the linear superposition of individual Green’s functions weighted by the load function. Make sense? This was an a-ha moment for me so let it sink in a second.

Verification

The above arguments were a bit hand-wavy, as mathematicians say. But they were important because they provided a heuristic for what we’re trying to accomplish. We can now formally solve the original differential equation and see how it aligns with the results above.

Start by integrating both sides of this equation twice to get the following double integral:

$$\begin{eqnarray}-Tu(x) = \int_0^x \int_0^\xi P(s)dsd\xi + C_2x + C_1\end{eqnarray}$$

Then switch the order of integration. The region of integration is $0 \le s \le \xi \le x$, so after the switch $\xi$ runs from $s$ to $x$. I have to comment that this switch is not obvious, so it may help to graph the region and compare the two integrands. A great discussion of this can be found in [SL]§9.1. The main motivation is that the load function is now constant with respect to the inner integration variable. (In the last step, the dummy variable $s$ is renamed $\xi$.)

$$\begin{eqnarray}=  \int_0^x \int_{s}^x P(s)d\xi ds + C_2x + C_1\end{eqnarray}$$ $$\begin{eqnarray} = \int_0^x (x-\xi) P(\xi)d\xi  + C_2x + C_1\end{eqnarray}$$

The boundary conditions determine the values of the constants. $C_1=0$ since $u(0)=0$. Similarly, $u(l)=0$ gives

$$\begin{eqnarray}C_2 = {-1 \over l}\int_0^l(l-\xi)P(\xi)d\xi \end{eqnarray}$$

Now plug this in, split the region of integration between $(0, x)$ and $(x, l)$, and divide through by $-T$. Grouping terms yields:

$$\begin{eqnarray}u(x) = {1 \over Tl}
\int_0^x \xi(l – x) P(\xi)d\xi + {1 \over Tl}\int_x^l x(l – \xi) P(\xi)d\xi \end{eqnarray}$$

Hey check it out: the two integrands are just our Green’s function, evaluated on either side of $x$! More generally we have:

$$\begin{eqnarray} u(x) = \int_0^l G(x, \xi) P(\xi) d\xi \end{eqnarray}$$

Which by the symmetry property is also

$$\begin{eqnarray} u(x) = \int_0^l G(\xi, x) P(\xi) d\xi \end{eqnarray}$$

And so we’ve arrived at the same result as with our hand-wavy approach earlier. We could have also obtained a more general expression that includes nonzero boundary conditions if we just left the values as $u(l)$ and $u(0)$. Here’s how the equation turns out with nonhomogeneous boundary conditions. Referring back to our original problem diagram, this solution is now incorporating the Dirichlet solution too.

$$\begin{eqnarray} u(x) = u(0)\frac{l – x}{l} + u(l)\frac{x}{l} + \int_0^l G(x, \xi) P(\xi) d\xi \end{eqnarray}$$

Indeed these new boundary terms have a relationship with $G(x,\xi)$ too (hint: its derivative) but we’ll save that for later.

What’s Next?

This foundational material can be generalized to $n$-dimensional problems, which is the realm of partial differential equations. The math gets harder but fortunately the conceptual ideas remain the same. You will often see such problems in introductory courses in electrostatics and heat transfer.

If you enjoyed this post please click the like button, share it, or leave a comment. It would mean a lot considering how much effort went into writing it. Thanks!

References

[MF] Fowler, M. Physics 152 class notes, University of Virginia. 2002.

[GBG] Greenberg, M. Applications of Green’s Functions in Science and Engineering. Dover. 1971, 2015.

[SG] Stakgold, I., Boundary Value Problems of Mathematical Physics, Macmillan. 1967.

[SJ] Johnson, S., Math 18.303 class notes. MIT. 2011.

[SL] Shearer, M. Levy, R., Partial Differential Equations: An Introduction to Theory and Applications. Princeton University Press. 2015