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

Mathematica.

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 producing a boundary-value problem of the form

             (10.1)
with
    u(a, ) = f()          (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 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, ) = R(r) Q(), we readily find

   

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 and + 2 , we must have:

   Q( + 2 ) = Q() (PBC)

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

Q(,m) = a(m) cos(m ) + b(m) sin(m ),                          (10.3)

   R'' + (1/r) R' - (m²/r²) R = 0 (10.5)

Equidimensional equations can be solved by making the guess R = r^a 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:

                           (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:

              (10.7)

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

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

Solution.

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:

a₀ + Sum( a_m cos(m ) + b_m sin(m ) = (cos())²

    cos(2 ) = 2 cos²() - 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    

or in other words that a₀ = 1/2, a₂ = 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:

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

              (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, ),

then we find that

    T'/T = (Grad² 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 - . The solutions for T are of the form

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

where 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).

Solution.

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

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

the eigenvalue need not be zero. The constant 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(),

leads to the same eigenvalue equation for Q() 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

with the knowledge that Q() = sin(2 m ) or cos(2 m ), 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' + (- m²/r²) 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 r^m or r^-m, but the latter are unphysical because the temperature can't diverge at r=0. The regular solutions, behaving roughly like r^m near r=0, are called regular Bessel functions and, with a standard choice of the overall constant, denoted J_m(^1/2 r). Mathematica uses the notation BesselJ. They are normalized so that J_mx) ~ x^m/(2 m!) as x -> 0:

In:=
  Plot[BesselJ[0,x],{x,0,25}]

Out=

In=
  Plot[BesselJ[1,x],{x,0,25}]
Out=

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 so 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

   j_m,n = the n-th positive value of x for which J_m[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} ]

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

In=
  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} ]

Out=
  {j[1, 1] -> 3.83171, j[1, 2] -> 10.1735, j[1, 3] -> 10.1735, 
   
    j[1, 4] -> 19.6159}

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

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 x/L or 2n x/L which appear in Fourier series, there is an expression involving j_m,n.

Orthogonality.

              (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 d.

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, ) found above can be written either with sines and cosines, or with exponential functions, and I shall make the latter choice:

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

Here, m can be positive, negative, or 0. With the normalization, the orthogonality relationship is:    

          (10.12)

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:

   

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

   

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

          (10.14)

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?

Solution.

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

   

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:

   

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

   

The coefficients c_0n 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):

The solution is thus:

   

          (10.15)

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

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

   

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

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

   

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:

              (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, ,z),

then we find that

             (10.17)

and the solutions for T are of the form

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

where 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).

Solution.

On the ends we have Neumann boundary conditions of the form

          (NBC)

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

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

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

    - = (² V)/V + Z''/Z,

or,

    - Z''/Z = + ² V)/V = constant =: .

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

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

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

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

_mnk =(j_mn/A)² + (k/L)².

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, ) = 2 - 2 r² + r² sin(2).

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

   u(0,r, ) = J₀(r) ², - <= <= .

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

   u(0,r, ) = (2 - 2 r²) z + r² sin(2 ) 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 cm²/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, ,z) = 1-r²
   u_t(0,r, ,z) = r sin() z(4-z).

• chapter XI
• chapter IX
• Table of Contents
• Evans Harrell's home page
• Jim Herod's home page