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 functions is and how it can help solve certain types of differential equations, namely inhomogeneous 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 all of its coefficients are continuous functions of only 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 written in operator form is $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 its constant terms are zero, 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.

The solution to any linear combinations of linear differential equations is the corresponding linear combination of their individual solutions. This is called superposition. 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 degree of the differential equation. For example, a second degree 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. Of course, if both the differential equation and the boundary conditions are homogeneous then the only solution is the trivial solution, $u=0$.

The domain is the set of points in the closure of the interior region $\Omega$ on which $u$ is defined. Its boundary is 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 handles the first of these problems: nonhomogeneous linear differential equation with homogeneous boundary conditions. The second problem is the well-known Dirichlet Problem. If there are no inhomogeneous 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 proposed in an 1828 paper on electrostatics by amateur physicist George Green. In this paper, he made the following proposition: what if the potential produced by a continuous charge distribution could be somehow represented as an infinite sum of point sources? 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 length, the concentration (or pressure) it produces at that point is infinite: it occupies an area of zero (i.e., force per unit area). Note that the effect on the system from that point source is not infinite. And it’s not limited to the point itself. Its effects lingers as you move away from the source (see figure).

This is the idea behind a Dirac Delta function, denoted $\delta(\xi)$, where $\xi$ is the location of the point source. You should note its appearance in the diagram above too: $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 flux (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) && u \in \Omega
\\\phantom{-}
G(x, \xi) = 0 && u \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. Green’s functions are 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$) and its boundary geometry. 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. For example, 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 those toy geometries can be translated 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 the formulas will be provided not derived to keep things moving along. The mathematical detail will be provided at the end.

A cable of uniform density and length of $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 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) \,\, \,\, u \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(\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 P(s)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 force of one unit across the entire 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}$$

Any linear function, $Ax + B$, will satisfy this differential equation except at $\xi$. 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 (which here just means connected), despite the differential discontinuity at one point. 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 of $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 (a.k.a., self-adjoint): $G(x, \xi) = G(\xi, x)$. This 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}$, at the same location, 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 solution is a linear differential equation (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$, 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_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 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. To do so, observe the change in the integrand values that results. 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 integrand.

$$\begin{eqnarray}=  \int_0^x \int_{\xi}^x P(\xi)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,

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

Now plug this in and split the region between $(0, \xi)$ and $(\xi, l)$. 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! More generally we have:

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

Which by the symmetry property is also

$$\begin{eqnarray} u(x) = \int_0^l G(x, \xi) 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 non-zero boundary conditions if we just left the values as $u(l)$ and $u(0)$. Here’s how the equation turns out with inhomogeneous boundary conditions. Referring back to our original problem diagram, this solution is now incorporating the Dirichlet solution too.

$$\begin{eqnarray} u(x) = u(0)(l – x) + u(l)x + \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