3 Green’s method

3.1 Motivation for Greens 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 eigen exapnsion 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 \(<f,w_k>\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.

Lets 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 involve 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. \]

Lets 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(x,b) - \frac{\partial a(x)}{\partial x}f(x,a). \]

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

So we have \[\begin{align} y(x)& = \int_{0}^{1}g(x,\xi)f(x)dx = \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 Greens function itself are easily checked:

1. The GF satisfies \(Lg=0\) if \(x\neq\xi\)
2. \(g(x,\xi)\) satisfies the boundary conditions as a function of \(x\).
3. \(g\) is continuous on the whole interval \([0,1]\)
4. \(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 Greens 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(x)dx,\\ \frac{d y}{d x} &=\int_{0}^{1}\frac{\partial g(x,\xi)}{\partial x}f(x)dx,\\ \frac{d^2 y}{d^2 x} & =\int_{0}^{1}\frac{\partial^2 g(x,\xi)}{\partial x^2}f(x)dx \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 Greens 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 \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 \(\eqref{u1}\) 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 \(\eqref{delta}\) 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 \[ y(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\,\,\,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

First we solve 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 x =\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. \] (the integral domain is vanishingly small and will vanish if the integrand is continuous).

Once the Green’s function is determined, the solution to the Homogeneous BVP is given by \[ y(x)=\int_a^bg(x,\xi)f(\xi)\;d\xi. \] Why homogeneous? We see the point is the Greens function (and its derivatives) going to zero on the boundary means \(y\) will also. This is not generally true for inh

3.4 Three dimension Green’s functions.

3.4.1 Green’s Function in 2D for the Laplacian.

We would like to solve Poisson’s equation in 2D: \[ \nabla^2 u(\vec{x}) = f(\vec{x}). \] for \(\vec{x}\in\mathbb{R}^2\): the free space problem. To do so we solve the Green’s problem in this domain: \[ \nabla^2 G(\vec{x}, \vec{\xi}) = \delta(\vec{x} - \vec{\xi}) \] in free space, with \(\vec{x},\vec{\xi}\in\mathbb{R}^2\). Note the derivatives in the Laplacian are with respect to \(x\), and also that it is translation invariant. If we make the change in coordinates \(\vec{y} = \vec{x} - \vec{\xi}\) then we have: \[ \nabla^2_{\vec{y}} G(\vec{y}) = \delta(\vec{y}). \] This allows us to centre the problem at the origin and adopt polar coordinates \((r,\theta)\), centered at origin. The Laplacian in polar coordinates for a radially symmetric function is: \[ \nabla^2_{\vec{y}} G = \frac{1}{r} \frac{d}{dr} \left( r \frac{dG}{dr} \right) \] Why radially symmetric? The \(\delta\) function must be independent of the angular parameter \(\theta\). As we saw in chapter 2 the solution to this equation would have a Fourier series form in the \(\theta\) coordinate, so we cannot have such behaviour.

Away from the singularity at \(r = 0\), the homogeneous equation is: \[ \frac{1}{r} \frac{d}{dr} \left( r \frac{dG}{dr} \right) = 0. \] Integrating once: \[ r \frac{dG}{dr} = C_1,\rightarrow \frac{dG}{dr} = \frac{C_1}{r}. \] Integrating again: \[ G(r) = C_1 \ln r + C_2. \] We set \(C_2 = 0\), we have freedom to do this as Poisson’s equation is only defined up to a constant with Dirichlet B.C’s., we get: \[ G(r) = C_1 \ln r. \] To determine \(C_1\), integrate over a disk of radius \(\epsilon\), this is the equivalent of the jump conditions: \[ \oint_{\partial B_\epsilon} \nabla G \cdot dS = 1. \] Since \(\nabla G = \frac{C_1}{r} \hat{r}\), we get: \[ 2\pi C_1 = 1 \quad \Rightarrow \quad C_1 = \frac{1}{2\pi}. \] Thus, on substituting back \(\vec{y} = \vec{x}-\vec{\xi}\) and using the fact that \(r = |\vec{y}|= |\vec{x}-\vec{\xi}|\), we obtain the Green’s function in 2D to be: \[ G(\vec{x}, \vec{\xi}) = \frac{1}{2\pi} \ln |\vec{x} - \vec{\xi}|. \]

Then the solution to the full equation is \[ u(\vec{x}) = \int_{\mathbb{R}^2}G(\vec{x}, \vec{\xi})f(\vec{\xi})\mathrm{d}V. \]