Partial Differential Equations PowerTool

by Dr. Jim Herod

Section 7.3: Heat Equation on a Disk

     We consider the diffusion equation on a disk. The physical situation is that of either a long circular beam, or a disk with insulated top and bottom, but with a fixed temperature on the sides, and an initial temperature throughout the disk. This can be written as

          1/k  =  u,

with u(t, a,  ) = f(  ),   0 < t,     <  <  

and    u(0, r,  ) = g( r,  ),     < a,      <  <  .

     To do this problem, we first solve the steady state problem

          0 =  v,  v(a,  ) = f(  ),    <  <  .

Then solve the transient problem

          1/k  =  u,

with u(t, a,  ) = 0,  0 < t,    <  <  

and    u(0, r,  ) = g( r,  ) - v(r, )  ,  < a,    <  <  .

     Finally, the solution for the original problem is the sum of these two.

     We have already studied how to solve the first problem in Sections 6.4 and 6.5. Here is how to think of solving the second problem. We write out the partial differential equation in unabbreviated form to facilitate separation of variables.

                        = .

(We have dropped the constant of diffusion for purposes of illustration.) Separation of variables leads to these three differential equations with boundary conditions:

(1)    '' +  = 0,    = ,    '(  ) =  '( ),

(2)   T ' = -  T,  and

(3)   r ( r R ' ) ' -  R = -    R,  R(a) = 0,  R is bounded on [0, a).

We know how to solve (1) and (2). Equation 3 leads to a study of Bessel's functions.

     To begin this study, we start with an example that is independent of . Here is the problem we solve

                 =  u,

with   u(t, a,  ) = 0,  0 < t,    <  <  

and    u(0, r,  ) = g( r),     < a,    <  <  .

Is it clear that the solution of this equation will be independent of  ? Thus, the partial differential equation reduces to

                   = .

Now we have the two differential equations

(1) T ' +  T = 0, and

(2) r ( r R ' ) ' +    R = 0, with R(a) = 0, and R bounded on [0, a).

     We examine the second differential equation. Expanded, it is

            R '' + r R ' +    R = 0.

We ask Maple to solve this equation. Commonly, the solutions are introduced in elementary differential equations in a discussion of series. We are, as you know, assuming that  > 0. Maple would be assuming only that it is complex, unless we tell it to do otherwise.

>	assume(lambda>0);

>	dsolve(r^2*diff(R(r),r,r)+ r*diff(R(r),r)+lambda^2*r^2*R(r)=0,R(r),series);

Of course, Maple knows the solutions in the full details, without truncated series. We ask for solutions, again, and do not ask that the solutions be as series.

>	dsolve(r^2*diff(R(r),r,r)+ r*diff(R(r),r)+lambda^2*r^2*R(r)=0,R(r));

It seems we have two solutions. To understand the Bessel functions, we draw their graphs.

>	plot([BesselY(0,r),BesselJ(0,r)],r=0..10,color=[BLACK,RED]);

There are several observations to make:

(1) BesselJ(0, r) is one at zero and BesselY(0, r) is unbounded at zero.
(2) The two solutions oscillate, are not periodic, and have intertwining zeros.
(3) These form eigenfunctions corresponding to the eigenvalues -  for the differential operator

          L(R) = 1/r  ( r R ' ) '.

Here's a check of this last.

>	R:=r->BesselY(0,lambda*r); 1/r*(diff(r*diff(R(r),r),r))+lambda^2*R(r); simplify(%);

>	R:=r->BesselJ(0,lambda*r); 1/r*(diff(r*diff(R(r),r),r))+lambda^2*R(r); simplify(%);

     It is of value to compare these two functions with the cosine and sine functions, which define the eigenfunctions for the problem  '' +    = 0. The comparison helps to crystallize the structure in your memory.  Since we are interested in only bounded solutions on the disk, we set aside BesselY.  

     Suppose a > 0. Take J(r) to be BesselJ

Comparison 1:  There is an infinite increasing sequence  ,    ,    ,  ...  such that J(  a) = 0. For the cosine function, these numbers were  .

Comparison 2:  The numbers  -  are eigenvalues and R(r) = J(  r) is an eigenfunction for the operator L(R)  defined above with R(a) = 0. The function cos(  x ) was an eigenfunction for  '' with  (a) = 0.

Comparison 3:  Orthogonality persists in the sense that

            = 0 if    and     =  J '  .

Comparison 4:  If f(r) is sectionally smooth on (0, a) then

          f(r) =  where  = <f(r), J(  r) >/

In this case, the dotproduct is defined by     =  with the associated norm.

It seems appropriate to do an example. We take the boundary conditions to be zero so that we will not have to do the steady state problem first. We take the initial condition to be  so that the initial condition is independent of . With a = 1, we will have a continuous function on the unit disk.

      Here is a graph of the initial condition.     

>	f:=r->1-r^2;

>	cylinderplot([r,theta,f(r)],theta=-Pi..Pi,r=0..1,axes=NORMAL, shading=zhue);

>

     Your expectation is that the initial distribution will decay to the steady state u(r,  ) = 0. Here are the details.

The first thing to do is to get the zeros of R(r) = BesselJ(r). We can create a program to compute these. Go back to the graph and see how far the zeros are apart. Then, solve for a zero in an interval.

>	zero[0]:=0; for n from 1 to 4 do      zero[n]:=fsolve(BesselJ(0,x)=0,x,zero[n-1]+1..zero[n-1]+4); end do; n:='n';

Pretty good. No surprise. These zeros are so important that Maple already knows them.

>	for n from 1 to 10 do     zero[n]:=evalf(BesselJZeros(0,n)); end do; n:='n';

>

To make further comparisons, we plot the first six of the eigenfuntions cos(  x) and the first six of the eigenfunctions BesselJ(0,  x).

>	plot([seq(cos(n*Pi*x),n=0..5)],x=0..1,color=[black,red,green,blue,yellow, brown]);

>	plot([seq(BesselJ(0,zero[n]*x),n=0..5)],x=0..1,color=[black,red,green,blue,yellow, brown]);

Let's check the orthogonality. In order to meet the boundary conditions that R is bounded and R(a) is zero, we need the positive zeros of the BesselJ(0, x) function. What you should expect is that we don't get exactly zero, for each zero is only an approximation.

>	for n from 1 to 3 do   for m from 1 to 3 do      print(int(BesselJ(0,zero[n]*x)*BesselJ(0,zero[m]*x)*x,x=0..1));   end do; end do; n:='n':m:='m':

Now is time to compute the coefficients.

>	for n from 1 to 10 do    A[n]:=int(f(x)*BesselJ(0,zero[n]*x)*x,x=0..1)/                int(BesselJ(0,zero[n]*x)*BesselJ(0,zero[n]*x)*x,x=0..1); end do; n:='n':

The next two graphs are so close, I offset the first one by 0.007.

>	plot([f(x)+0.007,sum(A[n]*BesselJ(0,zero[n]*x),n=1..10)],x=0..1,        color=[black,red]);

We now create the solution.

>	u:=(t,r)->sum(A[n]*BesselJ(0,zero[n]*r)*exp(-zero[n]^2*t),n=1..10);

>	cylinderplot([r,theta,u(1,r)],r=0..1,theta=-Pi..Pi,axes=NORMAL, shading=zhue);

>	animate3d([r,theta,u(t,r)],r=0..1,theta=-Pi..Pi,t=0..1,             coords=cylindrical, frames=30, shading=zhue);

We remark in passing that the following is an alternate way to do a cylinder plot that does not require a call of the package "plots".

>	plot3d([r,theta,u(0,r)],r=0..1,theta=-Pi..Pi,coords=cylindrical, shading=zhue);

>

EMAIL: herod@math.gatech.edu   or   jherod@tds.net

URL: http://www.math.gatech.edu/~herod

Copyright ©  2003  by James V. Herod

All rights reserved