Linear Methods of Applied Mathematics

Evans M. Harrell II and James V. Herod*

*(c) Copyright 1994,1995,1996, 2000 by Evans M. Harrell II and James V. Herod. All rights reserved.

version of 16 September 2000, text with


If you wish to print a nicely formatted version of this chapter, you may download the rtf file, which will be interpreted and opened by Microsoft Word or the pdf file, which will be interpreted and opened by Adobe Acrobat Reader.

(Some remarks for the instructor).

Some of the calculations of this chapter are available in a Maple worksheet or in a Mathematica notebook.

X. PDEs on a disk

Let us look again at the potential equation of Laplace, but this time in polar coordinates. This equation would apply to the equilibrium temperature distribution in a thin disk insulated on its two flat faces, or to a long cylinder along the z-axis, provided that the temperature does not depend on z. A temperature could be prescribed on the outside edge as a function of the angular variable theta producing a boundary-value problem of the form

   Grad^2 u = partial ^2 u / partial r^2 + (1/r) partial  u/partial r + (1/r^2) partial ^2 u/partial theta^2 = 0.          (10.1)
    u(a, theta) = f(theta)          (10.2)

The PDE (10.1) is what you get if you change variables using the 2-variable chain rule. I shall leave this to the exercises, but notice at least that (10.1) is dimensionally consistent: r has units of length and partial / partial r of length-1, whereas theta is dimensionless. Hence each term in (10.1) has dimensions length-2.

If we try to find a product solution for (10.1) of the form

u(r, theta) = R(r) Q(theta), we readily find

   0 = R''[r]/R[r] + (1/r) R'[r] /R[r] + Q''[theta]/ Q[theta]

The equation for Q is the familiar one that gives sines and cosines, but what are the boundary conditions?
Answer: periodic BC. Since there is no physical difference between theta and theta + 2 pi, we must have:

   Q(theta + 2 pi) = Q(theta) (PBC)

In order to arrange periodicity, we take solutions in the form

where m is an integer, or, equivalently, Q(theta,m) = c(m) exp(i m theta) + c(-m) exp(-i m theta).              &# 160;           (10.4) The function R satisfies the equidimensional, or Euler, equation:

   R'' + (1/r) R' - (m2/r2) R = 0 (10.5)

Equidimensional equations can be solved by making the guess R = ra and figuring out what the constant a has to be. In this case, a = m, forced by the periodicity of the solutions (10.3), so the general solution is of the form:

   a_0 + Sum for m=1 to infinity of r^m (a_m cos[m theta] + b_m sin[m theta])             + Sum for m=1 to infinity of r^{-m} (a_{-m} cos[m theta] + b_{-m} sin[m theta])           (10.6)

If we are studying the equilibrium temperature (or potential, or equilibrium displacement of a membrane) on a disk, the solutions containing r-m, m>0, are unphysical, so in this case the general solution will be of the form:

   a_0 + Sum for m=1 to infinity of r^m (a_m cos[m theta] + b_m sin[m theta])           (10.7)

where the sum now runs from m=1 to infinity. I have written a0 separately for the same reasons as with Fourier series.

Model Problem X.1. Solve (10.1) with the boundary condition    u(1, theta) = (cos(theta))2.


To find a solution in the interior, we need to determine the coefficients in (10.7). This is done by substituting r = 1 and comparing:

In other words, the coefficients a[m] and b[m] are none other than the usual Fourier coefficients for the function f, here (cos(theta))2. Rather than doing integrals here, the easiest way to find these Fourier coefficients is to use the double-angle formula in the form:

    cos(2 theta) = 2 cos2(theta) - 1

Remember - if you can find Fourier coefficients by some trick, the answer is perfectly correct, since the Fourier series is uniquely determined. This is great thing about uniqueness theorems in mathematics - they allow you find an answer by the easiest method available, without worrying that it might be a different answer from the one you would get from a harder but more standard procedure.

We conclude that    (cos[theta])^2 = 1/2 + (1/2) cos[2 theta]

or in other words that a0 = 1/2, a2 = 1/2, and all other coefficients are 0.

That is the solution on the boundary of the disk (r=1); inside the disk the temperature is:

   (cos[theta])^2 = (1/2) + r^2 cos[2 theta].

Next, let us look again at a time-dependent problem, such as the heat equation on a disk, in polar coordinates. (There is a Mathematica notebook parallel to this text, but emphasizing the wave equation rather than the heat equation.)

Setting k=1 for convenience, the heat equation in takes on the form:

   u_t = partial ^2 u / partial r^2 + (1/r) partial  u / partial r + (1/r^2) partial ^2 u / partial theta^2.           (10.8)

As before, we begin by searching for product solutions. Recall from Chapter IX that it is convenient to separate vector variables one at a time in the product solution. If

   u = T(t) V(r, theta),

then we find that

    T'/T = (Grad2 V)/V          (10.9)

and we can make the usual argument that the left side is independent of space while the right side is independent of time, hence both equal some constant, called - lambda. The solutions for T are of the form

   T(t) = const. Exp(-lambdat),

where lambda should be determined from the boundary conditions in space.

Model Problem X.2. Suppose that a thin disk of radius A is insulated on its faces while the round edge is held at temperature 0. Find the eigenfunctions and eigenvalues for the spatial part of the separated equation (10.9).


On the edge we have Dirichlet boundary conditions of the form    u(t,A, theta) = 0           (DBC)

This two-dimensional problem, resembles (10.1) except that in the eigenvalue equation,

   - Grad^2 V = lambda V,

the eigenvalue lambda need not be zero. The constant lambda may be zero, in which case the disk is at equilibrium and because of the boundary conditions at r=A, the equilibrium temperature is 0 throughout. This trivial solution is not considered an eigenfunction, since it is not useful in building a series for a general solution.

Separating variables again by writing

   V = R(r) Q(theta),

leads to the same eigenvalue equation for Q(theta) as before. The boundary conditions are also the same, periodic, so the solutions are the same sines and cosines as in (10.3), indexed by the
integer m.

Something new occurs for the other part of the solution, R. Evaluating

   - Grad^2 V/V

with the knowledge that Q(theta) = sin(2 m theta) or cos(2 m theta), we get an ordinary differential equation for R. In place of the Euler, equation 10.5), which was solved with power functions, we get Bessel's equation,

   R'' + (1/r) R' + (lambda- m2/r2) R = 0,              (10.10)

which no longer has elementary solutions. Just as with the equidimensional equation, near r=0 the solutions can behave like either rm or r-m, but the latter are unphysical because the temperature can't diverge at r=0. The regular solutions, behaving roughly like rm near r=0, are called regular Bessel functions and, with a standard choice of the overall constant, denoted Jm(lambda1/2 r). Mathematica uses the notation BesselJ. They are normalized so that Jmx) ~ xm/(2 m!) as x -> 0:


plot of J_0 Bessel function
plot of J_1 Bessel function

The oscillations are damped, and not quite as evenly spaced as those of the sines and cosines, but, as you might expect, there are many similarities between these functions and the sines and cosines. With the aid of a computer and software, Bessel functions should not be considered much more exotic or hard to use than sines and cosines.

In order to satisfy the boundary condition at A, we must fix lambdaso that one of these nodes happens at r=A. It turns out that the nodes of the Bessel functions are tabulated functions, and they are traditionally denoted

   jm,n = the n-th positive value of x for which Jm[x] = 0.

They are also quite easy for Mathematica to compute. There is a Mathematica package which tabulates them, called BesselZeros.m, but let's proceed on our own:. Looking at the graph, we choose reasonable guesses for these numbers, 2,5,9,15, etc., and let the FindRoot command locate the nearest actual root to our guess:

In:= FindRoot[{BesselJ[0,j[0,1]] == 0,BesselJ[0,j[0,2]] == 0, \
    BesselJ[0,j[0,3]] == 0,BesselJ[0,j[0,4]] == 0}, \
    {j[0,1],2}, {j[0,2],5}, {j[0,3],9},{j[0,4],15} ]

  {j[0, 1] -> 2.40483, j[0, 2] -> 5.52008, j[0, 3] -> 8.65373, 
    j[0, 4] -> 14.9309}

  FindRoot[{BesselJ[1,j[1,1]] == 0,BesselJ[1,j[1,2]] == 0, \
    BesselJ[1,j[1,3]] == 0,BesselJ[1,j[1,4]] == 0}, \
    {j[1,1],4}, {j[1,2],10}, {j[1,3],10},{j[1,4],20} ]

  {j[1, 1] -> 3.83171, j[1, 2] -> 10.1735, j[1, 3] -> 10.1735, 
    j[1, 4] -> 19.6159}
And so forth. The important point is that today Bessel functions and their zeroes jmn should be regarded as known, easily accessible numbers, nearly as accessible as pi and values of the sine.

The eigenfunctions and eigenvalues for Model Problem X.2 are (as usual, up to a constant, nonzero multiple):

J_m(j_{mn} r/A), with lambda{mn} = (j_{mn}/A)^2

Before continuing with the heat equation on the disk, a few more remarks about Bessel functions are in order.

One of the ways in which Bessel functions resemble sines and cosines is that they are the source of some complete sets, with which we can define Fourier-Bessel series. This is not an accident, of course; the eigenfunctions of any eigenvalue problem which satisfies the technical condition of self-adjointness form a complete set by the spectral theorem of J. von Neumann. Moreover, the eigenfunctions are orthogonal. (If the eigenvalues are different this is automatically true, and even if there are several eigenfunctions with the same eigenvalue, they can be recombined as an orthogonal set.)

Here is how it works for Bessel functions:

Just as we have different Fourier series, which use sines, cosines, or both, there are different Fourier-Bessel series. Each series uses the J's with a fixed index (like m) in place of the expressions n pix/L or 2n pix/L which appear in Fourier series, there is an expression involving jm,n.


   Integrate[BesselJ[m, r  j[m,n]/A] BesselJ[m, r  j[m,p]/A] r, {r,0,A}] =(A^2 /2) (BesselJ[m+1, j[m,n]])^2       delta[n,p].           (10.11)

You may at first think that the factor of r in this integral is peculiar, but it is really inherited from the two-dimensional area element in polar coordinates, which is

   d Area = r dr dtheta.

The alert reader will perhaps recall that when inner products were introduced we learned that the integrals could incorporate weights like the extra r here.

The product solutions for V(r, theta) found above can be written either with sines and cosines, or with exponential functions, and I shall make the latter choice:

   {BesselJ[Abs[m], r  j[m,n]/A]  Exp[I m theta] },

are orthogonal in the sense of the two-dimensional inner product

   <f,g>, with an integral r dr d theta.

Here, m can be positive, negative, or 0. With the normalization, the orthogonality relationship is:    <BesselJ[Abs[m], r j[m,n]/A] Exp[I m theta], BesselJ[Abs[q], r j[q,theta]/A] Exp[-I q
theta]> = Delta[m,p] Delta[n,theta] (A^2 Pi ) (BesselJ[m+1, j[m,n]])^2.


Note that if we were to consider the wave equation rather than the heat equation, with Dirichlet boundary conditions on a circle, the spatial dependence of the normal modes would be the same as the eigenfunctions of the Laplace operator as found in Model Problem X.2. Several of these normal modes may be viewed in a Mathematica notebook.

The following Fourier-Bessel series are useful, and result from identities involving the derivatives of Bessel functions:

   r^p = 2 A^p Sum[BesselJ[p, r j[p,n]/A] /(j[p,n] BesselJ[p+1, j[p,n]]),  {n,1,}]

(valid for each fixed m, 0 <= r < A).

   1 - (r/A)^2 = 8 Sum[BesselJ[0, r j[m,n]/A] /(j[0,n]^3 BesselJ[1, j[0,n]]),,  {n,1,}]

(valid for 0 <= r < A).           (10.13)

Including the time-dependence for the product solutions solving the heat equation yields:    u[t,r,theta,m,n] = BesselJ[Abs[m], r j[m,n]/A] Exp[I m theta] Exp[-\lambda[m,n] t].


Let us work a specific initial-boundary-value problem.

Model Problem X.3. Given the Dirichlet boundary conditions as in Model Problem X.2 , suppose L=A=1 and that the initial temperature distribution is 100. What is the temperature at a later time?


The general solution is a superposition of the functions umn of (10.14),

   sum of c_mn u_mn(t,r,theta)

This is a double sum in n and m, analogous to the multiple Fourier series of Chapter IX. The coefficients are determined as in all the earlier cases by setting t=0 and arranging to satisfy the initial conditions:

   100 =  Sum[c[m,n]  BesselJ[Abs[m], r  j[m,n]/A]  Exp[I m theta].

Since the left side is independent of theta, only the terms with m=0 contribute, and we have a Fourier-Bessel series in one variable:

   100 =  Sum[c[0,n]  BesselJ[Abs[m], r  j[m,n]/A], {n,1,}].

The coefficients c0n could be calculated as integrals by using the orthogonality relationship for Bessel functions (10.11). This is, however, the hard way, since Fourier-Bessel series have a uniqueness theorem just as for Fourier series. Indeed, the powerful theorem of Chapter III applies in every respect to Fourier-Bessel series with only minor changes to accommodate the different orthogonality relationship. As was noted above, if the answer is unique, then it is right no matter how you find it.

In this case, the coefficients can be read off simply by comparison with the series given in (10.13):

    c[0,n] =  200 /(j[0,n] BesselJ[1, j[0,n]]).

The solution is thus:

   Sum[ (200 /(j[0,n] BesselJ[1, j[0,n]])) BesselJ[0, r j[0,n]/A] Exp[-l[0,n] t].


There is no dependence on theta, and the only value of m which contributes is m=0. These are accidental circumstances.

Model Problem X.4. As a second example, given the same BC as in X.3 and L=A=1, suppose that the initial temperature distribution is r sin(theta) + 1-r2. What is the temperature at a later time?

The general solution is still

   sum of c_mn u_mn(t,r,theta)

Set t=0 and check the initial conditions. The portion 1-r2 is orthogonal to all exp(i m theta) except for m=0, and by (10.13) the coefficients in its Fourier-Bessel expansion are

   c[0,n] = 8 /(j[0,n]^3 BesselJ[1, j[m,n]]).

There is no contribution from this term when m is not 0.

On the other hand, r sin(theta) contributes only terms with m=1, since
    sin(theta) = -(i/2)( exp(i theta) exp(-itheta) ),
and from (10.13) with p=1, we get:

   [1,n] = 2  /(j[1,n] BesselJ[1, j[1,n]]).

We get the full solution when we substitute this into (10.15). Only terms with m=0 and 1 are present in this solution.

Let us close the chapter with a three-dimensional problem, the heat equation on a cylinder. Setting k=1 for convenience, the heat equation takes on the form:

   u_t = partial ^2 u/partial r^2 + (1/r) partial  u/partial r + (1/r^2) partial ^2 u/partial theta^2 + partial ^2 u/partial z^2.           (10.16)

Precisely as before, we begin by searching for product solutions, separating vector variables one at a time. If

   u = T(t) U(r, theta,z),

then we find that

   -T'/T  = - Grad^2 U/U = lambda,          (10.17)

and the solutions for T are of the form

   T(t) = const. exp(-lambdat),

where lambda should be determined from the boundary conditions in space.

Model Problem X.5. Suppose that the ends of a cylinder of length L and radius A are insulated and the curved side is held at temperature 0. Set up the one-dimensional eigenvalue problems which make up the spatial part of the separated equation (10.17).


On the ends we have Neumann boundary conditions of the form

   partial u[t,r,theta,0]/partial z = partial u[t,r,theta,L]/partial z = 0,       (NBC)

while on the side we have Dirichlet boundary conditions of the form

   u(t,A, theta,z) = 0           (DBC)

The z part also separates off easily; if U = V(r, theta) Z(z), then we get:

    -lambda = (grad2 V)/V + Z''/Z,


    - Z''/Z = lambda + grad2 V)/V = constant =: mu.

The eigenvalue problem for Z gives us familiar eigenvalues and eigenfunctions:

    Z(z,k) = cos(k pi z/L),
mu(k) = (k pi/L)2.

The remaining part V is just like the V of Model Problem same set of Bessel functions:

   J_m(j_{mn} r/A).

The Z-dependence contributes to the eigenvalue, however, so that

Model Problem X.6 Refer to Prof. Herod's cheesecake recipe (with or without Maple). Eggs congeal at about 140 degrees. All the ingredients of the cheese cake before cooking are about 46 degrees. When the ingredients are mixed and placed in the hot oven, how long does it take to get the center of the cooking cake to 140 degrees?

Exercises X

X.1. Derive (10.1) by changing variables.

X.2. Solve the IBVP as in X.3, but with the initial conditions

   u(0,r, theta) = 2 - 2 r2 + r2 sin(2theta).

X.3. Solve the IBVP in Model Problem X.3, but with the initial conditions

   u(0,r, theta) = J0(r) theta2, - pi <= theta <= pi.

X.4. Solve Model Problem X.5, with A = L = 1 with the initial conditions

   u(0,r, theta) = (2 - 2 r2) z + r2 sin(2 theta) sin(z).

X.5. How a graduate student cooks a meal. Instead of cooking a steak, as the mathematician did in the previous chapter, an impoverished student with an untrained palate might instead cook a meal by putting a can of Spam directly in the oven. Assume that the thermal properties of Spam lead to a value of k=100 cm2/hr in the heat equation, and that a can is 4 cm high and 8 cm in diameter. It starts off uniformly at temperature 20 C and that the oven is held at temperature 200 C. There is good thermal contact on the entire surface of the cylinder of Spam. Calculate the temperature distribution of the Spam at all later times.

X.6. Suppose that Spam is an elastic material satisfying the wave equation with c = 1000 cm/sec, and that it is firmly encased in the can described above, which enforces zero Dirichlet BC on all faces, as in Problem IX.6.

a) Find the normal modes of vibration (in, say, the vertical direction).

b) Solve the IBVP with IC

   u(0,r, theta,z) = 1-r2
   ut(0,r, theta,z) = r sin(theta) z(4-z).

Link to