Linear Methods of Applied Mathematics

Evans M. Harrell II and James V. Herod*

version of 7 May 2000

In this chapter we solve some partial differential equations with spherical symmetry in three dimensions. Many of the calculations and supplementary material are available in the form of a Maple worksheet or Mathematica notebook.

(Some remarks for the instructor).

Conceptually, this chapter will be much like the previous one, but the geometry will now be three-dimensional. If there is a radial symmetry in a three-dimensional physical problem, the best coordinate system to use is probably the spherical system, with variables

r = distance from the origin

= colatitude

= longitude.

The colatitude is the angle measured down from the pole rather than up from the equator, which would be the usual geographic coordinate of latitude. Since we use radians for all angles, the relationship is

colatitude = /2 - latitude,

and the colatitude runs from 0 to . The longitude runs from 0 to 2 . You may wish to review the spherical coordinate system at this stage in your vector calculus textbook. (A summary is to be found in an appendix.)

Let us again begin with the potential equation of Laplace. A temperature could be prescribed on the outside of a ball as a function of the angular variables , producing a boundary-value problem of the form

with u(a, , ) = f( ,) (11.2)

To determine the temperature in the interior from measurements on the boundary, we must solve for u. If we try to find a product solution for (11.1) of the form

u(r, ) = R(r) Q( , ) ,

we readily come up with

This equation has some features of the multidimensional problem of Chapter IX and some of those of the problems on disks of Chapter X. To separate variables, we look at the last term. If we solve for it, we see that

equals a function of r alone; since it does not depend on r, it must be a constant. Thus we have a two-dimensional eigenvalue problem of the form

Equivalently, after multiplying through by - r^{2}Q,

(11.4)

Notice that the variable explicitly appears only in the term

If Q = ( ) (), we can further separate variables. Making this substitution and dividing by Q leads to

and all the dependence on is isolated in the term Since the rest of the equation does not involve , we conclude that this term is a constant, which leads us to our most familiar eigenvalue problem,

What about boundary conditions? The variable is periodic - if you increase your longitude by 2 radians, you will come completely around the earth to your starting point. This means that we have periodic boundary conditions as in Chapter X. This problem is so familiar that we can write down the solution from memory:

that is, m can be a positive or negative integer or zero. (Sines and cosines
would serve just as well, of course, but then m >= 0.) The eigenvalue
associated with
_{m} is m^{2}. Looking back at
(11.5) and
substituting
''/
= - m^{2}, we get an ordinary
differential equation for
, which a step or two of algebra puts into
the form

This is probably not striking the reader as the friendliest equation ever seen,
and even hours of substitution cannot reduce it to a familiar equation,
although some simplification is achieved by substituting z =
cos(),
P(z) =
(). The result is known as *
Legendre's
equation*,

and for general values of m and u its solutions are new transcendental
functions, just as the
Bessel functions were new functions defined as the
solutions of a new equation in the previous chapter. (A good, inexpensive
source of further information about special functions and equations like
Bessel's equation and
Legendre's
equation is the Dover book by N.N. Lebedev,
*Special Functions and their Applications*.) The variable z runs from -1
to +1, and Legendre's equation becomes singular at those points, which means
that some solutions are bounded at the end points while others diverge. In
effect we have a boundary condition at two end points, and the boundary
conditions determine the eigenvalues u. The simplest case is when m=0, and in
this case you may be astonished to know that the regular solutions of
Legendre's
equation are the Legendre polynomials P*l*, *l* = 0, 1,
...!
The eigenvalues are of the form u = *l *(*l *+ 1). (Exercise 1 provides some insight to the solution of Legendre's equation.) When m is not 0,
the eigenvalues are still of the form

u = *l *(*l *+ 1), with the restriction that *l* >= |m|,
and the bounded
solutions are the *associated
Legendre polynomials*,
defined for m>0 by

(11.7)

The usual Legendre polynomials result from setting m=0, and when m<0,

Fortunately, these functions are all known to Mathematica, where the notation is

While Mathematica can take care of doing calculations with special functions like the Legendre polynomials, you may find it helpful to read someLegendreP[l,x] LegendreP[l,m,x]

(11.8)

Some useful formulae are:

(11.10)

and

(11.11)

The one remaining variable to consider is r, and when Q is a spherical harmonic, (11.3) can be simplified to

This is again an
equidimensional equation,
so it has solutions of the form R =r^{a},
where a is a root of the polynomial

0 = a(a-1) + 2 a - *l* (*l* + 1),

that is a = *l* or -*l* - 1. The general solution for R is a
linear combination
of r^{l} and r^{-l} ^{-
1}. The latter solution is singular at r=0, so it can be excluded if the
differential equation is operative at the origin. (See
Exercise 4 for a
situation where this is not the case.) The general solution to
(11.1) can then
be written as

The coefficients c*l* m in the linear combination will be determined by
the
boundary conditions.

**Model Problem XI.1**. Solve
(11.1) on the interior of the ball of radius
1 centered at the origin, with the boundary condition

u(1,
,) =
(cos())^{2}

on its surface.

Solution.

The boundary data are azimuthally symmetric (independent of ), so the solution will also be independent of . Equivalently, all terms with m!=0 can be dropped from the general solution (11.12). This means that the spherical harmonics are just multiples of our old friends, the Legendre polynomials, with the variable z=cos(). The general solution with this simplification is

where the coefficients b_{l} differ from the coefficients c_{l 0} by a
normalization factor. That factor is unimportant, since at this stage we are
just as ignorant about the values of the c's as of the b's, but in case it is
reassuring,

The boundary corresponds to r=1, and in terms of the variable z, the function u
on the boundary is z^{2}. In order to determine the b's we need to
represent z^{2} as a
Legendre series, just as in
Chapter II. The
coefficients are obtained by projecting z^{2} onto the various Legendre
polynomials. Alternatively, because of the uniqueness theorem for Legendre
series, we can be confident that if we find the coefficients with a shortcut,
we have the correct answer. Recalling that P_{0}(z) = 1 and

we can easily determine that

z^{2} = (2/3) P_{2}(z) + (1/3) P_{0}(z).

In other words, b_{0} = 1/3 and b_{2} = 2/3, while all other coefficients are 0. The
solution to the boundary-value problem is:

u(1, ,) = 1 for 0 <= < , -1 for <= < 2 .

Solution.

The general solution is again as in (11.12), and we compare with the boundary conditions at r=1. We need to find the coefficients in the series

The function on the right side equals 1 or -1, depending on the value of . To determine the coefficients we multiply both sides of this expression by

and integrate by
sin()d d. Because of the
orthonormality of the spherical harmonics, the left side will reduce to the one
term, c*l'* m. Let us do the
integration first, since it is
exactly like the integrals we did in
Chapter IV
when finding the
Fourier series
for the square pulse. We find that

when m' is odd, and 0 when m' is even. For the integrations it is somewhat more efficient to use the variable z = cos(). we shall use Mathematica to calculate several of the coefficients:

And so on.In:= c[l_,m_] := (4/(I m)) Sqrt[(2 l + 1) (l-m)! /(4 Pi (l + m)!)] \ Integrate[LegendreP[l,m,z], {z,-1,1}] In:= c[1,1] 3 Pi Out= I Sqrt[----] 2 In:= c[3,1] 3 I 7 Pi Out= --- Sqrt[----] 16 3 In:= c[3,3] 5 I 7 Pi Out= --- Sqrt[----] 16 5 In:= c[5,1] 15 I 11 Pi Out= ---- Sqrt[-----] 64 30 In:= c[5,3] 35 I 11 Pi Out= ---- Sqrt[-----] 128 35 In:= c[5,5] 21 I 11 Pi Out= ---- Sqrt[-----] 128 7 In:= c[7,1] 175 I 15 Pi Out= ----- Sqrt[-----] 2048 14

Next, let us look again at a time-dependent problem, such as the heat equation,

As usual, we begin by looking for the
normal modes, product solutions with the
time-dependence separated: u(t,**x**) = T(t) U(**x**). Upon
substitution and dividing by u we get

so if

then we have T'(t) = - u k T(t). This is an elementary equation with solution

T(t) = T(0) exp(-u k t).

The values of u have to come from (11.13), which would be the Laplace's equation (11.1) if u=0. When u>0, we can still rely on our analysis of Laplace's equation up to a point, especially if the equation is defined on a domain with spherical symmetry.

The values of u will be determined as eigenvalues, once we have imposed boundary conditions.

**
Model Problem XI.3**. A mathematician has recently received tenure and
can now afford a
better cut of meat, the round steak.
Having too little spare
time to enroll in Cordon Bleu classes, he still cooks by the unsophisticated
method of pulling the round steak - an exact sphere 10 cm. in radius - directly out of
the freezer, at 0 C, and putting it into a preheated oven, at 200C.

We wish to find the temperature at the center of the steak after, say, an hour.
Since the steak is of higher quality than in
Chapter IX, let's assume that the
thermal constant k = 25 cm^{2}/hr.

Solution.

The boundary condition,

u(t,10, ,) = 200

is nonhomogeneous, so we modify the problem by redefining v = u - 200. The
initial condition for v will be v(0,**x**) = -200, and v will still solve
the wave equation. We shall solve for v, and then add 200 at the end. If we
seek product solutions v = T(t) U(**x**), equation
(11.13) for U will become
a well-posed
eigenvalue problem when we supplement it with the boundary
condition (now homogeneous):

U(10, ,) = 0 (11.14)

Instead of going through the entire separation procedure again, let's speed matters up by noting that the angular part of this equation is identical to that of (11.1). We can surmise from this that the normal modes will be of the form

Indeed, if we calculate the effect of the Laplace operator on this, we find

Because the spherical harmonics are the eigenfunctions of the angular part of the Laplace operator, they are still a common factor when the product solution is acted upon by the operator. The eigenvalue equation for u boils down to an ordinary differential equation for the function R:

(11.15)

Equation (11.15) should remind you of
Bessel's equation; the only difference is
the factor of 2 in the second term. As with
Bessel's equation, the term with u
prevents the equation from being equidimensional, and thus from being solved by
a simple r^{a}. It is in fact a form of
Bessel's equation, except
that when the extra factor of 2 is scaled out, the index gets messed up.
According to Mathematica:

In:= DSolve[r^2 R''[r] + 2 r R'[r] - (l (l+1) - mu r^2) R[r] == 0, R[r],r] Out= {{R[r] -> 1 2 > (BesselJ[Sqrt[- + l + l ], Sqrt[mu] r] C[1] + 4 1 2 > BesselY[Sqrt[- + l + l ], Sqrt[mu] r] C[2]) / Sqrt[r]}} 4

The subscript of the
Bessel
function simplifies to *l* + 1/2. As in the
previous chapter, we reject the Y Bessel functions because they are singular at
r=0. You may now be worrying that we are about to encounter still new
complications when investigating Bessel functions with fractional indices, but
take heart - when the index is half an odd integer, the Bessel functions reduce
to more familiar functions:

and, in general,

Here are some plots for J_{1/2}, J_{3/2}, and J_{5/2}.

Plot of J_{1/2}

Plot of J_{3/2}

Plot of J_{5/2}

As with the
Bessel functions of integer index, these functions oscillate and
damp down, and we can determine the eigenvalues, at least numerically. If
*l* =0, then it is clear that J_{1/2}(z) = 0 when z = n
, but otherwise
it is best to find them numerically. We use Mathematica and begin the
numerical calculations by estimating the positions of the
zeroes from the
graphs of the Bessel functions:

In order to satisfy the boundary condition at r=10, the eigenvalue u must be such thatIn:= FindRoot[{BesselJ[3/2,j[3/2,1]] == 0,BesselJ[3/2,j[3/2,2]] == 0, \ BesselJ[3/2,j[3/2,3]] == 0,BesselJ[3/2,j[3/2,4]] == 0}, \ {j[3/2,1],4}, {j[3/2,2],7}, {j[3/2,3],11},{j[3/2,4],15} ] Out= {j[3/2,1] -> 4.493409457909064175, j[3/2,2] -> 7.725251836937706892, j[3/2,3] -> 10.90412165942889983, j[3/2,4] -> 14.0661939127779873} FindRoot[{BesselJ[5/2,j[5/2,1]] == 0,BesselJ[5/2,j[5/2,2]] == 0, \ BesselJ[5/2,j[5/2,3]] == 0,BesselJ[5/2,j[5/2,4]] == 0}, \ {j[5/2,1],6}, {j[5/2,2],9}, {j[5/2,3],12},{j[5/2,4],16} ] {j[5/2,1] -> 5.76345919689110109, j[5/2,2] -> 9.095011330476355156, j[5/2,3] -> 12.32294097056643961, j[5/2,4] -> 15.51460301083520084}

where j_{l +1/2} is the n-th zero of a
Bessel function, which was
called j[l + 1/2,n] in the Mathematica code. Solving for the eigenvalues,

Although the eigenvalue could in principle depend on the index m as well, it does not happen to in this problem. The eigenfunctions of the Laplace operator are thus:

(11.17)

and the general solution to the heat equation in the ball is

(11.18)

The coefficients are determined by comparison with the initial condition, which is:

The standard way to proceed would be to use the orthogonality properties of the
spherical harmonics and the
Bessel functions. The spherical harmonics are
chosen to be orthonormal with respect to the angular measure
sin()
d d,
as we know. Since the eigenfunctions are orthogonal when
integrated over the volume of the ball, by
r^{2 }sin()
d d dr,
the orthogonality relationship of these
Bessel functions
will incorporate the extra r^{2}:

unless n=n' .

For this problem there is an easier way. If we note that there is no dependence on or , we can tell that the initial condition is orthogonal to all the spherical harmonics except for

since this spherical harmonics is a constant function, the best approximation
to the constant initial condition is proportional to it. The effect of this is
that the indices *l* and m are set to 0. The same argument does not
apply to the r variable, since none of the
Bessel functions is a constant. The
Fourier Bessel series we need is:

because of (11.16). If we multiply this by r, it becomes a Fourier sine series, of the form

We can find these coefficients by using the standard formula for Fourier coefficients,

(It is better to work directly with the coefficients bIn:= b[n_] = (2/10) Integrate[(-200 r) Sin[n Pi r/10], {r,0,10}] Out= n 4000 (-1) ---------- n Pi

so

To calculate the temperature the center of the steak, we encounter a 0/0 form, so it is useful to recall the limit sin(x)/x -> 1 as x -> 0. The numerical value from the first 7 terms in the series is thus:

In:= N[200 + Sum[(-1)^n 400 Exp[-n^2 Pi^2 /4], {n, 1, 7}]] Out= 166.099

**XI.1**. Noticing that
Legendre's equation is unchanged by changing z to -
z, we can conclude that a fundamental pair of solutions can be found, one of
which is even and the other odd. Suppose that we try to solve the equation
with a power series in z:

a) Substitute this into the equation and find the *recurrence relation*
expressing the coefficient an+2 in terms of an. There will be no mixing of
even and odd terms, because of the symmetry mentioned above.

At least one way to be sure that the power series converges to a function which
remains finite as

z -> +/-1 is to arrange for the series to terminate, that
is, for all the coefficients to become zero after a_{N} for some N. This will
make the solution a polynomial. It is not obvious that these constitute all
the solutions for this equation which remain bounded for -1<=z<=1, but
that is true.

b) For simplicity, let m=0, and find the conditions which make the series
terminate. Use this to derive the eigenvalue u = *l *(*l* + 1).
Verify that the first few
Legendre polynomials satisfy the Legendre equation
with these values of u.

**XI.2**. In
Model Problem XI.1, suppose you do not jump to the conclusion
that the terms with m different from 0 can be dropped. Use the orthogonality properties of
the functions exp(im)
to show more carefully that those coefficients must
be 0.

**XI.3**. Solve
Model Problem XI.2, but change the orientation of the
hemispheres so that the boundary conditions are

u(1, , ) = 1 for 0 <= < /2, -1 for /2 <= <= .

**XI.4**. Suppose that the ball in question is hollow, so that 1 < r
< 2, and that boundary conditions are posed on both the inner and outer
surfaces, i.e.,

u(1, ,) and u(2, ,) are both specified.

a) Find the general solution appropriate to this problem. Is there any reason
to reject the functions

r ^{-l - 1} when separating
variables?

b) For definiteness, solve Laplace's equation subject to

u(1, ,) = 0,

u(2,
,) = sin^{2}() when for 0 <=
<
/2,

u(2,
,
) = - sin^{2}()
when for
/2 <=
<=
.

**XI.5**. Find the general solution to the
heat equation with k= 3 on

a) the hemisphere, 0 <= r < 4, 0 <= <= /2, with mixed BC:

u(t,4, ,) = 1

u_{
} (t,r,
/2,) = 0

a) the hemisphere, 0 <= r < 4, 0 <= <= /2, with mixed BC:

u_{r}(t,4,
,) = 0

u(t,r, /2,) = 0

Hint for both parts: first take care of the nonhomogeneous BC by modifying the problem, and then find the eigenvalues and eigenfunctions for the Laplace operator on the hemisphere with the homogeneous BC.

**XI.5**. Find the solution to the initial-boundary-value problem

0 <= r <= 4, 0 <= <= /2

u(t,4, ,) = 1

u_{
}(t,r,
/2,) = 0

u(0,r, ,) = 1 + cos()

**XI.6**. Find the general solution to the heat equation with k= 3 on the
orange slice,

0 <= r <= 4, 0 <= <= /6, with Dirichlet boundary conditions:

u(t,4, ,) = 0

u(t,r, ,0) = u(r, , /6) = 0

**XI.7**. Find the general solution of the
wave equation in a sphere of
radius 10 cm, with characteristic speed c=1000 cm/sec.

a) Suppose that 0 Dirichlet boundary conditions are imposed on the sphere.

b) Suppose that Neumann boundary conditions are imposed on the sphere.

Onward to chapter XII

Back to chapter X

Return to Table of Contents

Return to Evans Harrell's home page