3 Green’s method – MATH 2811 Mathematical Methods II Weeks 6-10

3.1 Motivation for Green’s method

In this section we will devise an alternative approach to viewing and solving linear BVP, using the so-called Green’s function. Along the way we will encounter one of the most fundamental “functions” in mathematics, the Dirac delta function, the understanding of this mathematical object is critical in many fields of mathematics.

3.1.1 Form of the eigenfunction expansion solution

Consider the form of the final solution obtained through the eigenfunction expansion approach. With a little rearranging we can write the eigenfunction expansion solution found in the previous chapter as: \[ y(x) = \sum_{k=1}^\infty \frac{\langle f, w_k \rangle}{\lambda_k \langle y_k, w_k \rangle}y_k(x) \]

This requires all $\lambda_k\neq0$. The case of zero eigenvalue has two subcases: One where $\langle f,w_k\rangle\neq 0$ and one where this inner product is zero. In the former sub case, the boundary value problem has no solution, in the latter, the solution is not unique as you can add any multiple of $y_k$ (the eigenfunction that belongs to the zero eigenvalue) to the solution. This observation is directly linked to the Fredholm Alternative as we saw in the last chapter.

Let’s assume this is not an issue for now. Let $n_k = \langle y_k, w_k \rangle$ (normalisation), then: \[\begin{eqnarray*} y(x) &=& \sum_{k=1}^\infty \frac{1}{\lambda_k n_k} \left(\int_a^b f(t) w_k(t) \mathrm{d} t \right)y_k(x)\\ &=& \int_a^b \left( \sum_{k=1}^\infty \frac{1}{\lambda_k n_k} w_k(t) y_k(x) \right)f(t) \mathrm{d} t\\ &=& \int_a^b g(x,t) f(t) \mathrm{d} t \end{eqnarray*}\] where \[ g(x,t) = \sum_{k=1}^\infty \frac{w_k(t) y_k(x)}{\lambda_k n_k}. \tag{3.1}\]

Thus, we have constructed a solution to $Ly=f$ in the form \[ y(x)=\int_a^bg(x,t)f(t)\;dt. \tag{3.2}\] The function $g(x,t)$ is called the Green’s function (GF), and Equation 3.1 is an eigenfunction expansion of $g(x,t)$.

Of course, if we knew the Green’s function, we would have the solution without any need for the expansion, i.e. no need for the eigenfunctions. The goal in this section is to understand the properties of the GF and how to construct it.

Observe that if $L = L^*$, then $w_k = y_k$ and: In this case $g(x,t)=g(t,x)$, and we have the important connection between a self-adjoint operator and a symmetric Green’s function.

3.1.2 Inverse of differential operator

A nice way to think of the Green’s function is in terms of inverting the differential operator. Think about the familiar equation $\mathbf{A}\vec{x}=\vec{b}$ from linear algebra, to be solved for the unknown vector $\vec{x}$. The solution is given by \[\vec{x}=\mathbf{A}^{-1}\vec{b},\] i.e. we find the solution by multiplying the inverse of the linear operator (matrix) by the inhomogeneous term. Once you know the inverse operator, you can solve the problem for any given vector $\vec{b}$. In the context of BVP’s, $L$ is a differential operator, so it stands to reason that the inverse operator involves integration, hence the form of Equation 3.2. Constructing the Green’s function is analogous to finding the inverse of the matrix, once we have $g$ we can write down the solution Equation 3.2 for any forcing function $f(x)$.

3.1.3 An example

There are numerous ways to construct a Green’s function. We’ve already seen one: the eigenfunction expansion. Another way that you’ve probably seen before is via variation of parameters. This approach gives the Green’s function in a piecewise form.

Let’s take a simple example and look at the behaviour. Consider the BVP: \[ \begin{split} & Ly\equiv-y''=f(x), 0<x<1\\ &y(0)=y(1)=0 \end{split} \] We will temporarily magic the Green’s function out of thin air. \[ g(x,\xi)=\left\{\begin{array}{ll} (1-\xi)x\quad&0<x<\xi\\ (1-x)\xi\quad&\xi<x<1. \end{array}\right. \]

Let’s check this, so do so we use the Leibniz rule for differentiating integrals: \[ \frac{d}{dx}\int_{a(x)}^{b(x)}f(x,\xi)d\xi = \int_{a(x)}^{b(x)}\frac{d f(x,\xi)}{dx} + \frac{d b(x)}{dx}f(b,\xi) - \frac{\partial a(x)}{\partial x}f(a,\xi). \]

Tip 3.1: Show that $y(x)=\int_0^1 g(x,\xi)f(\xi)d\xi$ is a solution to this equation.

So we have \[\begin{align} y(x)& = \int_{0}^{1}g(x,\xi)f(\xi)d\xi = \int_0^{x}(1-x)\xi f(\xi) d\xi + \int_{x}^{1}(1-\xi)xf(\xi) d\xi,\\ \frac{d y}{dx} &= -\int_0^{x}\xi f(\xi) d\xi +(1-x)x f(x) + \int_{x}^{1}(1-\xi)f(\xi) d\xi - (1-x)x f(x) = -\int_0^{x}\xi f(\xi)d\xi +\int_{x}^{1}(1-\xi) f(\xi)d\xi,\\ \frac{d^2 y}{dx^2} &= -xf(x) - (1-x)f(x) = -f(x). \end{align}\] as required. Also note as the Green’s function vanishes at the boundaries it satisfies the required boundary conditions.

The following properties regarding the Green’s function itself are easily checked:

The GF satisfies $Lg=0$ if $x\neq\xi$
$g(x,\xi)$ satisfies the boundary conditions *as a function of $x*$.
$g$ is continuous on the whole interval $[0,1]$
$g$ is differentiable everywhere except at $x=\xi$, where it suffers a jump in the derivative.

These properties are in fact always true of the GF of a second order linear operator.

Note, however, that the function $y(x)$ satisfying $Ly=f$ is continuously differentiable assuming continuously differentiable $f$, meaning that the integration with $f(x)$ smooths out the discontinuity in $g$. To make sense of this, and to build some physical intuition, we shall need the notion of the delta function.

To give you a flavour, we pretend that the Green’s function is something that can be safely differentiated (rather than carefully differentiated by splitting the domain integral as in the above demonstration) i.e. we want to do this:

\[\begin{align} y(x)& = \int_{0}^{1}g(x,\xi)f(\xi)d\xi,\\ \frac{d y}{d x} &=\int_{0}^{1}\frac{\partial g(x,\xi)}{\partial x}f(\xi)d\xi,\\ \frac{d^2 y}{d^2 x} & =\int_{0}^{1}\frac{\partial^2 g(x,\xi)}{\partial x^2}f(\xi)d\xi \end{align}\]

For this to be a solution we need some “function” $\delta(x,\xi)$ satisfying \[ \int_{0}^1\delta(x,\xi)f(\xi)d \xi = f(x),\quad -\frac{\partial^2 g(x,\xi)}{\partial x^2} = \delta(x,\xi). \] We will come to call this the Green’s equation (actually it will later have total derivatives not partial ones).

But we should be (initially) skeptical the first derivative of the Green’s function is a step function, it doesn’t have a well defined derivative in the usual sense you saw in analysis I last year. So how can this $\delta$ exist?

3.2 Green’s function via delta function

To fix the context, consider stationary heat conduction in a rod:

\[ -y''(x)=f(x)\quad 0<x<1,\quad y(0)=0,\;\; y(1)=0 \]

where $y(x)$ is the temperature field and $f(x)$ is a given heat source density.

3.2.1 Delta function

The function $f(x)$ describes any heat added or removed from the system by the outside world. As a simple scenario, consider a point heat source, say located at the middle of the rod. Physically, this would correspond to applying heat at a single point only. How would we describe such a situation mathematically? What should we use for the function $f(x)$?

The notion of a point source is described by the ``delta function’’ $\delta$, characterised by properties

\[ \delta(x)=0\quad\forall x\neq0,\quad\quad\int_{-\infty}^\infty\delta(x)\;dx=1. \tag{3.3}\]

The first property captures the notion of a point function. The second property constrains the area under the curve (which you might think of as infinitely thin and infinitely high). This is an idealized point source at $x=0$, a point source at $x=a$ would be given by $\delta(x-a)$.

The problem is that no classical function satisfies Equation 3.3 (think: any function that is non-zero only at a point is either not integrable or integrates to zero).

3.2.2 Approximating the delta function

One way around this is to replace $\delta$ by an approximating sequence of increasingly narrower functions with normalized area, i.e. $f_n(x)$ where \[ \int_{-\infty}^{\infty} f_n(x)\mathrm{d}x=1\quad\forall n,\quad \lim_{n\rightarrow \infty} f_n(x)=0\quad\forall x\neq 0. \]

Example: ``hat’’ functions \[ f_n(x)=\left \{ \begin{array}{cc} 0\quad &\text{for } \left\vert x\right\vert>1/n\\ n/2\quad &\text{for } \left\vert x\right \vert \leq 1/n\end{array}\right. \tag{3.4}\] You can verify the $f_n(x)$ approach $\delta(x)$ as $n\to\infty$.

3.2.3 Properties of delta function

We have defined $\delta$ by Equation 3.3. We can use the approximating functions to obtain further properties.

Sifting property What happens when $\delta$ is integrated against another function?

Let $f(x)$ be a continuous function, and $F(x)=\int^xf(s) ds$ its antiderivative. Now consider approximating sequences:

\[\begin{align} \int_{-\infty}^{\infty}\delta(x-a)f(x)\mathrm{d}x &=\lim_{n\rightarrow\infty}\int_{-\infty}^{\infty}f_n(x-a)f(x)\mathrm{d}x \end{align}\]

and if $f_n$ are the hat functions Equation 3.4

\[\begin{align} &=\lim_{n\rightarrow\infty}\int_{a-1/n}^{a+1/n}\frac{n}{2}f(x)\mathrm{d}x=\lim_{n\rightarrow\infty}\frac{F(a+(1/n))-F(a-(1/n))}{2/n}\\ &=\lim_{s\rightarrow0}\frac{F(a+s)-F(a-s)}{2s}=F'(a)=f(a). \end{align}\]

Thus, we have \[ \int_{-\infty}^{\infty}\delta(x-a)f(x)\mathrm{d}x=f(a) \quad \text{ if $f$ is continuous at $a$. } \] In particular, \[ \int_{-\infty}^{\infty}\delta(x)f(x)\mathrm{d}x=f(0) \quad \text{ if $f$ is continuous at $x=0$. } \]

Thus, the delta function can be seen to sift out the value of a function at a particular point.

Antiderivative of $\delta(x)$. The antiderivative of the delta function is the so-called Heaviside function,

\[ \int_{-\infty}^x \delta(s)\mathrm{d}s= H(x)\equiv\left \{\begin{array}{ll} 0\quad &x<0\\ 1\quad &x>0.\end{array}\right. \tag{3.5}\]

Note that Equation 3.5 follows by integrating the sequence of approximating functions and showing that the limit is the Heaviside function. That is, if $H_n(x)=\int_{-\infty}^x f_n(s) ds,$ then $\lim_{n\rightarrow\infty} H_n(x)=H(x)$. (We leave this detail as an exercise!)

3.2.4 Point heat source

Let’s return to the heat conduction BVP with a point heat source of unit strength at the centre of the rod:

\[ -y''(x)=\delta(x-1/2), \quad0<x<1\label{u1}\quad y(0)=y(1)=0. \tag{3.6}\]

Since $\delta(x-1/2)=0\;\;\forall x\neq1/2$, this implies \[ -y''(x)=0, \quad0<x<1/2,\;1/2<x<1. \tag{3.7}\]

We can easily solve Equation 3.7 in each of the two separate domains $[0,1/2)$ and $(1/2,1]$ and then apply the BC of Equation 3.6. But be careful: there are two constants of integration for each domain, meaning four unknown constants total, and only two boundary conditions.

As you might expect (since $\delta(x-1/2)$ has vanished from Equation 3.7, the extra two conditions come in at $x=1/2$. To derive the extra conditions, imagine integrating equation Equation 3.6 across $x=1/2$: \[ \int_{1/2-}^{1/2+}-y''(x)\;dx=\int_{1/2-}^{1/2+}\delta(x-1/2)\;dx, \] where $1/2-$ ($1/2+$) signifies just to the left (right) of 1/2. Using property Equation 3.3 of the delta function, we have \[ -y']_{1/2-}^{1/2+}=1\quad\Rightarrow\quad y'(1/2+)-y'(1/2-)=-1. \tag{3.8}\] That is, the presence of the delta function defines a jump condition on $y'$.

Here, $y(\xi-)=\lim_{x\uparrow \xi}y(x)$, and $y(\xi+)=\lim_{x\downarrow \xi}y(x)$}

The other extra condition needed comes as a requirement that $y(x)$ is continuous across the point source, that is \[ y]_{1/2-}^{1/2+}=0. \tag{3.9}\]

More on this condition below. Solving Equations Equation 3.7, (with B.C’s Equation 3.6) along with extra conditions Equation 3.8 and Equation 3.9, we obtain the solution \[ u(x)=\left\{\begin{array}{ll} \;\;\;\frac{x}{2}\quad&0<x<1/2\\ -\frac{x}{2}+\frac{1}{2}\quad&1/2<x<1. \end{array}\right. \]

3.2.5 Green’s function construction

To motivate the construction of the Green’s function, consider the heat conduction problem with an arbitrary heat source:

\[ -y''(x)=f(x), \quad0<x<1,\quad y(0)=y(1)=0. \tag{3.10}\]

Imagine now describing $f$ by a distribution of point heat sources with varying strength; that is at point $x=\xi$ we imagine placing the point source $f(\xi)\delta(x-\xi)$.

The idea of the Green’s function is to introduce such an extra parameter $\xi$, and consider the system

\[ -g''(x,\xi)=\delta(x-\xi), \quad0<x<1\label{g1}\,\,\,g(0,\xi)=g(1,\xi)=0. \tag{3.11}\]

Note that prime denotes differentiation with respect to $x$, while $\xi$ is more like a place-holding variable. So, we have replaced $f(x)$ by a delta function, in order to solve for the Green’s function $g(x,\xi)$.

We have seen how to solve Equation 3.11, in the last section. The Green’s function is

\[ g(x,\xi)=\left\{\begin{array}{ll} (1-\xi)x\quad&0<x<\xi\\ (1-x)\xi\quad&\xi<x<1. \end{array}\right. \]

This is the solution I showed differentiating correctly at the start of this chapter in sec 3.1.3.

How to get back to the solution of Equation 3.10? For each $\xi$, the Green’s function gives the solution if a point heat source of unit strength were placed at $x=\xi$. Conceptually, then, to get the full solution we must ``add up’’ the point sources, scaled by the value of the heat source at each point: \[ y(x)=\int_0^1g(x,\xi)f(\xi)\;d\xi. \tag{3.12}\] To verify that this is indeed a solution, we can plug Equation 3.12 into Equation 3.10: \[ -y''(x)=\int_0^1-g''(x,\xi)f(\xi)\;dx=\int_0^1\delta(x-\xi)f(\xi)\;dx=f(x)\;\;\checkmark \]

This is the proposed derivative construction I mentioned was on shaky ground in section 3.1.3. We have now shown, if we are careful and stay inside integrals it is actually valid.

3.3 General linear BVP

We now consider a general $n$th order linear BVP with arbitrary continuous forcing function,

\[ Ly(x)=a_ny^{(n)}(x)+a_{n-1}y^{(n-1)}(x)+\dots+a_1y'(x)+a_0y(x)=f(x) \tag{3.13}\] for $a<x<b$, where each $a_i=a_i(x)$ is a continuous function, and moreover $a_n(x)\neq0$ $\forall x$. Along with Equation 3.13 are $n$ boundary conditions, each a linear combination of $y$ and derivatives up to $y^{(n-1)}$, evaluated at $x=a, b$. For instance, in the case $n=2$, the general form is:

\[ \begin{split} &B_1y\equiv \alpha_{11}y(a)+\alpha_{12}y'(a)+\beta_{11}y(b)+\beta_{12}y'(b)=\gamma_1\\ &B_2y\equiv \alpha_{21}y(a)+\alpha_{22}y'(a)+\beta_{21}y(b)+\beta_{22}y'(b)=\gamma_2. \end{split} \]

3.3.1 General Green’s Function

To solve Equation 3.13 with homogeneous BC \[B_iy=0,\;\;i=1\dots n-1,\] we first determine the Green’s function by solving

\[ \begin{split} &Lg(x,\xi)=\delta(x-\xi),\quad a<x<b\\ &B_ig=0. \end{split} \tag{3.14}\]

As before, \[Lg(x,\xi)=\delta(x-\xi)\] implies \[Lg(x,\xi)=0\quad\text{ on }a<x<\xi,\;\;\xi<x<b,\] i.e. we have a homogeneous problem to solve on two separate domains. As before, we require extra conditions, which come by integrating $Lg(x,\xi)=\delta(x-\xi)$ across $x=\xi$:

\[ \int_{\xi-}^{\xi^+}a_ng^{(n)}(x,\xi)+\dots+a_0g(x,\xi)\;d\xi=\int_{\xi-}^{\xi^+}\delta(x-\xi)\;dx. \]

The right hand side clearly integrates to one. If we were to perform an integration by parts on the first term of the left hand side, we would obtain \[ a_n(x)g^{(n-1)}(x,\xi)]_{\xi-}^{\xi+}+\int_{\xi-}^{\xi^+}(a_{n-1}-a_n')g^{(n-1)}+\dots+a_0g(x,\xi)\;dx=1. \] This equation is balanced by setting a jump condition on the $n-1$st derivative: \[ g^{(n-1)}(x,\xi)]_{\xi-}^{\xi+}=1/a_n(\xi), \] and taking all lower derivatives to be continuous across $x=\xi$: \[ g^{(j)}(x,\xi)]_{\xi-}^{\xi+}=0,\quad j=0,1,\dots n-2. \]

Once the Green’s function is determined, the solution to the BVP is given by \[ y(x)=\int_a^bg(x,\xi)f(\xi)\;d\xi. \]

3.4 Three-dimensional Green’s functions in (pretty) arbitrary domains.

We seek $G(x,\xi)$ such that \[ -\nabla_x^2 G(x,\xi)=\delta(x-\xi),\qquad G(x,\xi)\to0\ \text{as }|x|\to\infty. \] By translation and rotational symmetry, $G$ depends only on the distance $r=|x-\xi|$, i.e. $G(x,\xi)=\Phi(r)$ for some radial function $\Phi$.

Step 1: Solve away from the source.
For $r>0$ we have $\nabla_x^2 G=0$, so in spherical coordinates the radial Laplacian gives \[ \Phi''(r)+\frac{2}{r}\,\Phi'(r)=0. \] Integrating once: $r^2\Phi'(r)=A$ (constant), hence $\Phi'(r)=\dfrac{A}{r^2}$ and \[ \Phi(r)=-\frac{A}{r}+B. \] The decay condition $G\to0$ as $r\to\infty$ forces $B=0$, so for $r>0$ \[ G(x,\xi)=\Phi(r)=\frac{C}{r}\quad\text{with }C=-A. \]

Step 2: Fix the constant by normalization.
Integrate the defining equation over a ball $B_\varepsilon(\xi)$ of radius $\varepsilon$ centered at $\xi$ and use the divergence theorem: \[ \int_{B_\varepsilon(\xi)}\!\!-\nabla_x^2 G\,dV =\int_{\partial B_\varepsilon(\xi)}\!\!-\frac{\partial G}{\partial n}\,dS =\int_{B_\varepsilon(\xi)}\!\!\delta(x-\xi)\,dV=1. \] Since $G=C/r$, its radial derivative is $\dfrac{\partial G}{\partial r}=-\dfrac{C}{r^2}$. On the sphere $|x-\xi|=\varepsilon$, the outward normal derivative is $\dfrac{\partial G}{\partial n}=\dfrac{\partial G}{\partial r}=-\dfrac{C}{\varepsilon^2}$, hence \[ \int_{\partial B_\varepsilon(\xi)}\!\!-\frac{\partial G}{\partial n}\,dS =\int_{\partial B_\varepsilon(\xi)}\!\!\frac{C}{\varepsilon^2}\,dS =\frac{C}{\varepsilon^2}\cdot 4\pi \varepsilon^2 =4\pi C. \] Therefore $4\pi C=1$, so $C=\dfrac{1}{4\pi}$ and \[ G(x,\xi)=\frac{1}{4\pi\,|x-\xi|}. \]

Uniqueness.
Any other solution differing from $G$ by a harmonic function vanishing at infinity must be identically zero, so this $G$ is unique.

3.4.1 Arbitrary bounded domain: Green’s representation and boundary integrals

Let $\Omega\subset\mathbb{R}^3$ be a bounded domain with smooth boundary $\partial\Omega$.
For \[ -\nabla^2 u(x)=f(x)\ \text{in }\Omega, \] Green’s representation formula (with outward normal $n_\xi$ at $\xi\in\partial\Omega$) reads \[ u(x)=\int_{\partial\Omega}\!\Big(G(x,\xi)\,\frac{\partial u}{\partial n_\xi}(\xi)\;-\;u(\xi)\,\frac{\partial G}{\partial n_\xi}(x,\xi)\Big)\,dS_\xi \;+\;\int_{\Omega} G(x,\xi)\,f(\xi)\,dV_\xi, \] where $G$ is the Green’s function appropriate to the boundary condition.

Dirichlet problem ($u=g$ on $\partial\Omega$): choose the Dirichlet Green’s function $G_D$ solving \[ -\nabla_x^2 G_D(x,\xi)=\delta(x-\xi),\qquad G_D(x,\xi)=0\ \text{for } x\in\partial\Omega. \] Then \[ u(x)=-\int_{\partial\Omega} g(\xi)\,\frac{\partial G_D}{\partial n_\xi}(x,\xi)\,dS_\xi \;+\;\int_{\Omega} G_D(x,\xi)\,f(\xi)\,dV_\xi,\qquad x\in\Omega. \] Hence $u$ is obtained from known boundary data $g$ and the volume source $f$.
Neumann problem $\big(\frac{\partial u}{\partial n}=h$ on $\partial\Omega\big)$: choose the Neumann function $G_N$ solving \[ -\nabla_x^2 G_N(x,\xi)=\delta(x-\xi),\qquad \frac{\partial G_N}{\partial n_x}(x,\xi)=0\ \text{for } x\in\partial\Omega, \] together with the usual compatibility $\int_\Omega f\,dV=\int_{\partial\Omega} h\,dS$. Then \[ u(x)=\int_{\partial\Omega} G_N(x,\xi)\,h(\xi)\,dS_\xi \;+\;\int_{\Omega} G_N(x,\xi)\,f(\xi)\,dV_\xi,\qquad x\in\Omega, \] which uses only the known Neumann data $h$ and $f$.

Remark (arbitrary shapes).
For a general $\partial\Omega$, $G_D$ or $G_N$ is rarely available in closed form. In practice one represents $u$ using the free-space kernel and unknown boundary densities, e.g. the single-layer potential \[ u(x)=\int_{\partial\Omega} \frac{1}{4\pi|x-\xi|}\,\sigma(\xi)\,dS_\xi \quad\text{or}\quad u(x)=\int_{\partial\Omega} \frac{\partial}{\partial n_\xi}\!\left(\frac{1}{4\pi|x-\xi|}\right)\mu(\xi)\,dS_\xi, \] and enforces the boundary condition to obtain an integral equation for $\sigma$ or $\mu$ on $\partial\Omega$.
This reduces the PDE in an arbitrary domain to a boundary-only problem.