Unit 5: Boundary Value Problems for Partial Differential Equations

Chapter 29: Separation of Variables in Non-Cartesian Coordinates

Section 29.2: Laplace's equation in a cylinder

Copyright

Copyright * 2001 by Addison Wesley Longman, Inc.

All rights reserved. No part of this publication may be reproduced, stored in a retrieval system, or transmitted, in any form or by any means, electronic, mechanical, photocopying, recording, or otherwise, without the prior written permission of the publisher. Printed in the United States of America.

Initializations

> restart;

> with(linalg):
with(plots):
with(plottools):
with(student):
with(DEtools):

Warning, the protected names norm and trace have been redefined and unprotected

Warning, the name changecoords has been redefined

Warning, the names adjoint and translate have been redefined

> libname:="C:/Program Files/Maple 6.01/Fps",libname:

> with(FPS):

>

Problem Statement and Formulation

Find the steady-state, bounded temperatures inside a solid cylinder of radius and height , if the base is insulated, and the curved lateral surface and the top are maintained at temperatures of zero and , respectively.

Because none of the data depends on the angle , the temperature will be a function of the height and radius . Thus, the temperature function, which can be written as in cylindrical coordinates, satisfies Laplace's equation inside the cylinder. The insulation condition on the bottom of the cylinder, a homogeneous Neumann condition, is expressed by the vanishing of the temperature gradient, a partial derivative in the -direction, on that face. The complete boundary value problem determining therefore consists of the partial differential equation

that is,

> q := expand(laplacian(u(r,z),[r,theta,z], coords=cylindrical)) = 0;

>

and the boundary conditions

bounded

>

Statement and Analysis of the Solution

The solution is given by

=

where is the Bessel function of order zero, and is the k th zero of .

The functions satisfy Bessel's equation

that is,

> Bessel_eqn := x^2*diff(y(x),x,x) + x*diff(y(x),x) + (x^2-nu^2)*y(x) = 0;

>

and are called Bessel functions of order . As seen in Section 16.2, two linearly independent solutions of Bessel's equation are and , as we see from the Maple solution

> dsolve(Bessel_eqn,y(x), output=basis);

>

Only one of these solutions, , is bounded at the origin. The Bessel function of the second kind, , has a logarithmic singularity at .

Series expansions for and , the Bessel functions of order zero and one, respectively, can be found with the FormalPowerSeries command from the FPS package. The series expansion of is

> FormalPowerSeries(BesselJ(0,x),x);

>

while the series expansion of is

> FormalPowerSeries(BesselJ(1,x),x);

>

In Figure 29.4, is the black curve, and is the red curve.

> plot([BesselJ(0,x),BesselJ(1,x)],x=0..10, color=[black,red], scaling=constrained, xtickmarks=5, ytickmarks=2, labels=[x,``], labelfont=[TIMES,ITALIC,12]);

>

It should be recalled from Section 16.2 that if , are the roots of , then the functions are, on the interval , orthogonal with respect to the weight function . Thus, for , we have

>

Example 29.2

For the cylinder described above, let , , and . Then, the first ten eigenvalues are the first ten zeros of , listed as in Table 16.1 in Section 16.2.

For the convenience of the reader executing these calculations in Maple, we reconstruct these first ten eigenvalues as follows.

As seen in Section 16.2, is asymptotic to , so the eigenvalues , the first ten roots of , are computed by

> for k from 1 to 10 do
m[k] := fsolve(BesselJ(0,2*x)=0, x, (4*k-1)*Pi/8-.1 .. (4*k-1)*Pi/8+.1);
od;

>

The coefficients , given by

can be evaluated with numeric integration, yielding, for example, as

> evalf(Int(r*BesselJ(0,m[1]*r),r=0..2))/2/ cosh(3*m[1])/BesselJ(1,2*m[1])^2;

>

Alternatively, because , the integral in the numerator can be evaluated exactly via the formula

if the change of variables is made in the integral in the numerator of , as summarized by

In fact, the integral

> unassign('k'):
X := Int(r*BesselJ(0,lambda[k]*r),r=0..2);

>

under that change of variables, becomes

> Y := changevar(x=lambda[k]*r,X,x);

>

Maple knows the value of this integral, and gives

> value(Y);

>

In fact, Maple can also evaluate the original integral, giving

> value(X);

>

Consequently, the expression for

> Ak := Int(r*BesselJ(0,lambda[k]*r),r=0..2)/2/ cosh(3*lambda[k])/BesselJ(1,2*lambda[k])^2;

>

simplifies to

> AAk := value(Ak);

>

that is, to

Hence, the first ten coefficients are

> for k from 1 to 10 do
A[k] := evalf(subs(lambda[k]=m[k],AAk));
od;

>

The coefficients rapidly decrease in magnitude, so the convergence is rapid, and the partial sum

that is,

> U10 := add(A[k]*BesselJ(0,m[k]*r)*cosh(m[k]*z),k=1..10);

>

is adequate, except on the boundary where a much slower convergence is evident, as seen in Figure 29.5.

> p1 := plot3d(U10,r=0..2,z=0..3, color=red):
display(p1, scaling=constrained, labels=[r,z,u], labelfont=[TIMES,ITALIC,12], axes=frame, tickmarks=[2,3,2], orientation=[-40,60]);

>

The surface represents the temperatures on the rectangle R = { , }, half a plane section through the axis of the cylinder. It is impossible to draw a graph of the temperatures inside the cylinder. Every representation of these temperatures will require some visual compromise.

In Figure 29.6, the temperatures are graphed on the rectangle R , shown supporting the temperature surface inside the cylinder.

> p2 := cylinderplot(2,theta=Pi/2..2*Pi,z=0..3):
p3 := display(rotate(p1,-Pi/2,0,-Pi/2)):
display([p2,p3], orientation=[30,70], axes=frame, tickmarks=[5,5,4], labels=[` x`,` y`,`z `], labelfont=[TIMES,ITALIC,12]);

>

The rectangle rotates about the axis of the cylinder, supporting the same temperatures for all values of the angle .

Figure 29.7 gives another view of the temperatures inside the cylinder, the isothermal surface on which .

> PP := implicitplot3d(U10=.1,r=0..2,theta=0..2*Pi,z=0..3, coords=cylindrical, scaling=constrained, orientation=[35,60], color=yellow, style=patchcontour, axes=frame, tickmarks=[5,5,4], grid=[20,20,20]):
PP1 := spacecurve({[[0,0,0],[2,0,0]],[[0,0,0],[0,-2,0]],[[0,0,0], [0,0,1]]}, color=black, linestyle=2):
PP2 := textplot3d([-.1,-.1,.1,`O`]):
display([PP,PP1,PP2],labels=[` x`,`y `,`z `], labelfont=[TIMES,ITALIC,12], orientation=[-45,65]);

>

In fact, the level surfaces , are each bowl-shaped, with the same circle bounding the open top. This corresponds to the discontinuity in temperatures at the junction of the top of the cylinder with the curved lateral wall of the cylinder. At this juncture, all level surfaces intersect.

Executing the following Maple commands will generate an animation of this sequence of level surfaces through the interior of the cylinder.

> for k from 1 to 9 do
P||k := implicitplot3d(U10 = .1*k, r=0..2, theta=0..2*Pi, z=0..3, coords=cylindrical, color=COLOR(HUE,k/10)):
od:
display([P||(1..9)], insequence=true, orientation=[25,80], style=patchcontour, scaling=constrained, lightmodel=light3);

>

Solution by Separation of Variables

Assuming a separated solution of the form

> U := R(r)*Z(z);

>

and substituting it into the boundary data and leads to

=>

and

Z'(0) = 0 => Z'(0) = 0

Application of the separation assumption to Laplace's equation

> q;

>

results in

> q1 := eval(subs(u(r,z)=U,q));

>

As usual, division by

> q2 := expand(q1/U);

>

leads to a separation of variables if the term containing is brought to the righthand side via

> q3 := q2 - (op(3,lhs(q2))=op(3,lhs(q2)));

>

Introducing the Bernoulli separation constant leads to

=

Overlooking the abuse of notation whereby prines denote differentiation with respect to both and , we have the two ODEs

and

which follow from

> q4 := lhs(q3) = mu;
q5 := rhs(q3) = mu;

>

and the rearrangements

> q6 := numer(normal(lhs(q4) - rhs(q4))) = 0;
q7 := -numer(normal(lhs(q5) - rhs(q5))) = 0;

>

The -equation is a modified form of Bessel's equation. In fact, the solution is

that is,

> dsolve(q6,R(r));

>

The appearance of suggests renaming the separation constant to so that the two ordinary differential equations become

and

that is,

> q8 := subs(mu=-lambda^2,q6);
q9 := subs(mu=-lambda^2,q7);

>

Then, the solution of the -equation is

that is

> dsolve(q8,R(r));

>

Since we know , the zero-order Bessel function of the second kind, is not bounded at the origin, its coefficient is taken as zero, and the solution for is a multiple of .

To transform the -equation to Bessel's equation

that is,

> Bessel_eqn;

>

make the change of variables so that becomes = . Done in Maple, we obtain

> q10 := Dchangevar(r=x/lambda,R(r)=y(x),q8,r,x);

>

Factoring and "cancelling" via

> expand(q10/lambda);

>

gives precisely Bessel's equation of order zero.

Incidentally, the change of variables was accomplished by a careful use of the chain rule starting with the definitions and

=

Then,

and

=

so that

R'' + R' + R = 0

becomes

y'' + y' + y = 0

and hence

y'' + y' + y = 0

when is replaced by .

Since , the boundary condition requires the solution of for the eigenvalues .

The -equation

that is,

> q9;

>

now inherits the known eigenvalues, and becomes

that is,

> unassign('k');
q11 := subs(lambda=lambda[k],q9);

>

which has exponential solutions. The solution satisfying the Neumann boundary condition Z'(0) = 0 is

> q12 := dsolve({q11,D(Z)(0)=0},Z(z));

>

which is a multiple of a hyperbolic cosine, as seen by

> simplify(convert(q12,trig));

>

Consequently, a single eigensolution is

that is,

> Uk := A[k]*BesselJ(0,lambda[k]*r)*cosh(lambda[k]*z);

>

and the general solution, a sum over all such eigensolutions, is

that is,

> U := Sum(Uk,k=1..infinity):
u(r,z) = U;

>

Application of , the final nonhomogeneous Dirichlet boundary condition at the top of the cylinder, leads to the equation

that is,

> f(r) = subs(z=h,U);

>

from which the coefficients are determined via multiplication by and integrating with respect to over the interval . The orthogonality of the eigenfunctions with respect to the weighting function causes all terms but one on the right side to vanish. The one surviving term is the one for which . Hence, what results is the equation

> q13 := Int(r*f(r)*J[0](lambda[j]*r),r=0..sigma) = A[j]*cosh(lambda[j]*h)*Int(r*J[0](lambda[j]*r)^2,r=0..sigma);

>

from which is determined to be

> isolate(q13,A[j]);

>

In Section 16.2, the equivalence of the integral in the denominator with ( ) was explored in Maple..

>