10 Green's functions for PDEs In this final chapter we will apply the idea of Green's functions to PDEs, enabling us to solve the wave equation, diffusion equation and Laplace equation in unbounded domains. We will also see how to solve the inhomogeneous (i.e. forced) version of these equations, and uncover a relationship, known as Duhamel's principle, between these two classes of problem. In our construction of Green's functions for the heat and wave equation, Fourier transforms play a starring role via the 'differentiation becomes multiplication' rule. We derive Green's identities that enable us to construct Green's functions for Laplace's equation and its inhomogeneous cousin, Poisson's equation. We conclude with a look at the method of images one of Lord Kelvin's favourite pieces of mathematical trickery. 10.1 Fourier transforms for the heat equation Consider the Cauchy problem for the heat equation = D2 t (10.1) on Rn [0, ), where D is the diffusion constant and where obeys the conditions |Rn {0} = f (x) and lim (x, t) = 0 t . |x| (10.2) Taking the Fourier transform w.r.t. the spatial variables we have t) (k, t) = D|k|2 (k, t (10.3) 0) = f(k). This equation has a unique solution satisfying the with initial condition (k, initial conditions, given by t) = f(k) eD|k|2 t . (k, (10.4) Our solution (x, t) itself is the inverse Fourier transform of the product of f(k) with the 2 Gaussian eD|k| t , and the convolution theorem tells us that this will be the convolution of the initial data f (x) = F 1 [f] with the inverse Fourier transform of the Gaussian. On the third problem sheet you've already shown that, for a Gaussian in one variable k22 2 2 F[ea x ] = e 4a (10.5) a and therefore, setting a2 = 1/4Dt and treating t as a fixed parameter in performing the 2 2 Fourier transforms, we find F 1 [eDtk ] = ex /4Dt / 4Dt in one dimension, or ! " 1 |x|2 1 D|k|2 t F [e= exp =: Sn (x, t) (10.6) 4Dt (4Dt)n/2 in n dimensions, since the n-dimensional Gaussian is a product of n Gaussians in independent variables, each of which may be (inverse) Fourier transformed separately. The function Sn (x, t) is known as the fundamental solution of the diffusion equation. - 126 - Putting everything together, the general solution to our Cauchy problem is the convolution " ! # 1 |x y|2 (x, t) = (f Sn )(x, t) = dn y (10.7) f (y) exp 4Dt (4Dt)n/2 Rn Notice that, provided f is integrable, this function indeed obeys 0 as |x| . Also |x|2 note the appearance of the dimensionless similarity variable 2 = 4Dt that we previously introduced. For example, suppose instead we take the initial condition itself to be a Gaussian f (x) = normalised so that & Rn $ a %n/2 2 0 ea|x| , (10.8) f (x) dn x = 0 . In this case, our general solution (10.7) gives ! " # |x y|2 a %n/2 2 exp a|y| dn y (x, t) = 0 4 2 Dt 4Dt Rn ! " $ a %n/2 # (1 + 4aDt)|y|2 + 2x y |x|2 exp dn y = 0 4 2 Dt 4Dt Rn ' ( (2 ) ! "# $ a %n/2 ( a|x|2 (1 + 4aDt) (( x ( = 0 exp exp y dn y ( 2 4 Dt 1 + 4aDt 4Dt 1 + 4aDt ( n $ = 0 ! a 1 2 1 + 4aDt "n/2 ! exp R a|x|2 1 + 4aDt "# 2 Rn where in going to the last line we've substituted * ! " 1 + 4aDt x z= y and thus 4Dt 1 + 4aDt Thus finally (x, t) = 0 ! a/ 1 + 4aDt "n/2 e|z| dn z , (10.9) n d y= ! ! 4Dt 1 + 4aDt a|x|2 exp 1 + 4aDt " . "n/2 dn z . (10.10) Thus an initial Gaussian retains a Gaussian form, with its squared width (1 + 4aDt)/a spreading linearly with t (recall linear growth of variance for diffusing probabilistic processes) while the total area remains constant (c.f. the conservation law we obtained for the heat equation in chapter 4) and the peak at x = 0 drops as t1/2 . This is as we saw in chapter 4, see in particular figure 8. As a limiting case, suppose the initial profile Gaussian (10.8) becomes more and more sharply peaked around x = 0, retaining the same total amount of heat 0 . Then in the limit f (x) = 0 n (x) and our solution (10.7) becomes simply ! " # 0 |x y|2 (n) (x, t) = (y) exp dn y 4Dt (4Dt)n/2 Rn (10.11) = 0 Sn (x, t) . - 127 - In other words, the fundamental solution is the solution (up to a constant factor) when the initial condition is a -function. For all t > 0, the -pulse spreads as a Gaussian. As t 0+ we regain the function as a Gaussian in the limit of zero width while keeping the area constant (and hence unbounded height). A striking property of this solution is that || > 0 everywhere throughout Rn for any finite t > 0, no matter how small. Thus, unlike for the wave equation, disturbances propagate via the heat equation arbitrarily rapidly the presence of the spike of heat at the origin when t = 0 affects the solution at arbitrarily great distances (even if only by exponentially suppressed terms) immediately. Mathematically, this is a consequence of the fact that the heat equation is parabolic, and so has only one family of characteristic surfaces (in this case, they are the surfaces t = const.). Physically, we see that the heat equation is not compatible with Special Relativity; once again this is because it is really just a macroscopic approximation to the underlying statistical mechanics of microscopic particle motion. 10.1.1 The forced heat equation In the previous section we solved the homogeneous equation with inhomogeneous boundary conditions. Here instead suppose solves the inhomgeneous (forced) heat equation D2 = F (x, t) t on Rn [0, ), but now with homogeneous boundary conditions: |Rn {0} = 0 and lim = 0 t . |x| (10.12) (10.13) Note that the general case the inhomogeneous equation with inhomogeneous boundary conditions can be reduced to these two cases: We can write the solution as = h + f where h satisfies the homogeneous equation with the given inhomogeneous boundary conditions while f obeys the forced equation with homogeneous boundary conditions. (Such a decomposition will clearly apply to all the other equations we consider later.) Turning to (10.12), we seek a Green's function G(x, t; y, ) such that G(x, t; y, ) D2 G(x, t; y, ) = (t ) (n) (x y) (10.14) t and where G(x, 0; y, ) = 0 in accordance with our homogeneous initial condition. Given such a Green's function, the function # # (x, t) = G(x, t; y, ) F (y, ) dn y d , (10.15) 0 Rn solves (at least formally) the forced heat equation (10.12) using the defining equation (10.14) of the Green's function and the definition of the -functions. To construct the Green's function, again take the Fourier transform of (10.14) w.r.t. the spatial variables x, recalling that F[ (n) (x y)] = eiky . After multiplying through 2 by eDk t we obtain , + D|k|2 t 2 e G(k, t; y, ) = eiky+D|k| t (t ) (10.16) t - 128 - subject to the initial condition G(k, 0; y, ) = 0. This equation is easily solved, and we find that the Fourier transform of the Green's function is # t 2 ikyD|k|2 t G(k, t; y, ) = e eD|k| u (u ) du 0 0 t< (10.17) = 2 (t ) ikyD|k| e t> = (t ) eikyD|k| 2 (t ) , where (t ) is the step function. Finally, taking the inverse Fourier transform in the k variables we recover the Green's function itself as # (t ) 2 G(x, t; y, ) = eik(xy) eD|k| (t ) dn k . (10.18) n (2) Rn We recognise the integral here as the (inverse) Fourier transform of a Gaussian exactly as in (10.6), but where the resulting function has arguments x% = x y and t% = t . Thus we immediately have G(x, t; y, ) = (t ) Sn (x y, t ) where (10.19) ! " |x% |2 1 exp Sn (x , t ) = 4Dt% (4Dt% )n/2 % % is the same fundamental solution we encountered in studying the homogeneous equation with inhomogeneous boundary conditions. 10.1.2 Duhamel's principle The fact that the same function Sn (x, t) appeared in both the solution to the homogeneous equation with inhomogeneous boundary conditions, and the solution to the inhomogeneous equation with homogeneous boundary conditions is not a coincidence. To see this, let's return to (10.7) but now suppose we impose |Rn { } = f (x) (10.20) at time t = (rather than at t = 0). A simple time translation of (10.7) shows that for times t > , this problem is solved by # h (x, t) = f (y) Sn (x y, t ) dn y . (10.21) Rn Thus, the initial data f (x) at t = has been propagated using the fundamental solution for a time t up to time t. On the other hand, using the causality constraint from the step function, the solution (10.15) for the forced problem takes the form / # t .# n f (x, t) = F (y, ) Sn (x y, t ) d y d . (10.22) 0 Rn - 129 - Now suppose that for each fixed time t = , we view the forcing term F (y, t) as an initial condition (!) imposed at t = . The integral in square brackets above represents the effect of this condition propagated to time t as in (10.21). Finally, the time integral in (10.22) expresses the solution f to the forced problem as the accumulation (superposition) of the effects from all these conditions at times earlier than t, each propagated for time interval t up to time t. The upper limit t of this integral arose from the step function (t ) in the Green's function (10.19) and expresses causality: the solution at time t depends only on the cumulative effects of 'initial' conditions applied at earlier times < t. Conversely, an 'initial condition' applied at time cannot influence the past (t < ). The relation between solutions to homogeneous equations with inhomogeneous boundary conditions and inhomogeneous equations with homogeneous boundary conditions is known as Duhamel's principle. Below, we'll see it in action in the wave equation, too. 10.2 Fourier transforms for the wave equation In this section, we'll apply Fourier transforms to the wave equation. Suppose F : Rn [0, ) is a given forcing term and that solves the forced wave equation 2 c 2 2 = F t2 (10.23) on Rn (0, ), subject to the homogeneous initial conditions that |Rn {0} = 0 , t |Rn {0} = 0 and lim = 0 , |x| (10.24) and where c is the wave (phase) speed. As usual we seek a Green's function G(x, t; y, ) that obeys the simpler problem ! 2 " 2 2 c G(x, t; y, ) = (t ) (n) (x y) (10.25) t2 where the forcing term is replaced by -functions in all time and space variables, and where the Green's function obeys the same boundary conditions G(x, 0; y, ) = 0 , t G(x, 0; y, ) = 0 and lim G(x, t; y, ) = 0 |x| (10.26) as itself. Once again, if we can find such a Green's function then our problem (10.23) is solved by # # (x, t) = F (y, ) G(x, t; y, ) dn y d , (10.27) 0 Rn at least formally. Proceeding as before, we Fourier transform the defining equation for the Green's function wrt the spatial variables x, obtaining 2 G(k, t; y, ) + |k|2 c2 G(k, t; y, ) = eiky (t ) . t2 - 130 - (10.28) and t G vanish at t = 0. Treating Again, this is subject to the conditions that both G k as a parameter, we recognize this equation as an initial value problem in the variable t essentially the same as the example given previously at the end of section 7.5. As we found there, The solution is sin (|k|c(t )) G(k, t; y, ) = (t ) eiky , c|k| (10.29) as we found there. To recover the Green's function itself we must compute the inverse Fourier transform # (t ) sin |k|c(t ) n G(x, t; y, ) = d k (10.30) eik(xy) n (2) c|k| Rn Unlike the case of the heat equation, here the form of the Green's function is sensitive to the number of spatial dimensions. To proceed, for definiteness we'll pick n = 3, which is physically the most important case. We introduce polar coordinates for k chosen so that |k| = k and k (x y) = k|x y| cos That is, we choose to align the z-axis in k-space along whatever direction x y points. Thus, in three spatial dimensions the Green's function becomes # (t ) sin kc(t ) 3 G(x, t; y, ) = eik|xy| cos d k 3 (2) ck R3 / # .# 1 (t ) ik|xy| cos = e d(cos ) sin kc(t ) k dk (10.31) (2)2 c 0 1 . # / (t ) 1 ik|xy| = e sin kc(t ) dk 2ic|x y| 2 where in going to the third line we used that fact that result of integrating over cos leaves us with an even function of k. We recognise the remaining integral (the contents of the square bracket in the last line) as the inverse Fourier transform of sin kc(t ), where the variable conjugate to k is |x y|. From (8.43) we have that F 1 [sin ka] = i ((x + a) (x a)) 2 (10.32) and therefore our Green's function is (t ) [(|x y| + c(t )) (|x y| c(t ))] 4c|x y| 1 = (|x y| c(t )) 4c|x y| G(x, t; y, ) = (10.33) where we note that the first -function has no support in the region where the step function is non-zero, and that in the second line the argument of the -function already ensures t > 0 so the step function is superfluous. The fact that only this term contributes to our -function can be traced to the boundary conditions we imposed in (10.26). The Green's function here is often called the retarded Green's function. - 131 - Using this Green's function, our solution can finally be written as # # 1 F (y, ) (x, t) = (|x y| c(t )) dn y d 4c Rn |x y| 0 # 1 F (y, tret ) n = d y 2 4c Rn |x y| (10.34) where the retarded time tret is defined by ctret := ct |x y| . (10.35) The retarded time illustrates the finite propagation speed for disturbances that obey the wave equation: the effect of the forcing term F at some point y is felt by the solution at a different location x only after time t tret = |x y|/c has elapsed. This is how long it took the disturbance to travel from y to x, and is in perfect agreement with what we expect from our discussion of characterstics for the wave equation. Above we solved the forced wave equation in 3 + 1 dimensions. As an exercise, I recommend you show that the general solution to the forced wave equation with homogeneous boundary conditions in 1 + 1 dimensions is # t# (c(t ) |x y|) (x, t) = F (y, ) dy d 2c 0 R (10.36) # # 1 t x+c(t ) = F (y, ) dy d . 2c 0 xc(t ) (Really have a go! it's not entirely trivial.) Thus, the retarded Green's function in 1+1 dimensions is a step function, which in the second line we used to restrict the region over which we take the y-integral. Using the step function this way means that the order of the remaining integrations is important if we are to capture the correct, -dependent domain of influence (i.e. y [x c(t ), x + c(t )]) of the forcing term. We can also see Duhamel's principle in operation for the wave equation here. Comparing to equation (9.45), we see that for t > 1 2c # x+c(t ) F (y, ) dy xc(t ) is simply d'Alembert's solution to the homogeneous (unforced) wave equation in the case that the (inhomogeneous) initial data is given by (x, ) = 0 and t (x, ) = F (x, ) , (10.37) with these conditions imposed when t = . Hence the solution (10.36) of the forced equation with homogeneous initial data can again be viewed as a superposition of influences from the forcing term F used as an initial condition (here on the time-derivative of ) for each < t. - 132 - 10.3 Poisson's equation Poisson's equation for a function : Rn C is the forced Laplace equation 2 = F (x) (10.38) where as always the forcing term F : Rn C is assumed to be given, and the minus sign on the rhs of (10.38) is conventional. We will again solve this equation by first constructing a Green's function for it, although the solution via Green's functions is less straightforward than for either the heat or wave equations in particular, we will not obtain it via Fourier transforms (the Green's function itself will turn out not to be absolutely integrable). We also note that, since there is no time variable (the equation is elliptic), the solutions will have no notion of causality. 10.3.1 The fundamental solution The fundamental solution Gn to Poisson's equation in n-dimensions is defined to be the solution to the problem 2 Gn (x; y) = (n) (x y) , (10.39) where the forcing term is chosen to be just an n-dimensional -function. Since the problem rotationally symmetric about the special point y, the fundamental solution can only depend on the scalar distance from that point: Gn (x; y) = Gn (|x y|) . (10.40) Integrating both sides of (10.39) over a ball Br = {x Rn : |x y| r.} of radius r centred on y we obtain # # 2 1= Gn dV = Br Br n Gn dS (10.41) where the integral of the -function gives 1, and we've used the divergence theorem to reduce the integral of 2 Gn to an integral over the boundary sphere Br = {x Rn : |xy| = r}. When n = 3, this is the surface of a 'standard' sphere, while it is a circle when n = 2 and a 'hypersphere' for n > 3. In each case, the outward pointing normal n points in the radial direction, so n Gn = dGn /dr. Thus we obtain # # dGn n1 n1 dGn r dn = r 1= dn (10.42) dr S n1 Br dr where we have used the fact that G is a function only of the radius, and where dn is the measure on a unit S n1 ; for example sin d d when n=3 dn = . d when n=2 - 133 - The remaining angular integral just gives the total 'surface area' An of our (n 1)- dimensional sphere. In particular, A2 = 2 is the circumference of a circle, whereas A3 = 4 is the total solid angle of a sphere in three dimensions39 . Then dGn 1 = . dr Arn1 (10.43) This ode is easily solved to show that the fundamental solution is |x y| + c1 n=1 1 n=2 Gn (x; y) = 2 log |x y| + c2 1 1 + cn n 3 An (n 2) |x y|n2 (10.44) When n 3 we can fix the constant cn to be zero by requiring limr G(r) = 0; however, this condition cannot be imposed in either one or two dimensions, and some other condition (usually suggested by the details of a specific problem) must instead be imposed. Note that when n > 1 these fundamental solutions satisfy the Laplace equation for all x Rn {y}, but that they become singular at x = y. Thus we shall need to exercise special care when using them as Green's functions for the Poisson equation throughout Rn . 10.3.2 Green's identities Green's identities40 establish a relationship between the fundamental solutions Gn and solutions of Poisson's equation. First, suppose Rn is a compact set with boundary , and let , : R be a pair of functions on that are regular throughout . Then the product rule and divergence theorem give41 # # () dV = 2 + () () dV # (10.45) = (n ) dS where n is the outward pointing normal to . This is Green's first identity. Simply interchanging and we likewise find # # # 2 () dV = + () () dV = (n ) dS , (10.46) and subtracting this equation from the previous one we arrive at # # 2 2 dV = (n ) (n ) dS , (10.47) ! You may enjoy showing that An = 2 (n1)/2 /((n 1)/2) in general, where (t) = 0 xt1 ex dx is Euler's Gamma function. We won't need this result, though. 40 Green's identities are sometimes called Green's theorems, but they're not to be confused with the two-dimensional version of the divergence theorem! 41 I'm using the notation dV to denote the standard Cartesian measure dn x on , and dS to denote the measure on the (n 1)-dimensional boundary of . 39 - 134 - which is Green's second identity. In order to make use of this result, we'd like to apply it to the case that = Gn , our fundamental solution. However, Gn (x, y) is singular at y (at least when n > 1) so it's not clear that the derivation of (10.47) is valid, because the divergence theorem usually requires functions to be regular throughout . Let's give a more careful derivation. Let B$ be a ball of radius ) - 1 centred on the dangerous point y, and let be the region = Br B$ = {x Rn : ) |x y| r} . (10.48) Now Gn is perfectly regular everywhere within , so we can apply Green's second identity with = Gn to find # # 2 2 Gn Gn dV = Gn 2 dV # (n Gn ) Gn (n ) dS = (10.49) # # = (n Gn ) Gn (n ) dS + (n Gn ) Gn (n ) dS Srn1 S!n1 where the first equality follows since 2 Gn = 0 in , and we've included contributions from both boundary spheres (with radii r and )) in the final line. On the inner boundary a sphere of radius ) the outward-pointing unit normal is n = r and we have 1 1 An (n 2) )n2 1 1 = . An )n1 Gn |inner bdry = n Gn |inner bdry (10.50) The measure on an (n 1)-dimensional sphere of radius ) is dS = )n1 dn where again dn is an integral over angles. Thus the final term in the last line of (10.49) is given by # # Gn (n ) dS = ) (n ) dn . (10.51) S!n1 Since is regular at y by assumption, the value of the remaining integral is bounded, so this term vanishes as ) 0. On the other hand, the penultimate term in the final line of (10.49) becomes # # 1 (n Gn ) dS = dn = (10.52) An S!n1 where is the average value of on the small sphere surrounding y. As the radius of this sphere shrinks to zero we have (y), the value of at the centre of the sphere. Putting all this together we find # # (y) = Gn 2 dV + (n Gn ) Gn (n ) dS (10.53) as our new form of Green's second identity involving Gn , where now we take the boundary of to be just the large sphere (the radius of the inner sphere having been shrunk to zero). - 135 - Finally, we take to be a solution to Poisson's equation 2 = F , whereupon we obtain Green's third identity # # (y) = Gn (x, y) F (x) dV + [(x) (n Gn (x, y)) Gn (x, y) (n (x))] dS (10.54) where the integrals are over the x variables. This is a remarkable formula. Recalling that y can be any point in Rn , we see that (10.54) describes the solution throughout our domain in terms of the solution on the boundary, the known function Gn and the forcing term. Also notice that (unlike the previous cases we have considered in the course) the Green's function is here providing an expression for the solution with inhomogeneous boundary conditions. Provided and n tend to zero suitably fast at asymptotically large distances, the boundary terms in (10.54) can be neglected. In this case, Green's third identity also shows that Gn is the free-space Green's function for the Poisson equation; this also follows formally from the definition 2 Gn = (n)(x y), since # # 2 2 (n) (x y) F (x) dV = F (y) (10.55) (y) = (y Gn (x, y)) F (x) dV = Rn Rn where the Laplacian is wrt y. There's something puzzling about equation (10.54) however. Setting F = 0, (10.54) provides an expression for the solution of the Laplace equation 2 = 0 in the interior of a domain in terms of the values of both and n on the boundary. But we already know that if is a compact domain with boundary , then | alone is enough to specify a unique solution to Laplace's equation (Dirichlet boundary condtions), while n| alone gives a solution unique up to a constant (Neumann). How then should we interpret the Green's identity (10.54), involving both | and n| ? The point is that whilst (10.54) (with F = 0) is a valid expression obeyed by any solution to Laplace's equation, it's not a constructive expression we can't use it to solve for in terms of boundary data, because the boundary values of and n cannot be specified independently. 10.3.3 Dirichlet Green's functions for the Laplace & Poisson equations In view of this curious 'over-determined' property of equation (10.54), we now finally ask how can the Dirichlet42 problem for the Laplace and Poisson equations be solved using Green's function techniques, in domains with boundaries? We do so by constructing a 42 It's also interesting to ask about the Neumann problem (where n is specified on the boundary), but in this course we'll discuss only the Dirichlet problem. The reason is that the Neumann problem is somewhat more complicated than the Dirichlet problem since (unlike the Dirichlet case) a consistency condition must be satisfied by the boundary data if a solution is to exist at all. Indeed, suppose we specify n | = g for some choice of function g : R. It turns out that we cannot pick g freely, since " " " " g dn1 S = n dn1 S = 2 dn1 x = F dn x , by the divergence theorem. Thus, in the Neumann case, the boundary data must be correlated with the integral of the forcing function F (the rhs of Poisson's equation), integrated over the whole domain. This consistency condition complicates the construction of a Green's function. - 136 - Dirichlet Green's function for the Laplacian operator on (with x, y ). This is defined to be the function G(x; y) := Gn (x; y) + H(x, y) where H is finite for all x (including at x = y), and where H satisfies the Laplace equation throughout . For G to be the Dirichlet Green's function, the boundary value of G should be chosen to cancel that of the fundamental solution Gn , so that G| = 0. Note that G also satisfies the Laplace equation whenever x .= y. We'll see a method for constructing such a Dirichlet Green's function in the following section (though in general it's a difficult problem finding H requires solving Laplace's equation with non-trivial boundary conditions). For now, let's just suppose G exists. If so, then given G we can use it to solve Poisson's equation on with Dirichlet boundary conditions for . To see this, suppose 2 = F inside , with | = g for some given function g : R. Substituting Gn = G H into Green's third identity (10.54) gives # # (y) = [ (n [G H]) [G H] (n )] dS F (x)[G H] dV . (10.56) Putting = H into Green's second identity (10.47) gives # # (n H) H (n ) dS = F H dV , (10.57) where we've used the fact that 2 H = 0 throughout , and that 2 = F . Comparing this with (10.56) shows that all the H terms in (10.56) cancel. We thus obtain # # (y) = (x) (n G(x, y)) G(x, y) (n (x)) dS F (x) G(x, y) dV . (10.58) Finally, having chosen H so as to ensure G| = 0 and using our boundary condition that | = g, we get # # (y) = g(x) nG(x; y) dS F (x) G(x; y) dV . (10.59) At last, this expression is constructive; the solution throughout is given purely in terms of the known boundary conditions and the Green's function. If instead we had a Neumann problem with n | = g being specified, then instead of H we'd seek a harmonic function H % to make n G| = 0. This condition defines H % only up to an additive constant. Then equation (10.58) gives us # # (y) = g(x) G(x; y) dS F (x) G(x, y) dV (10.60) as the solution to the Neumann problem. 10.4 The method of images For domains with sufficient symmetry we can sometimes use the elegant method of images (also called the reflection method) to construct the required Green's function. Indeed, this method can also be used for the Green's functions of the forced heat and wave equations - 137 - too, as well as for homogeneous equations. The key concept is to match the boundary conditions by placing a extending the domain beyond the region of interest, and placing a 'mirror' or 'image' source or forcing term in the unphysical region. We'll content ourselves with just illustrating how this works by looking at a few examples. (You can find more on the final problem sheet.) 10.4.1 Green's function for the Laplace's equation on the half-space Let = {(x, y, z) R3 : z 0} and suppose : R solves Laplace's equation 2 = 0 inside , subject to (x, y, 0) = h(x, y) and lim = 0 . |x| (10.61) Before we can use the formula of equation (10.59), we must construct a Green's function which vanishes on . As well as vanishing on the (x, y)-plane, we interpret the boundary condition on the asymptotic value of in the natural way by requiring G also vanishes as |x| . We'll set x = (x, y, z) and y := x+ 0 = (x0 , y0 , z0 ) in terms of Cartesian coordinates, with z0 > 0. We know that the free space Green's function G3 (x, x+ 0)= 1 1 4 |x x+ 0| satisfies all the conditions except that G3 (x; x0 )|z=0 = 1 1 .= 0 , 2 4 [(x x0 ) + (y y0 )2 + z02 ]1/2 (10.62) so that the homogeneous boundary condition G|z=0 = 0 does not hold for G3 (x; x+ 0 ). We thus need to 'cancel' the nonzero boundary value of G3 by adding on some function that solves Laplace's equation throughout our domain. Now let x 0 be the point (x0 , y0 , z0 ); in + other words, x is the reflection of x in the boundary plane z = 0. The location x / , 0 0 0 so the Green's function G3 (x, x0 ) is regular everywhere within , and so obeys Laplace's equation everywhere in the upper half-space. It's also clear from (10.62) that + G3 (x; x 0 )|z=0 = G3 (x; x0 )|z=0 . Thus the Dirichlet Green's function we seek is G(x; x0 ) := G3 (x; x+ 0 ) G3 (x; x0 ) 1 1 = . + + 4 |x x0 | 4 |x x 0| (10.63) Note that this expression obeys the asymptotic condition lim|x| G(x, x0 ) = 0 as well as the boundary condtion G(x; x0 )|z=0 = 0. The intuition behind this solution is as follows. Recall that the Green's function + + 2 (3) G3 (x; x+ 0 ) obeyed G3 (x; x0 ) = (x x0 ) and thus describes the contribution (x) - 138 - Figure 16. The contribution to (10.63) from the image source in the unphysical region z < 0 cancels that of the source in the physical region along the plane z = 0. You can see more examples at this Wikipedia page, from where I took the picture. receives from a point source of unit strength placed at x0 + . On the other hand, if we were solving the problem on R3 rather than the upper half-space, then G3 (x; x 0 ) would represent the contribution to (x) from a point source of equal magnitude but opposite sign placed at the mirror location x 0 . The combined effect of these two sources cancels out at any point x for which |x x+ 0 | = x x0 |, in other words along the plane z = 0 (see figure 16). This is equivalent to the boundary condition we wanted. To apply the formula (10.59), set F = 0 (we're considering Laplace's equation here) and note also that there is no contribution from the far field since 0 asymptotically by our boundary condition. The outward normal from the domain at z = 0 points in the negative z-direction, so the only contribution to (10.59) comes from the lower boundary: ( G (( (n G)|z=0 = z (z=0 ! " 1 z + z0 z z0 (10.64) = 3 3 4 |x x |x x+ 0| 0| z=0 53/2 z0 4 = (x x0 )2 + (y y0 )2 + z02 . 2 Therefore we get z0 (x0 , y0 , z0 ) = 2 # R2 4 h(x, y) (x x0 )2 + (y y0 )2 + z02 - 139 - 53/2 dx dy (10.65) as the solution of the Dirichlet problem for the Laplace equation in the half space. The method of images can also be applied to (some) more complicated domains. For example, suppose = {(x, y, z) R3 : a z b} is the region between two planar boundaries. We can view each of the two boundaries as mirrors and then any special point x0 gives rise to a double infinity of multiply reflected images not only must x0 itself be reflected in each mirror plane, but each of the image points obtained in this way must then be reflected again in the other mirror. If the mirror charges are taken with alternating plus and minus signs, we can construct a Green's function that vanishes on both boundaries. 10.4.2 Method of images for the heat equation We can also apply the method of images to the Green's functions that we developed for the forced wave and heat equations (with the initial conditions all being zero). For example, a chimney will provide a source of smoke F (x, t) from times t > 0 after the fire is lit. If this was a problem in R3 [0, ) then from equation (10.15) the density of smoke at some point x R3 at time t would be ! " # t# 1 |x y|2 ? (x, t) = F (y, ) exp d3 y d (10.66) 3/2 4D(t ) 3 (4D(t )) 0 R using the Green's function (t ) S3 (x y; t ) for the heat (or diffusion) equation on R3 [0, ). However, this is clearly incorrect in the present situation because the smoke does not diffuse into the ground (z < 0). The condition that there be no flow of smoke across the plane z = 0 amounts to the Neumann condition that z |z=0 = 0. We can ensure this is the case by modifying the Green's function to become 4 5 G(x, t; x0 , ) = (t ) S3 (x x+ 0 ; t ) + S3 (x x0 ; t ) (10.67) where x 0 = (x0 , y0 , z0 ). Again, the second term here obeys the homogeneous diffusion equation (t D2 )S3 (x x 0 ; t ) = 0 for all points x above ground. One can also easily check that the sum on the rhs of (10.67) indeed obeys the Neumann condition z G|z=0 = 0. Then the distribution of smoke at time t > 0 will be given by 4 5 # t# S3 (x x+ 0 ; t ) + S3 (x x0 ; t ) + (x, t) = F (x0 , ) d3 x0 d (10.68) 3/2 3 (4D(t )) 0 R in the region z > 0. We can imagine that the second term here describes the smoke emitted by a 'mirror chimney' the reflection of the real chimney in the boundary plane. The unmodified Green's function (for x R3 ) allows smoke to flow into the region z < 0, but this is exactly compensated by smoke flowing up into z > 0 from the mirror chimney. Note that we take the same source term F (x+ 0 , ) = F (x0 , ) in order for the smoke distributions to mirror eachother. I want to emphasize that the second source of smoke in the z < 0 region is (of course) entirely fictitious. It's just a device we've introduced in order to help solve our problem - 140 - on the upper half space with the condition that no smoke flows across z = 0 in either direction. We wanted to find the smoke distribution at points x in the upper half space. Our solution (10.68) is indeed valid in the region z > 0, and obeys the boundary condition |z=0 = 0 for all times. The uniqueness theorem for solutions of the heat equation on a domain with Dirichlet boundary conditions on then assures us that we've found the solution we wanted. We could also evaluate (10.68) at points with z < 0 and in fact, were this region to be filled with air rather than ground and were the mirror chimney to actually exist, our solution would be correct there. But this is clearly irrelevant to the physical situation. 10.4.3 Method of images for inhomogeneous problems For homogeneous diffusion or wave problems the method of images can be applied to the initial data itself, adapting our previous solutions to be applicable in the presence of extra boundaries. For example consider the wave equation for (x, t) on the half-line x 0, subject to the Neumann boundary condition that x (0, t) = 0, and with initial conditions (x, 0) = b(x) and and where b(x) = - t (x, 0) = 0 where x 0 , 1 x0 a x x0 + a 0 else (10.69) (10.70) is the 'top-hat' function. (Physically, the Neumann condition at x = 0 physically models small amplitude waves in the vicinity of a frictionless wall.) d'Alembert's general solution (9.45) shows that the solution on whole line x R would be a pair of boxes, each of height half and width 2a, one traveling to the left and one to the right with speed c. This solution clearly does not satisfy the Neumann boundary condition at x = 0 (when the left moving box passes over the origin). To overcome this, we consider an initial 'image box' on the negative x axis, obtained by reflecting the given initial conditions in x = 0. Taking this image box to have the same sign as the original box, d'Alembert gives us the solution (x, t) = b(x ct) + b(x + ct) b(x ct) + b(x + ct) + 2 2 (10.71) on the whole line. This describes two pairs of moving boxes, each member of each pair having height 1/2. The solution satisfies the initial condition (x, 0) = b(x) in the region x 0 and satisfies the Neumann boundary condition at x = 0 for all time. Thus, (x, t) restricted to the region x 0 is the solution to our original problem. It's instructive to picture how this solution evolves in time (especially in the region x 0). Let BR and BL denote the right and left moving boxes corresponding to the original initial box, and IR and IL the boxes from the initial 'image box'. IL just moves off to the left and never appears as part of the solution for x 0. BR just moves off to the right and is always present in the x 0 region. At time t = (x0 a)/c the leading - 141 - edges of both BL and IR reach x = 0. As they pass through each other (for a time 2a/c) the solution piles up to double height where they overlap and then BL goes off into x < 0 while IR emerges fully into x 0. Thereafter, the solution in x 0 is just two boxes of heights half, with centres separated by constant distance 2x0 (the initial separation of IL and BR ) moving off to the right at speed c. If we view the solution near x = 0 just on the positive side, we see BL \"hit the wall\\f