XI. Great Balls of PDEs

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 7 May 2000

XI. Great Balls of PDEs

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
    theta = colatitude
    phi = 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 theta , phi producing a boundary-value problem of the form

(11.1)

with u(a, theta , phi ) = f( theta , phi ) (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, theta ) = R(r) Q( theta , phi ) ,

we readily come up with

(11.3)

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²Q,

(11.4)

Notice that the variable phi explicitly appears only in the term

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

(11.5)

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

What about boundary conditions? The variable phi 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 Phi _m is m². Looking back at (11.5) and substituting Phi ''/ Phi = - m², we get an ordinary differential equation for theta , 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( theta ), P(z) = Theta ( theta ). The result is known as Legendre's equation,

(11.6)

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 Pl, 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


	LegendreP[l,x]
	LegendreP[l,m,x]

While Mathematica can take care of doing calculations with special functions like the Legendre polynomials, you may find it helpful to read some further discussion of the Legendre polynomials when doing the exercises. The term spherical harmonic refers to an eigenfunction of the Laplace operator, which is a normalized product of Legendre polynomials with exp(i m phi

(11.8)

Some useful formulae are:

(11.9)

(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

(11.12)

The coefficients cl 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, theta , phi ) = (cos( theta ))²

on its surface.

Solution.

The boundary data are azimuthally symmetric (independent of phi ), so the solution will also be independent of phi . 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( theta ). 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². In order to determine the b's we need to represent z² as a Legendre series, just as in Chapter II. The coefficients are obtained by projecting z² 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₀(z) = 1 and

we can easily determine that

z² = (2/3) P₂(z) + (1/3) P₀(z).

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

Model Problem XI.2. A classic problem is electrostatics is to find the potential inside a sphere when one hemisphere is kept at potential +1 and the other at 0. This can be solved in two different-looking ways, depending on the orientation of the hemispheres. Here I will solve the problem in one orientation, leaving the other to the exercises (yours is the easier way). The electrostatic potential solves Laplace's equation (11.1), and we may take as our boundary conditions:

u(1, theta , phi ) = 1 for 0 <= phi < , -1 for <= phi < 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 phi . To determine the coefficients we multiply both sides of this expression by

and integrate by sin( theta )d theta d phi . Because of the orthonormality of the spherical harmonics, the left side will reduce to the one term, cl' m. Let us do the phi 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 theta integrations it is somewhat more efficient to use the variable z = cos( theta ). we shall use Mathematica to calculate several of the coefficients:

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

And so on.

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

(HE)

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

(11.13)

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²/hr.

Solution.

The boundary condition,

u(t,10, theta , phi ) = 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, theta , phi ) = 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:

(11.16)

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:=
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}

In order to satisfy the boundary condition at r=10, the eigenvalue u must be such that

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( theta ) d theta d phi , as we know. Since the eigenfunctions are orthogonal when integrated over the volume of the ball, by r²sin( theta ) d theta d phi dr, the orthogonality relationship of these Bessel functions will incorporate the extra r²:

unless n=n' .

For this problem there is an easier way. If we note that there is no dependence on theta or phi , 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,


In:=
b[n_] = (2/10) Integrate[(-200 r) Sin[n Pi r/10], {r,0,10}]

Out=
              n
     4000 (-1)
     ----------
        n Pi

(It is better to work directly with the coefficients b_n rather than the c_00n along with its additional coefficients, which just complicate matters.) When we bring back the t dependence, we have

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

Exercises XI

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 phi ) 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, theta , phi ) = 1 for 0 <= theta < /2, -1 for /2 <= phi <= .

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, theta , phi ) and u(2, theta , phi ) 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, theta , phi ) = 0,

u(2, theta , phi ) = sin²( phi ) when for 0 <= theta < /2,

u(2, theta , phi ) = - sin²( phi ) when for /2 <= phi <= .

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

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

u(t,4, theta , phi ) = 1

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

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

u_r(t,4, theta , phi ) = 0

u(t,r, /2, phi ) = 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 <= theta <= /2