Chapter 7
Transform Methods and Green's Functions

7.1 The Fourier Transform

Fourier Integral Formula (Figure 7.1)

The convergence of the limit in Equation (7.6) can be examined within Mathematica by constructing the Fourier series representation in (7.1) over the interval (-L/2,L/2) for increasing values of L.  We consider the function

                                               f[x]=Exp[-x^2]

which clearly satisfies the integrability condition in Equation (7.7).  The assembly of
the various terms in (7.6) proceeds in very much the same way as the calculations
in Example 4.2.  The Fourier coefficients in (7.3) (written as a[n,L] below) are determined and substituted into Equation (7.6), which is represented as a partial sum:

a[n_,L_]=Integrate[Exp[-t^2]Exp[I(2 n Pi t)/L],
     {t,-L/2,L/2}]

Do[A[k,1]=N[a[k,1]],{k,-20,20}]

f[x_,L_]=Sum[(1/L)A[k,L]Exp[-I(2 k Pi x)/L],
     {k,-20,20}];

Plot[f[x,1],{x,-2,2},PlotRange->{0,1},
     AxesLabel->{x,f}]

Do[A[k,2]=N[a[k,2]],{k,-20,20}]

Plot[f[x,2],{x,-2,2},PlotRange->{0,1},
     AxesLabel->{x,f}]

Do[A[k,10]=N[a[k,10]],{k,-20,20}]

Plot[f[x,10],{x,-2,2},PlotRange->{0,1},
     AxesLabel->{x,f}]

We have retained both the cosine and sine contributions in these expressions
despite f[x] being an even function.  The quantity n in f[x,n,L] must be large enough to insure convergence to the desired threshold.  We have used n=20 for the calculations performed here. 

The Fourier representations for L=1, L=2 and L=10 are shown in Figure 7.1.  As L increases, successively larger domains of the function are included in the periodic unit of the Fourier representation and the limit in (7.6) is clearly seen to converge to the original function.

Example 7.1

Fourier transforms can be computed with Mathematica both symbolically and numerically.  We consider as an example the function f[x]=0 if |x|>1 and f[x]=0 otherwise.  The Fourier transform of this function is computed symbolically with the following entry:

(Integrate[Cos[k x],x]/.x->1)-
      (Integrate[Cos[k x],x]/.x->-1)

The inverse transform is carried out with the command

F[k_]:=If[k!=0,(2 Sin[k])/k,2]

f[t_,n_]:=NIntegrate[F[k]Cos[k t]/(2Pi),
     {k,-n Pi,n Pi}]

This second of these commands is seen to generate the equivalent in Mathematica of the Fourier integral theorem in Equation (7.9) as the range of integation is increased to span the entire real line. The plots of the inverse transform integrated over an increasing region of k is shown in Figure 7.2.  The gross features of the function in (7.14) are obtained from the Fourier components with small values of k, which correspond to long wavelength variations of the function.  However, reproducing the discontinuities at x=1 and x=-1 requires the higher Fourier components to be included in the range of integration.  This is similar to the behavior found for Fourier series representations of functions with discontinuities.

f[x_]:=If[Abs[x]<=1,1,0]

Plot[{f[x,1],f[x]},{x,-2,2},PlotRange->{-0.1,1.2},
     AxesLabel->{x,f}]

Plot[{f[x,2],f[x]},{x,-2,2},PlotRange->{-0.1,1.2},
     AxesLabel->{x,f}]

Plot[{f[x,3],f[x]},{x,-2,2},PlotRange->{-0.1,1.2},
     AxesLabel->{x,f}]

A different type of behavior is seen for smooth functions.  Consider again the
function in (7.10).  The Fourier transform can be determined by similar commands
as those used above and the rate of convergence is displayed in Figure 7.3.  The
smooth form of (7.10) is recovered with a much smaller range of integration than
that used above, as expected from the rapid decay of the Fourier transform with
increasing k.

f[x_]:=Exp[-x^2]

F[k_]=(Integrate[f[x]Exp[I k x],x]/.x->Infinity)-
      (Integrate[f[x]Exp[I k x],x]/.x->-Infinity)

f[t_,n_]:=NIntegrate[F[k]Cos[k t]/(2Pi),
     {k,-n Pi,n Pi}]

Plot[{f[x,1/4],f[x]},{x,-2,2},PlotRange->{-0.1,1},
     AxesLabel->{x,f}]

Plot[{f[x,1/2],f[x]},{x,-2,2},PlotRange->{-0.1,1},
     AxesLabel->{x,f}]

Plot[{f[x,3/4],f[x]},{x,-2,2},PlotRange->{-0.1,1},
     AxesLabel->{x,f}]

7.2 The Laplace Transform and the Bromwich Integral

Example 7.1

The package LaplaceTransform.m (which may be found in the file called 'Supplementary') can be used to perform Laplace and inverse Laplace transforms both for particular functions and symbolically.  For example, the Laplace transform of e^[at] and the corresponding inverse transform are determined by the following commands:

<<LaplaceTransform.m

LaplaceTransform[Exp[a t],t,s]

InverseLaplaceTransform[%,s,t]

These packages can also be used to perform symbolic Laplace transforms as in 
Equations (7.29) and (7.30) and applied to solving differential equations such as (7.28) with general initial conditions (7.29).  The procedure is the same as in Equations  (7.30)-(7.34).  The Laplace transform Y[s] of y[t] is first defined as a generic function and then LaplaceTransform is applied to the differential equation (7.28):

LaplaceTransform[y[x_],x_,s_]=Y[s]

LaplaceTransform[#,x,s]& /@ (y''[x]+2y'[x]+2y[x]==1)

The Laplace-transformed equation is then solved for Y[s] as in Equation (7.32), and InverseLaplaceTransform is used to obtain the solution of the differential equation:

Y[s] /. Solve[%,Y[s]][[1]]

InverseLaplaceTransform[%,s,x] // Expand

which agrees with the solution obtained in (7.34). Additional examples are provided in the problems at the end of the chapter.

7.3 Fundamental Solution of the Diffusion Equation

Figure 7.5

G0[x_,x0_,t_]:=(1/Sqrt[4 Pi t])Exp[-(x-x0)^2/(4 t)]

Plot[{G0[x,0.5,10^(-10)],G0[x,0.5,0.005],
     G0[x,0.5,0.01],G0[x,0.5,0.1]},{x,0,1},
     PlotRange->{0,5},AxesLabel->{x,G}]

Example 7.3 (Figure 7.6)

The solution (7.51) of the initial-value problem of the diffusion equation can be constructed by carrying out the corresponding operations in Mathematica.  Consider the initial condition where f[x]=1, if |x| is less than or equal to 1, and f[x]=0 if |x| is greater than 1. The solution in (7.51) is obtained by constructing this initial function and the fundamental solution (7.46) and then performing the
convolution of these two functions using NIntegrate:

f[x_]:=If[Abs[x]<=1,1,0]

G0[x_,t_]:=(1/Sqrt[4 Pi t])Exp[-(x^2)/(4 t)]

u[x_,t_]:=NIntegrate[f[y] G0[x-y,t],{y,-1,1}]

The early stages of the time development of this solution is shown in Figure 7.6.
The sharpest features in the profile of f (near the discontinuities) are smoothed out first as the profile spreads out as it decays to zero.  This is to be expected from the behavior of the fundamental solution in Figure 7.5 and the time development of the individual Fourier components in (7.39).

Figure 7.6 is generated with Plot as follows:

Plot[{f[s],u[s,0.05],u[s,1]},{s,-3,3},
     PlotRange->{0,1},AxesLabel->{x,u}]

7.4 Fundamental Solution of Poisson's Equation

The integral in (7.64) can be evaluated by using Mathematica in two ways.  The integral can be evaluated symbolically by using Integrate:

Integrate[Sin[x]/x,{x,0,Infinity}]

Alternatively, the numerical evaluation of (7.66) proceeds by the following command sequence.  Notice that this procedure is different than that described in the book.  The commands shown below are much more efficient, since the each range of the integral is calculated only once.  The effects of this difference are more evident in version 2.1 than they were in version 1.2.

PartialIntegral[i_]:=NIntegrate[Sin[t]/t,
     {t,(i-1)^4,i^4},MaxRecursion->20,
     WorkingPrecision->20,AccuracyGoal->7]

OldValue=0;
     Do[
          NewValue=OldValue+PartialIntegral[i];
          OldValue=NewValue;
          Print[N[NewValue,7]],
     {i,1,15}]

To four decimal places the numerical value of this integral agrees with the exact result: Pi/2=1.570796.  Listing the value of the integral with an increasing upper bound shows the slow convergence, which is due to the oscillatory nature of the integrand.

7.6 Green's Functions in the Presence of Boundaries:
      The Method of Images

Figure 7.7

The two panels of this figure are generated by first constructing the fundamental solution and then taking the appropriate linear combination to satisfy the boundary condition at the origin.  The fundamental solution is entered as

G0[x_,x0_,t_]:=(1/Sqrt[4 Pi t])Exp[-(x-x0)^2/(4 t)]

The Green's function for an absorbing boundary at the origin is then given by

GA[x_,x0_,t_]:=G0[x,x0,t]-G0[x,-x0,t]

The top panel of Figure 7.7 is obtain as follows:

Plot[{GA[x,0.5,0.005],GA[x,0.5,0.01],GA[x,0.5,0.05],
     GA[x,0.5,0.15]},{x,0,1},PlotRange->{0,4},
     AxesLabel->{x,G},Ticks->{Automatic,Range[0,4,1]}]

The Greens' function for a reflecting boundary at the origin is given by

GR[x_,x0_,t_]:=G0[x,x0,t]+G0[x,-x0,t]

and the bottom panel of Figure 7.7 is constructed in a similar fashion as the top panel:

Plot[{GR[x,0.5,0.005],GR[x,0.5,0.01],GR[x,0.5,0.05],
     GR[x,0.5,0.15]},{x,0,1},PlotRange->{0,4},
     AxesLabel->{x,G},Ticks->{Automatic,Range[0,4,1]}]

Figure 7.8

The three panels of this figure are generated by following steps analogous to used to obtain the two panels of Figure 7.7.  The fundamental solution is first constructed and the appropriate linear combinations are taken in order to satisfy the boundary conditions at x=0 and x=1.  Only the nearest image solutions are needed for the times considered

G0[x_,x0_,t_]:=(1/Sqrt[4 Pi t])Exp[-(x-x0)^2/(4 t)]

Top panel:

GAA[x_,x0_,t_]:=G0[x,x0,t]-G0[x,-x0,t]-G0[x,2-x0,t]

Plot[{GAA[x,0.5,0.005],GAA[x,0.5,0.01],
     GAA[x,0.5,0.05],GAA[x,0.5,0.15]},{x,0,1},
     PlotRange->{0,4},Ticks->{Automatic,Range[0,4,1]},
     AxesLabel->{x,G}]

Middle panel:

GAR[x_,x0_,t_]:=G0[x,x0,t]-G0[x,-x0,t]+G0[x,2-x0,t]

Plot[{GAR[x,0.5,0.005],GAR[x,0.5,0.01],
     GAR[x,0.5,0.05],GAR[x,0.5,0.15]},{x,0,1},
     PlotRange->{0,4},Ticks->{Automatic,Range[0,4,1]},
     AxesLabel->{x,G}]

Bottom panel:

GRR[x_,x0_,t_]:=G0[x,x0,t]+G0[x,-x0,t]+G0[x,2-x0,t]

Plot[{GRR[x,0.5,0.005],GRR[x,0.5,0.01],
     GRR[x,0.5,0.05],GRR[x,0.5,0.15]},{x,0,1},
     PlotRange->{0,4},Ticks->{Automatic,Range[0,4,1]},
     AxesLabel->{x,G}]