Solving Partial Differential Equations using Integral Transforms

Up to this point, we have employed the method of separation of variables to solve partial differential equations. Now let's look at applying integral transforms, specifically Foruier transforms and Laplace transforms. As usual, our goal in applying transform methods is to reduce the problem to one that is easier to solve and then to apply the inverse transform to the solution of the easier problem to obtain the desired solution to the original problem. When dealing with ODEs, we typically applied a transform to obtain an algebraic equation to solve. Here we will see that the transform will produce a problem where the number of independent variables has again been reduced by one, so that we obtain an ODE to solve for the transform of the desired solution.

To be concrete, let us consider a particular example problem of unsteady heat conduction in a thin, infinitely long rod that is insulated along its length. Assuming that there is no significant temperature variation across a transverse section of the rod, the temperature u(x,t) is governed by the one-dimensional heat equation:

(Subscript[u, t])=(c^2)(Subscript[u, xx]) , u(x,0) = f(x) , -∞<x<∞ , 0≤t<∞ .

Let's consider various approaches to solving this problem.

Separation of variables

Note that if we apply separation of variables and look for a solution of the form u(x,t)=X(x) T(t), everything would go smoothly and we would obtain the following separated equations:

T' + (c^2)(k^2)T = 0 , X'' + (k^2)X = 0

where (k^2) is the separation constant.These equations should look familiar and lead to solutions of the form

(Subscript[u, n])(x,t)= ((Subscript[A, n])Cos(kx) + (Subscript[B, n]) Sin(kx)) Exp(-(c^2)(k^2) t).

This is when we would normally apply boundary conditions to detemine allowed values of k. However, we do not have such boundary conditions to apply in this problem. Instead of a discrete spectrum of allowed frequencies (or periods) for the eigenfunctions, we now have a continuous spectrum, and we should expect our solution to be in the form of an integral rather than a summation. This is exactly the same situation we ran into when we looked at ODEs with periodic inputs (described by Fourier series) and let the period become unbounded to arrive at the Fourier integral and Fourier transform. (Note that there are also finite Fourier transforms for handling problems on a finite domain, but we will not need them here.)

Fourier Transform

Let's take that as a hint to now apply the Fourier transform to our PDE. The first issue to deal with is to determine which of the two independent variables we should apply the transform to. Recall that the Fourier transform involves a doubly infinite domain (i.e., -∞<x<∞), so let's apply the Fourier transform to the spatial variable. Recall the Fourier transform pair:

ℱ{u(x,t)} =(1/(Sqrt[2π]))(∫_(-∞)^∞)u(x,t)(e^(-iωx))dx=(OverHat[u])(ω,t) 

(OverHat[u])(ω;t)=(1/(Sqrt[2π]))(∫_(-∞)^∞)(OverHat[u])(ω,t)(e^(iωx))dω

and the operational property of how the transform works on derivatives with respect to the variable being transformed. Each derivative produces a factor of iω, so

ℱ{(Subscript[u, xx])(x,t)} =-(ω^2) (OverHat[u])(ω,t).

When transforming the time derivative, the derivative can be pulled outside the integral:

ℱ{(Subscript[u, t])(x,t)} =(1/(Sqrt[2π]))(∫_(-∞)^∞)(Subscript[u, t])(x,t)(e^(-iωx))dx= (∂/(∂t))  (OverHat[u])(ω,t).

Applying the Fourier transform to the PDE then produces an ODE for the Fourier transform of the solution with respect to the spatial variable:

(∂/(∂t))  (OverHat[u])(ω,t)+(ω^2)(c^2) (OverHat[u])(ω,t) = 0.

If we just consider ω as a parameter, this is a familiar (linear, first order, constant coefficient) ODE that we know how to solve:

(OverHat[u])(ω,t) = C(ω) (e^(-(c^2)(ω^2)t)).

Evaluating at t=0, C(ω) = (OverHat[u])(ω,0) = (OverHat[f])(ω);  i.e. the arbitrary "constant" in the solution (which is really a function of the parameter ω) corresponds to the Fourier transform of the initial temperature distribution.

We can now write a formal solution by applying the inverse Fourier transform:

u(x,t) = (1/(Sqrt[2π]))(∫_(-∞)^∞)(OverHat[u])(ω,t) (e^(iωx)) dω= (1/(Sqrt[2π]))(∫_(-∞)^∞)(OverHat[f])(ω)  (e^(-(c^2)(ω^2)t))  (e^(iωx)) dω.

Concrete example: f(x) = δ(x)

To make things more specific (and so that we can plot our results), let's consider a specific initial temperature distribution. In particular, let's put a "unit of heat" at the origin so that the initial temperature distribution is a delta function. (This may not seem like an intuitive choice, but we will see that it is useful.)

Recalling how the δ function operates inside an integral (it extracts the value of the remaining part of the integrand where the argument of the δ function is zero) the transform of this initial distribution evaluates to:

(Subscript[(OverHat[f]), δ])(ω) = (1/(Sqrt[2π]))(∫_(-∞)^∞)δ(x) (e^(-iωx))dx  = (1/(Sqrt[2π]))

and the transform of the corresponding solution is:

(Subscript[(OverHat[u]), δ])(ω,t) = (1/(Sqrt[2π])) (e^(-(c^2)(ω^2)t))

and a function of this form is conveniently located as formula #9 in the table of Fourier transforms on p. 578 in the text which states that:

ℱ{(e^(-(ax^2)))} =(1/(Sqrt[2a]))(e^(-(ω^2)/4a))   provided that  a>0.

So if we identify (c^2)t=1/4a (and, accordingly, (1/(Sqrt[2a]))=(Sqrt[2(c^2)t])), we now have

(Subscript[(OverHat[u]), δ])(ω,t) = (1/(Sqrt[4π (c^2)t])) (e^(-(ω^2)/4(c^2)t))= (1/(Sqrt[4π (c^2)t]))ℱ{(e^(-(1/4(c^2)t)(x^2)))}.

So the solution for the how the temperature evolves over space and time is given by:

(Subscript[u, δ])(x,t) = (1/(Sqrt[4π (c^2)t]))(e^(-(x^2)/(4(c^2)t))).

Visualize the results.

We can now make some plots so show how the temperature distribution evolves. Let's begin by creating a function corresponding to the solution and plotting the spatial distribution at particular values of time.

solution[x_, t_, c_] := Exp[-x^2/(4*c^2*t)]/Sqrt[4*Pi*c^2*t]

Note that plotting the function exactly at t=0 will be problematic, so let's start at some small positive time (t=0.00025) and look at the distribution at intervals of 1/2 time unit.

Table[
Plot[solution[x,0.00025+t, 1], {x,-5,5}, PlotRange -> {0,1},
Frame -> True,
FrameLabel -> {"X","u(x,t)","Time = "<> ToString[t],""},
PlotStyle -> AbsoluteThickness[2]], 
{t,0,4,0.5}]

Show[GraphicsArray[Partition[%,2]]]

The initial impulse spreads and flattens over time until, asymptotically, the temperature approaches zero everywhere. Below, we see that if we decrease c (i.e., decrease the conductivity), the spreading/flattening takes longer.

Table[
Plot[solution[x,0.00025+t, 0.5], {x,-5,5}, PlotRange -> {0,1},
Frame -> True,
FrameLabel -> {"X","u(x,t)","Time = "<> ToString[t],""},
PlotStyle -> AbsoluteThickness[2]], 
{t,0,4,0.5}]

Show[GraphicsArray[Partition[%,2]]]

Alternatively, we can construct 3-dimensional plots of temperature as a function of time and space:

Table

Table[
Plot3D[solution[x,t,c], {x,-5,5},{t,.000025, 4},
PlotRange -> {0,1}, AxesLabel -> {"X", "Time", "u[x,t,c]"},
PlotLabel -> "c= "<>ToString[c]],
{c,1, 0.5,-0.5}]

More realistic initial distributions

The solution we have obtained that describes the diffusion of a "heat impulse" along a rod is referred to as the 1-D diffusion kernel. While we have not shown it yet, the diffusion in a 2-D plate or a 3-D solid looks very similar:

u(x,t) = (1/((((Sqrt[4(c^2) π t])))^D))(e^(-(r^2)/(4(c^2)t))) = (Subscript[ker, D])(x,t)

where D corresponds to the spatial dimension of the domain. We can use these diffusion kernels to express the solution from a more realistic initial temperature distribution using convolution. If the initial temperature distribution in a rod is u(x,0)=f(x), then the temperature evolves according to the solution:

u(x,t) = (∫_(-∞)^∞)f(ξ) (Subscript[ker, 1])(x-ξ,t)dξ = (∫_(-∞)^∞)f(ξ) ((1 )/(Sqrt[4(c^2) π t]))(e^(-(((x-ξ))^2)/(4(c^2)t)))dξ  .

Again, let's be particular and consider the case of a semi-infinite rod at temperature zero (i.e., u(x,t)=0 , -∞<x<0) that is brought into contact with another semi-infinite rod at temperature (i.e., u(x,t)=1 , 0<x) at t=0. At this point, we can obtain the solution directly by plugging in a unit step function for f(x) which effectively just changes the integration interval to o<x<∞. We can clean up the solution a bit by changing variables to

y=((x-ξ)/(Sqrt[4(c^2)t]))

and the solution becomes

u(x,t)=((-1)/(Sqrt[π]))(∫_(x/(Sqrt[4(c^2)t]))^(-∞))  (e^(-(y^2)))dy.

After a few more manipulations (use the minus sign to exchange the limits and break the integral into two pieces separated at x=0), we obtain

u(x,t) = (1/2)(2/(Sqrt[π]))(∫_(-∞)^0)(e^(-(y^2)))dy+(1/2)(2/(Sqrt[π]))(∫_0^(x/(Sqrt[4(c^2)t])))(e^(-(y^2)))dy.

The first integral evaluates to ((Sqrt[2])/π)by a handy trick (multiply it by another version of itself with x as the integration variable and you find that the square of the integral is the area of the positive quadrant of the unit circle). The second integral is a tabulated function known as the "error function" (abbreviated "Erf") because the integrand is a Gaussian, so the integral corresponds to the cumulative distribution function for the Gaussian or normal distribution. So, at long last, we can write the solution compactly as

u(x,t) = (1/2)+(1/2) Erf((x/(Sqrt[4(c^2)t]))).

Now let's create a function corresponding to our solution associated with unit step initial conditions and visualize how the temperature evolves over time. Again, let's start with a small positive value of time, but here let's use successively longer time increments:

(Subscript[u, step])[x_,t_, c_] := 1/2*(1+Erf[x/Sqrt[4*c^2*t] ])

Table[
Plot[(Subscript[u, step])[x,10^n, 2], {x,-10,10}, PlotRange -> {0,1}, AxesLabel -> {"X", "u[x,"<>ToString[N[10^n]]<>",2]"},
PlotStyle -> AbsoluteThickness[2]],
{n,-3,4}]

Show[ GraphicsArray[Partition[%,2]]]

Note that this solution does not suffer from the high frequency oscillations and overshoot (Gibb's phenomenon) that we observed when using series solutions  to deal with discontinuous initial conditions.

Laplace Transform

Now let's consider approaching the same problem using the Laplace transform. Recall that the Laplace transform operates on a semi-infinite interval, so we can aplly the Laplace transform to the time variable in our problem.

Again, we employ the Laplace transform as usual:

ℒ{u(x,t)} = (∫_(-∞)^∞)u(x,t)(e^(-st))dt = g(x,s)

with the usual operational property for time derivatives

ℒ{(Subscript[u, t])(x,t)} = sℒ{u(x,t)} - u(x,0) = s g(x,s) - u(x,0).

Now the transform is defined by an integral over the time variable, so the spatial derivatives can be pulled outside the integral:

ℒ{(Subscript[u, xx])(x,t)}=((d^2)/(dx^2))ℒ{u(x,t)}= ((d^2)/(dx^2))g(x,s).

Applying the Laplace transform to the PDE ((Subscript[u, t])=(c^2)(Subscript[u, xx])) and initial condition (u(x,0)=δ(x)) now produces a different ODE:

s g(x,s)-δ(x)=(c^2)((d^2)/(dx^2))g(x,s)

or in more traditional form:

((d^2)/(dx^2))g(x,s)-(s/(c^2))g(x,s) = ((-1)/(c^2))δ(x).

The complementary solution is straightforward:

(Subscript[g, c])(x,s) = A(s) (e^(-(x/c)(Sqrt[s])))+B(s)(e^(+(x/c)(Sqrt[s])))

and the particular solution can be found by variation of parameters. (Consider A and B as functions of x, plug into the ODE and require - in addition to satisfying the inhomogeneous ODE - that the second derivative terms vanish so the new equations for A and B are first order.) The corresponding "variable parameters" are:

A(x,s)={({{(1/(2c(Sqrt[s]))), x>0}, {0               , x<0}})

B(x,s)={({{0                 , x>0}, {((-1)/(2c(Sqrt[s])))  , x<0}})

Putting the pieces together gives the final version of the transform

g(x,s) =(1/(2c(Sqrt[s])))exp(-(Sqrt[s])|x|/c)

which canalmost be found in the table of Laplace transforms (actually it is the derivative of formula #39 in the table on p. 298 in the text with k=|x|/c )

ℒ{(e^(-k(Sqrt[s])))}= (k/(2(Sqrt[π (t^3)])))(e^(-(k^2)/4t))

Together with the rule for inverting derivatives of known transforms:

(ℒ^(-1)){F'(s)}=-t f(t)

we again obtain the solution:

(Subscript[u, δ])(x,t) = (1/(2c(Sqrt[π t])))(e^(-(x^2)/(4(c^2)t))).

Recap

So while we have focused on using separation of variables for solving linear partial differential equations, today's example has (hopefully) demonstrated the transform methods can also be applied. The Fourier transform can be applied to a doubly infinite domain (such as the spatial domain in the example) and the Laplace transform can be applied on a singly infinite domain (such as the time domain in the example). There are also numerous other transform methods that can be applied to various classes of problems: finite Fourier transforms for handling problems on a finite domain, Hankel transforms for problems in polar coordinates, and many others.

You may or may not run across these in the future but, in case you do, the same concepts developed here apply. The transform is applied in an attempt to obtain an easier problem (typically a PDE with one less independent variable or an ODE) whose solution can be inverted to obtain the solution of the original problem. The inverse transform can be accomplished using tables or by computing the integral definition of the inverse transform, and (time permitting) we may get to do an example or two of computing the inversion integral at the end of this class.