College of Engineering and Computer Science

Mechanical Engineering Department

Notes on Engineering Analysis
Larry Caretto January 21, 2009

Solution of the Diffusion Equation

Introduction and problem definition

The diffusion equation describes the diffusion of species or energy starting at an initial
time, with an initial spatial distribution and progressing over time. The simplest example
has one space dimension in addition to time. In this example, time, t, and distance, x,
are the independent variables. The solution is obtained for t 0 and 0 x xmax. The
dependent variable, species or temperature, is identified as u(x,t) in the equations
below. The basic diffusion equation is written as follows.

u 2u
2 [1]
t x

Here, , is a species or thermal diffusion coefficient with dimensions of length squared

over time. The initial condition specifies the value of u at all values of x at t = 0. This
initial condition is usually written as follows:

u(x,0) = u0(x) [2]

The boundary conditions at x = 0 and x = xmax, can be functions of time in general.

These are written as follows (with subscripts L and R denoting left and right).

u(0,t) = uL(t) u(xmax,t) = uR(t) [3]

There is actually a discontinuity in the initial conditions and the boundary conditions.
The initial conditions are assumed to hold for all times less than or equal to zero. The
boundary conditions hold for t > 0. This discontinuity only occurs at the boundaries. For
example, the initial condition may be that u = 0 for all x and the boundary conditions
may be that u = 100 at both boundaries for t > 0. This creates a discontinuity at t = 0 for
x = 0 and x = xmax. Fortunately, this discontinuity does not affect our ability to solve the
problem.

For the simplest cases we will consider both the initial and boundary conditions to be
constant. However, in the most general case the boundary conditions may be given by
any functional forms for u0(x). uL(t), and uR(t)

Basic solution process

The solution of the diffusion proceeds by a method known as the separation of

variables. In this method we postulate a solution that is the product of two functions,

Jacaranda (Engineering) 3333 Mail Code Phone: 818.677.6448

E-mail: lcaretto@csun.edu 8348 Fax: 818.677.7062
T(t) a function of time only and X(x) a function of the distance x only. With this
assumption, our solution becomes.

u(x,t) = X(x)T(t) [4]

We do not know, in advance, if this solution will work. However, we assume that it will
and we substitute it for u in equation [1]. Since X(x) is a function of x only and T(t) is a
function of t only, we obtain the following result when we substitute equation [4] into
equation [1].

u X ( x)T (t ) T (t )
X ( x)
t t t
[5]
2u 2 X ( x)T (t ) 2 X ( x)
2 T (t )
x x 2 x 2

If we divide the final equation through by the product X(x)T(t), we obtain the following
result.

1 1 T (t ) 1 2 X ( x)
[6]
T (t ) t X ( x ) x 2

The left hand side of equation [6] is a function of t only; the right hand side is a function
of x only. The only way that this can be correct is if both sides equal a constant. This
also shows that the separation of variables solution works. In order to simplify the
solution, we choose the constant1 to be equal to2. This gives us two ordinary
differential equations to solve.

1 1 dT (t ) 1 d 2 X ( x)
2 2 [7]
T (t ) dt X ( x) dx 2

dT (t )
The first equation becomes 2T (t ) . The general solution to this equation is
dt
2
known to be2 T (t ) Ae t . The second differential equation in [7] can be written as
d 2 X ( x)
2
2 X ( x) 0 . The general solution to this equation is X(x) = Bsin(x) +
dx
Ccos(x). Thus, our general solution for u(x,t) = X(x)T(t) becomes

u ( x, t ) T (t ) X ( x) Ae t B sin(x) C cos(x) e t C1 sin(x) C 2 cos(x)

2 2
[8]

1
The choice of 2 for the constant as opposed to just plain comes from experience. Choosing
the constant to have this form now gives a more convenient result later. If we chose the constant
to be simply , we would obtain the same result, but the expression of the constant would be
awkward.
2
As usual, you can confirm that this solution satisfies the differential equation by substituting the
solution into the differential equation.

Jacaranda (Engineering) 3333 Mail Code Phone: 818.677.6448

E-mail: lcaretto@csun.edu 8348 Fax: 818.677.7062
We first consider a simple set of boundary conditions where u(0,t) = 0 and u(x max,t) = 0.
In order to satisfy these boundary conditions on u for all times, we have to have
boundary conditions on X(x) that X(0) = X(xmax) = 0. With these homogenous boundary
conditions on X(x) and the form of the original differential equation for X(x), the solution
for X(x) is the solution of a Strum-Liouville problem. Accordingly, expect that we will be
able to obtain an infinite set of solutions for X(x) that will give a complete set of
orthogonal eigenfunctions.

If we substitute the boundary condition at x = 0 into equation [8], get the following result
because sin(0) = 0 and cos(0) = 1.

u (0, t ) 0 e t C1 sin(0) C 2 cos(0) C 2 e t

2 2
[9]

Equation [9] will be satisfied for all t if C2 = 0. With C2 = 0, we can apply the solution in
equation [8] to the boundary condition at x = x max.
2
u ( x max , t ) 0 C1e t sin(x max ) [10]

Equation [10] can only be satisfied if the sine term is zero. 3 This will be true only if xmax
is an integral times . If n denotes an integer, we must have

n
xmax n or [11]
xmax

We have an infinite number of solutions for X(x), as the integer n goes from 1 to infinity.
We can write an individual solution as Xn(x) = sin(nx/xmax). Because the differential
equation and boundary conditions for X(x), form a Sturm-Liouville problem, we know
that the solutions to this problem, Xn(x), are an infinite set of orthogonal eigenfunctions.
Since any integral value of n gives a solution to the original differential equations, with
the boundary conditions that u = 0 at both boundaries, the most general solution is one
that is a sum of all possible solutions, each multiplied by a different constant, C n. This
general solution is written as follows

n
u ( x, t ) Cn e nt sin(n x )
2
with n [12]
n 1 xmax

We still have to satisfy the initial condition that u(x,0) = u 0(x). Substituting this condition
for t = 0 into equation [12] gives

nx nx
u 0 ( x) u ( x,0) C n e 0 sin
xmax


Cn sin xmax


[13]
n 1 n 1

3
Of course, the boundary condition is also satisfied if C1 = 0, but this gives u(x,t) = 0, not a very
interesting solution.

Jacaranda (Engineering) 3333 Mail Code Phone: 818.677.6448

E-mail: lcaretto@csun.edu 8348 Fax: 818.677.7062
Since the eigenfunctions sin(nx/xmax) are a complete set of orthogonal eigenfunctions
in the interval x x xmax, we can expand any initial condition function, u 0(x) in terms of
these eigenfunctions. We can obtain the usual equation for eigenfunction expansion
coefficients of Sturm-Liouville solutions by using the orthogonality relationships for
integrals of the sine. If we multiply both sides of equation [13] by sin(mx/xmax), where
m is another integer, and integrate both sides of the resulting equation from a lower limit
of zero to an upper limit of xmax the upper and lower limits of the Sturm-Liouville
problem for X(x) we get the following result.
x max x max
mx mx nx

0
u0 ( x) sin
xmax
dx C
0 n 1
n sin
xmax
sin
xmax
dx
[14]
x max x max
mx nx mx
Cn sin sin dx Cm sin
2
dx
n 1 0 xmax xmax 0 xmax

In the second row of equation [14] we reverse the order of summation and integration
because these operations commute. We then recognize that the integrals in the
summation all vanish except for the one where n = m, because the eigenfunctions,
sin(nx/xmax), are orthogonal. This leaves only one integral to evaluate. Solving for C m
and evaluating the last integral4 in equation [14] gives the following result.
x max
mx
u0 ( x) sin
x
max
dx
2
x max
mx
Cm 0
x max
mx

xmax u 0 ( x) sin
x
max
dx [15]
sin dx
2 0

0
x
max

For any initial condition, then, we can perform the integral on the right hand side of
equation [15] to compute the values of Cm and substitute the result into equation [12] to
compute u(x,t). We recognize that equation [15] is the usual equation for eigenfunction
expansions in Sturm-Liouville problems.

As an example, consider the simplest case where u0(x) = U0, a constant. In this case
we find Cm from the usual equation for the integral of sin(ax).
x max x max x max
2 mx 2U 0 mx 2U 0 xmax mx
Cm
xmax u 0 ( x) sin
xmax
dx
xmax sin
xmax
dx
xmax

m
cos
xmax



0 0 0

4
Using a standard integral table, and the fact that the sine of zero and the sine of an integral
multiple of are zero, we find the following result:
x max x max
2 mx x x 2mx xmax xmax 2mxmax xmax
sin
x
max
dx max sin
2 4m

xmax

2

4m
sin
xmax
0
2
0 0

Jacaranda (Engineering) 3333 Mail Code Phone: 818.677.6448

E-mail: lcaretto@csun.edu 8348 Fax: 818.677.7062
 2U 0 xmax mxmax xmax m0 2U 0
Cm cos cos 1 cos m
xmax m x max m x
max m
[16]

Since cos(m) is 1 when m is even and 1 when m is odd, the factor of 1 cos(m) is
zero when m is even and 2 when m is odd. This gives the final result for C m shown
below.

4U 0
m m odd

Cm
0 m even

[17]

Substituting this result into equation [12] gives the following solution to the diffusion
equation when u0(x) = U0, a constant, and the boundary temperatures are zero.
2
4U e nt sin(n x) n
u ( x, t ) 0 with n [18]
n 1,3,5 n xmax

If we substitute the equation for n into the summation terms we get the following result.

n 2 2 t
2 nx
e x max
sin
4U 0
xmax [19]
u ( x, t )


n 1, 3, 5 n

Equation [19] shows that the ratio u(x,t)/U 0 depends on the dimensionless parameters
t/(xmax)2 and x/xmax. Thus we can compute a universal plot of the solution as shown in
Figure 1. Because the solution is symmetric about the midpoint x/L = 0.5, only the left
half is plotted here.

Although we can obtain reasonable solutions using equation [19] for large times, we
will not get a correct solution, within the limits of practical computers, if we try to
evaluate equation [19] for times near zero. Figure 2 shows the results of using equation
[19] to compute u(x,0)/U0 using various numbers of terms for t = 0. In this case we are
trying to use the sine series to evaluate X(x) = 1. We see that taking more terms makes
the solution more accurate away from x = 0, but no matter how many terms we take,
there are still strong oscillations near x = 0. This result is well known and is called the
Gibbs phenomena.

Nonzero boundary conditions

The use of zero boundary conditions to give us a Sturm-Liouville problem appears to be

a major limitation on this approach. We can solve a broader class of problems, still
using the same approach, by appropriate definitions of variables.

Jacaranda (Engineering) 3333 Mail Code Phone: 818.677.6448

E-mail: lcaretto@csun.edu 8348 Fax: 818.677.7062
First we consider the case were the boundary potentials, u(0,t) and u(x max,t) are both set
to a fixed value uB. In this case we can make the following simple variable
transformation from u(x,t) to v(x,t).

Jacaranda (Engineering) 3333 Mail Code Phone: 818.677.6448

E-mail: lcaretto@csun.edu 8348 Fax: 818.677.7062
 v(x,t) = u(x,t) uB [20]

This transformation gives v(0) = u(0,t) uB = uB uB = 0 and v(xmax) = u(xmax,t) uB = uB

uB = 0. Thus, v has homogenous boundary conditions.

Since uB is a constant, when we substitute equation [20] into equation [1], we get the
same equation in terms of v instead of u. Since the boundary conditions on v are that v
= 0 at x = 0 and x = xmax, the differential equation for u is that of equation [1], and
boundary conditions for v are the same as those used to derive the solution for u in
equation [12]. Accordingly, equation [12] is the solution for v(x,t). When we apply the
initial condition that u(x,0) = u0(x) to equation [12] for this revised problem, we get the
following equation in place of equation [13].

nx
u0 ( x) u ( x,0) v( x,0) u B Cn sin [21]
n 1 xmax

Starting with this equation, we can repeat the derivation of equation [15] for C m, which
gives the following result for this case where both x = 0 and x = x max boundaries have
the same common boundary value, uB.
x max x max

u 0 ( x) uB sin mx dx 2
2 mx cos m 1
Cm
xmax
0 xmax xmax
0
u 0 ( x) sin
xmax
dx uB
m
[22]

Thus, the solution for u(x,t) for the general initial condition, u(x,0) = u 0(x) and the
boundary conditions that u(0,t) = u(xmax,t) = uB is given by the following equation.

2 x max

u ( x) u sin x dx
2
u ( x, t ) u B 0 B n e nt sin(n x ) [23]
xmax
n 1 0

We can extend the procedure used above to consider the case where the two
boundaries are set at two different constant temperatures. That problem is specified by
the following equations, where uL and uR are constants.

u 2u
2 with u ( x,0) u0 ( x), u (0, t ) uL , u ( xmax , t ) u R [24]
t x

The solution to this problem proceeds by defining u(x,t) as the sum of two functions
v(x,t) and w(x), both of which satisfy the differential equation in [24].

u(x,t) = v(x,t) + w(x) [25]

Again, the function v(x,t) is defined to have a zero value at x = 0 and x = x max, to give us
the Sturm-Liouville problem, which provides us with a complete set of eigenfunctions.
This zero boundary condition on v(x,t) tells us what the boundary conditions on w(x)
must be set to satisfy the boundary conditions in equation [24].

Jacaranda (Engineering) 3333 Mail Code Phone: 818.677.6448

E-mail: lcaretto@csun.edu 8348 Fax: 818.677.7062
 w(0) = u(0,t) v(0,t) = uL 0 = uL w(xmax) = u(xmax,t) v(xmax,t) = uR 0 = uR [25a]

Setting w(0) = uL and w(xmax) = uR means that the solution for u(x,t) proposed in equation
[25] satisfies the boundary conditions of equation [24]. Of course the proposed solution
must also satisfy the differential equation. Since both v(x,t) and w(x) satisfy the
differential equation in [24], their sum, u(x,t) = v(x,t) + w(x), which is the solution that we
seek, will also satisfy the differential equation. We then examine the solution to the
differential equations for both v(x,t) and w(x).

Since w(x) does not depend on t, we have the following result when we substitute w(x)
into the differential equation.

w 2w d 2w
0 2 0 [26]
t x dx 2

That is, we have a simple, second-order, ordinary differential equation to solve for w. 5
Integrating this equation twice gives the following general solution for w(x).

w(x) = C1x + C2 [27]

Since this equation has to satisfy the boundary conditions that w(0) = u L and w(xmax) =
uR, we have the following solution for w(x) that satisfies the differential equation and
boundary conditions.

uR uL
w( x ) u L x [28]
xmax

Since v is defined to be zero at both boundaries, the solution for v(x,t) is the same as
the solution originally found for u(x,t) in equation [12]. Thus the solution for u(x,t) =
v(x,t) + w(x) is found by adding equation [12] for v(x,t) and equation [28] for w(x) to give.

nx uR uL
Cne t sin
2
u ( x, t ) v ( x, t ) w( x ) n u L x [29a]
n 1
xmax xmax

Setting t = 0 in this equation gives the initial condition for u(x,t).

nx uR uL
u ( x,0) v ( x,0) w(0) Cn sin xmax
u L
xmax
x [29b]
n 1

Applying the general initial condition that u(x,0) = u 0(x) gives the following result.

5
The w function does not depend on time. If the solution for v(x,t) is such that this function
approaches zero for long times, then u(x,t) = v(x,t) + w(x) becomes equal to w(x). Thus, w(x)
represents the steady state solution and we can regard the definition of u(x,t) = v(x,t) + x(x,t) as
dividing the solution into a transient part that decays to zero and a steady-state part.

Jacaranda (Engineering) 3333 Mail Code Phone: 818.677.6448

E-mail: lcaretto@csun.edu 8348 Fax: 818.677.7062

nx uR uL
u0 ( x) Cn sin xmax
u L
xmax
x [29c]
n 1

Repeating the derivation of equation [15] for C m, in the same manner used to find Cm for
the boundary condition of two equal, non-zero boundary values in equation [22], gives
the following result in this case.
x max
2 uR uL mx
Cm
xmax u 0 ( x) u L
xmax
x sin
xmax
dx
0
x max [30]
2 mx 2u 1 cos m 2 u R u L cos m
xmax u 0 ( x ) sin
x
max
dx L
m

m
0

Thus, the solution to the differential equation and boundary conditions in equation [24] is
given by the following equation, where equation [30] is used to find the values of C n.

uR uL
n
x Cn e nt sin(n x)
2
u ( x, t ) u L with n [31]
xmax n 1 xmax

2
We see that the e nt terms in the summation approach zero for large times leaving the
solution as uL + (uR uL)x/xmax which is the steady state solution, defined in equation [25]
as w(x).

As in the case of zero boundaries, the easiest problem to examine is the one where
u0(x) is a constant, U0. In this case equation [30] becomes
x max
2 mx 2u 1 cos m 2 u R u L cos m
Cm
xmax U 0 ( x) sin
xmax
dx L
m

m [32]
0
2U 0 u L 1 cos m 2 u R u L cos m

m m

Since cos m = 1 for even m and -1 for odd m, we have the following result for C m.

2 u R u L
m even
m
Cm [33]
4U 0 u L 2 u R u L
m odd
m m

We can substitute this result into equation [31] and divide by U 0 uL to get the following
dimensionless result.

Jacaranda (Engineering) 3333 Mail Code Phone: 818.677.6448

E-mail: lcaretto@csun.edu 8348 Fax: 818.677.7062

2
n
u ( x, t ) u L u R u L x 2 1 t nx

U 0 u: L xmax n 2, 4, 6, n
e xmax
sin
U 0 u: L max
x
[34]
2
n
2 u uL
1 t nx


2 R
U 0 u: L

n 1, 3, 5, n
e xmax
sin
xmax

In this equation the dimensionless ratio, [u(x,t) u L] / (U0 uL) depends on the
dimensionless distance x/xmax, dimensionless time, = t/(xmax)2, and the ratio of initial
and boundary conditions, (uR uL) / (U0 uL). A plot of this result for (uR uL) / (U0 uL)
= 0.5 is shown in Figure 3. This plot shows that the steady-state solution is a linear
profile with the expected boundary values for [u(x,t) u L] / (U0 uL); these are zero at x
= 0 and (uR uL) / (U0 uL) = 0.5 at x = xmax. If we changed the value of the parameter
(uR uL) / (U0 uL) from 0.5 to some other value, the plot would be similar, but the
location of the dimensionless parameter at x = x max would shift.

Figure 3 Solution of diffusion equation for nonzero boundaries

Cylindrical geometry

In a cylindrical coordinate system, 0 r R, the diffusion equation has the following

form.

u 1 u
r [35]
t r r r

Jacaranda (Engineering) 3333 Mail Code Phone: 818.677.6448

E-mail: lcaretto@csun.edu 8348 Fax: 818.677.7062
The most general initial and boundary conditions for the radial problem are

u(r,0) = u0(r) u/r|r=0,t = 0 u(R,t) = uR(t) [36]

We will first consider the conditions where u R is a constant. We will later see that we will
have to use an expansion of Sturm-Liouville eigenfunctions to match the initial
condition. Anticipating this, we will obtain a problem with zero boundary conditions by
splitting u(r,t) into two parts, as we did in equation [20].

u(r,t) = v(r,t) + uR [37]

The new function, v(r,t) satisfies the differential equation in [35] and has v/r = 0 at r =
0 and v = 0 at r = R. Since uR is a constant, the proposed solution in [37] will satisfy the
original differential equation and the boundary conditions from equations [35] and [36].
As before, we seek a solution for v(r,t) using separation of variables. We postulate a
solution that is the product of two functions, T(t) a function of time only and P(r) a
function of the radial coordinate, r, only. With this assumption, our solution becomes.

v(r,t) = P(r)T(t) [38]

We substitute equation [38] for u in equation [35]. Since P(r) is a function of r only and
T(t) is a function of t only, we obtain the following result.

u P(r )T (t ) T (t )
P (r )
t t t
1 u 1 P(r )T (t ) 1 P(r ) [39]
r r T (t ) r
r r r r r r r r r

If we divide the final equation through by the product P(r)T(t), we obtain the following
result.

1 1 T (t ) 1 1 P ( r )
r 2 [40]
T (t ) t P ( r ) r r r

The left hand side of equation [40] is a function of t only; the right hand side is a function
of r only. The only way that this can be correct is if both sides equal a constant. As
before, we choose the constant to be equal to2. This gives us two ordinary
differential equations to solve.

dT (t )
2T (t ) has the general solution T (t ) Ae t .
2
The first equation becomes
dt
P( r )
The second differential equation in [40] can be written as r 2 rP (r ) 0 . The
r r
general solution to this equation is P(r) = BJ0(r) + CY0(r), where J0 and Y0, are the

Jacaranda (Engineering) 3333 Mail Code Phone: 818.677.6448

E-mail: lcaretto@csun.edu 8348 Fax: 818.677.7062
Bessel functions of the first and second kind with zero order.6 Thus, our general
solution for u(x,t) = X(x)T(t) becomes

v(r , t ) T (t ) P (r ) Ae t BJ 0 (r ) CY0 (r ) e t C1 J 0 (r ) C2Y0 (r )

2 2
[41]

As r 0, Y0(r) -; to keep the solution finite, we require that C2 = 0. (We can also
show that the first derivative of J0 is zero at r = 0 so that we satisfy the zero gradient
condition in equation [36]. To satisfy the condition that v(R,t) = 0, we must satisfy the
following equation.
2
v( R, t ) 0 C1e t J 0 (R) [42]

Equation [42] defines a transcendental equation similar to that in equation [10] which
required sin(xmax) = 0. That equation has a simple solution since we know that the
sin(n) = 0 when n is an integer. We do not have such a simple result for the equation
that J0(R) = 0. However, the points at which J0 = 0 can be determined and are
available in tables. We use the symbol mn to indicate the mth point where Jn is zero.
That is we define mn by the following equation.

J n ( mn ) 0 for m 1,, [43]

For example, the first five points at which J0 is zero are for 10 = 2.40483, 20 = 5.52008,
30 = 8.65373, 40 = 11.79153, and 50 = 14.93092. There are an infinite number of such
values such that J0(mr), with mR = m0 provides a complete set of orthogonal
eigenfunctions that can be used to represent any other function over 0 r R. These
mn values are called the zeros of Jn.

We can now write the solution to our radial diffusion, from equation [38], as the sum of
v(r,t) and uR.

Cm e
2
u (r , t ) m t J 0 ( m r ) u R m R m0 [44]
m 1

6
Bessel functions, like sines and cosines, are solutions to differential equations. Although Bessel
functions occur less frequently than sines and cosines, they are functions that can be found in
tables or in computer function calculations. The general form of Bessels equation used for
d 2 y ( x) 1 dy ( x) x 2 2
obtaining series solutions is y ( x) 0 . This equation can be
dx 2 x dx x2
d dy 2 2
transformed into
dz dz z k z y 0 whose solution is AJ(kz) + BY(kz). We see that
z

P( r )
this second equation has the same form as r 2 rP (r ) 0 , provided that we set
r r
= 0. This gives the result above that the solutions are J 0 and Y0. The factor of r in the 2rP(r) term
is a weighting factor that must be included in the definition for the inner product of solutions to the
radial portion of the diffusion equation.

Jacaranda (Engineering) 3333 Mail Code Phone: 818.677.6448

E-mail: lcaretto@csun.edu 8348 Fax: 818.677.7062
We still have to satisfy the initial condition that u(r,0) = u 0(r). We can do this by using an
eigenfunction expansion.

u ( r ,0) u R u 0 ( r ) u R C m J 0 ( m r ) [45]
m 1

The values of Cm are found from the general equation for orthogonal eigenfunction
expansions which includes a weighting function. Here the weighting function p(r) is
equal to r. For any initial condition, u0(r), we find the values of Cm from the following
equation. (In the second equation below, the integral in the denominator has been
evaluated using integral tables and the fact that J 0(mR) = 0.)
R R

rJ 0 (m r ) u0 (r ) u R dr rJ 0 (m r ) u0 (r ) u R dr
Cm 0
R
0
[46]
R2
r J 0 (m r ) dr
2
J1 (m R) 2
0
2

If we consider the simple case where u0(r) is a constant equal to U0, we have the
following result for Cm.
R rR
rJ1 ( m r )
rJ 0 ( m r )U 0 u R dr
0 2U 0 u R m r 0
Cm
R 2
J1 ( m R ) 2 R 2
J1 ( m R ) 2

RJ1 ( m R)
2U 0 u R
m 2U 0 u R 2U 0 u R [47]
Cm
R J1 ( m R )
2 2 m RJ1 ( m R) m0 J1 ( m0 )

With this result for Cm our solution for u(r,t) with the constant initial condition is written as
shown below after some algebra.

2U 0 u R m2 0 R 2
t
r
u (r , t ) e J 0 m0 uR [48]
m 1 m 0 J1 ( m 0 ) R

This equation shows the important variables are r/R and t/R 2. We can rearrange this
equation to get a dimensionless form, plotted in Figure 4.
t
u (r , t ) u R
2 m2 0 2
r
e R
J 0 m0 [49]
U 0 uR m 1 m 0 J 1 ( m 0 ) R

Jacaranda (Engineering) 3333 Mail Code Phone: 818.677.6448

E-mail: lcaretto@csun.edu 8348 Fax: 818.677.7062
 Figure 4 Solutions to radial diffusion with constant initial condition

In Figure 4, the notation tau denotes the dimensionless time t/R 2. The results in
Figure 4 are similar to those in Figure 1 for the rectangular geometry. Here the plot
from 0 to R covers half the cylinder and is equivalent to the plot in Figure 1 that covered
only half of the one-dimensional region.

Mixed boundary condition

Boundary conditions on differential equations can be expressed in terms of the value of

u at the boundary, the gradient of u at the boundary, or a combination of the two. The
boundary condition which is an equation relating the value of u and its gradient at the
boundary is called a mixed boundary condition.

In this example we consider a problem of one dimensional heat conduction. Here we

will explicitly use the temperature, T, as the dependent variable (instead of u) and note
that in this case is called the thermal diffusivity. We consider the case where a slab,
initially at a constant temperature, T0, is placed in a medium at a different temperature,
T. Convective heat transfer occurs on each side such that the results in the slab are
symmetric about the center. In such a case, we can solve for only half the geometry. If
we assume that the slab has a thickness of 2L, where L x L, we solve for a region
between 0 and L using a symmetry boundary condition that the gradient is zero at the
center line. The convection boundary condition at x = L simply states that the
conduction heat flux in the solid at the boundary is the same as the convection heat flux
leaving the boundary. This is written as follows.

Jacaranda (Engineering) 3333 Mail Code Phone: 818.677.6448

E-mail: lcaretto@csun.edu 8348 Fax: 818.677.7062
 dT
qx L k h(TL T ) [50]
dx xL

For this problem, we will define the following dimensionless variables that will give us a
Sturm-Liouville problem for the x solution in the separation of variables result.

T T x t
[51]
T0 T L L2

With these definitions, the diffusion equation and the initial and boundary conditions
may be written in the following dimensionless form.

2 hL
( ,0) 1 0 1 0 [52]
2 0 1 k

We can perform our usual separation of variables solution to obtain the following
general solution.

( , ) e C1 sin( ) C2 cos( )
2
[53]

The derivative of this solution is shown below.

e C1 cos( ) C2 sin( )
2

[54]

The gradient / will be zero at = 0 only if C1 is zero. With C1 = 0, the boundary

condition at = 1 becomes

hL hL 2
1 e C2 sin( )
2
e C2 cos( ) 0 [55]
1 k k

We can rearrange this equation to obtain a transcendental equation for .

k cos( )
cot( ) [56]
hL sin( )

The roots of this equation can be visualized as the points were the plot of a straight line
through the origin with a slope of k/hL intersects with the plot of the cotangent. Such a
plot is shown in Figure 5 for hL/k = 1.

Jacaranda (Engineering) 3333 Mail Code Phone: 818.677.6448

E-mail: lcaretto@csun.edu 8348 Fax: 818.677.7062
 Figure 5 Eigenvalues for convection problem

In this problem, the eigenfunctions are cos(mx), but the values of m do not have a
simple equation such as they did in previous problems with simpler boundary
conditions. Instead the values of m for this problem, similar to those for the
eigenvalues for J0, are not evenly spaced and have to be found by numerical solutions.
The first few values of m for hL/k = 1 are7 1 = 0.860333589, 2 = 3.425618459, 3 =
6.437298179, 4 = 9.529334405, 5 = 12.64528722, 6 = 15.77128487, and 7 =
18.90240996. Changing the value of hL/k would change the slope of the line in Figure 5
and consequently change the intersection points. As m increases, the values of m
increase and the intersections come closer to an integer value of where the cotangent
goes to infinity. For example, 7 = 18.90240996 is close to 6 = 18.85.

Even though the eigenvalues, m, are not evenly spaced, the part of the separation-of-
variables problem that was solved for the dependence is a Sturm-Liouville problem, so
we know that we can expand any function, representing the initial condition, in terms of
these eigenfunctions.

As usual, the most general solution for is a sum of all the eigenfunctions that satisfy
the differential equations and the boundary conditions.

hL
( , ) Cme m cos(m )
2
m cot(m ) [57]
m 1 k

7
You can observe that these values of l are the intersections indicated by the red dots in Figure 5.

Jacaranda (Engineering) 3333 Mail Code Phone: 818.677.6448

E-mail: lcaretto@csun.edu 8348 Fax: 818.677.7062
If, for the moment, we regard T0 as some arbitrary measure of a variable initial
temperature profile, we can use an eigenfunction expansion to fit any initial condition,
(,0) = 0() = (T0(L) T)/(T0 - T).

hL
( ,0) 0 ( ) Cm cos(m ) m cot(m ) [58]
m 1 k

We have the usual equation for the eigenfunction expansion, where we can find the
integral in the denominator from integral tables.
1 1

0 ( ) cos(m )d 2m 0 ( ) cos(m )d
Cm 0
0
[59]
1
cos(m ) sin(m ) m
cos ( )d
2
m
0

In particular for a constant initial temperature, T0, 0 = 1, and Cm is found as follows.

1
1
sin(m )
2m cos(m )d 2m
m 2 sin(m ) [60]
Cm 0
0

cos(m ) sin(m ) m cos(m ) sin(m ) m cos(m ) sin(m ) m

Substituting this expression for Cm into equation [57] gives our solution for the
dimensionless temperature, , assuming a constant initial temperature.
2

2 sin(m )e m cos(m ) hL
( , ) m cot(m ) [61]
m 1 cos(m ) sin(m ) m k

The solution is plotted for hL/k = 1 in Figure 6 and for hL/k = 10 in Figure 7.
Both solutions are different from the earlier ones in which the boundary temperature
was set to a new value at zero time. Here, the boundary temperature decreases over
time because of convection heat transfer with the environment. With hL/k = 10, the
convection is ten times as strong as it is with hL/k = 1. This leads to a more rapid
change in temperature and a greater temperature gradient (i.e., greater heat transfer) at
the x = L boundary.

Gradient boundary condition

Consider the diffusion equation with zero gradient boundary conditions shown below.

u 2u u u
2 0 xL 0 u x,0 f ( x) [62]
t x x x 0, t x x L,t

Jacaranda (Engineering) 3333 Mail Code Phone: 818.677.6448

E-mail: lcaretto@csun.edu 8348 Fax: 818.677.7062
The same process that was used to solve equation [1] can be applied to equation [62] to
obtain the following solution.

Jacaranda (Engineering) 3333 Mail Code Phone: 818.677.6448

E-mail: lcaretto@csun.edu 8348 Fax: 818.677.7062
 2
n
t
nx
u ( x, t ) C 0 C n e L
cos [63]
n 1 L

The constants in equation [63] are given by the following equation.

L L
1 2 mx
L 0
C0 f ( x)dx For m 1 : Cm f ( x) cos dx [64]
L0 L

As t , the exponential terms will vanish and the only term that will be left in the
L
1
solution is the C0 term. This gives the result that u x, t C 0 f ( x)dx .
L
0
Physically, this result tells us what happens when there is no flux into the region.
(Recall that the flux of physical quantities described by the diffusion equation is
proportional to the gradient; thus a zero gradient means no flux.) With no flux added to
the region, the final potential at all points in the region is simply the average of the initial
condition.

Nonzero gradients on two sides of a one-dimensional geometry

This section should be omitted at first reading. Go directly to the heading Mixed
boundary condition in a cylinder on page 24. If we want to consider two nonzero
gradients we have a problem that does not have a steady-state solution unless both
gradients are the same. Even if both gradients are the same the usual notion of a
solution in terms of two variables, u(x,t) = v(x,t) + w(x) will not have a unique solution.
To see this, consider the following problem.

u 2u u u
2 0 xL g0 gL u x,0 f ( x ) [65]
t x x x 0, t x x L ,t

A solution of the form u(x,t) = v(x,t) + w(x) where v(x,t) satisfies the diffusion equation
with zero gradient boundary conditions and w(x) satisfies the equation d 2w/dx2 = 0 with
the boundary conditions that dw/dx = g0 at x = 0 and dw/dx = gL at x = L will satisfy the
differential equation. However, the equation for w cannot satisfy its boundary
conditions. To see this, we integrate the equation d 2w/dx2 = 0 two times to obtain the
general solution w = A + Bx.

Taking the first derivative of this general solution gives du/dx = B. This shows that du/dx
is a constant which must be the same at both boundaries so that B = g 0 = gL. However,
after we find the value for B, we have no other boundary condition for the constant A.
Thus, any solution u = A + g0x = A + gLx will satisfy the differential equation and the
boundary conditions only when g0 = gL.

This mathematical result parallels physical reasoning. In the diffusion equation

derivation the gradient of u represents a flux. If we have a one dimensional region in

Jacaranda (Engineering) 3333 Mail Code Phone: 818.677.6448

E-mail: lcaretto@csun.edu 8348 Fax: 818.677.7062
which we have a different flux on two sides, there cannot be a steady solution. How can
we find a solution to the problem posed in equation [65]?

Since this problem does not admit a steady state solution we can augment our usual
division of u(x,t) into a function v(x,t) that satisfies the diffusion equation with zero
boundary conditions and a function, w(x), that depends only on the space coordinate.
We have to include another function, (t) that provides the long term solution. We write
our proposed solution as follows.

u(x,t)=v(x,t)+ w(x) + (t) [66]

Here v satisfies the diffusion equation with zero gradient boundary conditions, and the
dw
function w(x) has the following gradient boundary conditions: dx g 0 and
x 0

dw
gL . We can show that the solution in equation [66] satisfies the original
dx xL

gradient boundary conditions.

u v (t ) w
0 0 g0 g0 [67]
x x 0,t x x 0 ,t x x 0,t x x 0 ,t

u v (t ) w
0 0 gL gL [68]
x x L ,t x x L ,t x x L ,t x x L ,t

The initial condition for v can de derived from the initial condition for u:

v ( x,0) u ( x,0) (0) w( x ) f ( x ) (0) w( x ) [69

If we substitute equation [66] into the diffusion equation and note that w(x) is a function
of x only and (t) is a function of time only, we obtain the following result.

u 2 u v w 2 v w v 2v 2w
2 0 [70]
t x t x 2 t t x 2 x 2

Since v satisfies the diffusion equation, the v terms in the last expression cancel leaving
the following relationship between and w.

(t ) 2 w( x )
[71]
t x 2

We see that equation [71] has a function of t on one side and a function of x on the
other side. Just as in the solution of an original partial differential equation by
separation of variables, the only way that a function of time can equal a function of x is if
both functions equal a constant. Thus, we can rewrite this equation in the same way as
we usually apply separation of variables.

Jacaranda (Engineering) 3333 Mail Code Phone: 818.677.6448

E-mail: lcaretto@csun.edu 8348 Fax: 818.677.7062
 (t ) 2 w( x )
A [72]
t x 2

This gives two ordinary differential equations whose solutions are known.

d
A At B
dt
[73]
d 2w A Ax 2
w Cx D
dx 2 2

We can apply the boundary conditions on w to evaluate A and C.

dw Ax
g0 C C g0
dx x 0 x 0
[74]
dw Ax g C g g0
gL C A L L
dx x L x L L L

With the constants just found we can write the solutions for (t) and w(x) as follows.

g L g0 gL g0 2
(t ) At B tB w( x) x g0 x D [75]
L 2L

The solution for v(x,t) is the solution to the diffusion equation with zero gradient
boundary conditions. This solution is an infinite series in the cosine of nx/L, which was
given in equation [63].
2 2
n n
t
nx t
nx
v ( x, t ) C n e L
cos C0 C n e
L
cos [76]
n 0 L n 1 L

The solution for u(x,t) = v(x,t) + w(x) + (t) is then found by combining equations [73]
and [76].
2
n
t
nx
u ( x, t ) v( x, t ) w( x) (t ) C0 Cn e L
cos
n 0 L [77]
A 2
x Cx D At B
2

Before substituting the results of equation [74] for A and C, we will consider the initial
condition on v: v(x,0) = f(x) (0) w(x) = f(x) B (A/)(x2/2) Cx D. Satisfying this
initial condition requires

Ax 2 nx
v ( x ,0 ) f ( x ) B Cx D C0 Cn cos [78]
2 n 1 L

Jacaranda (Engineering) 3333 Mail Code Phone: 818.677.6448

E-mail: lcaretto@csun.edu 8348 Fax: 818.677.7062
We can determine the coefficients Cn by an eigenfunction expansion. We use a
separate integration for the n = 0 eigenfunction, cos(0) = 1, because it has a different
normalization constant.
L L L
Ax 2 nx
f ( x) B
2
Cx D (1) dx C 0 (1)dx


n 1
C n cos
L
dx
0 0 0
L L L [79]
A x x 3 2 L nx
f ( x )dx Bx
2 3
C
2
Dx C0 L Cn
n sin L

C0 L
0 0 n 1 0

We thus find the following result for C0

L
1 A L2 L
C0
L f ( x ) dx B
2 3
C D
2
[80]
0

To get the remaining values of Cn, we multiply both sides of the initial condition equation
by cos(mx/L) and integrate from 0 to L. Applying the orthogonality relationship for the
cosine eigenfunctions gives the following result for m > 0.
L
Ax 2 mx
f ( x ) B
2
Cx D cos
L
dx
0
L L [81]
nx mx mx C m L
C n cos
L
cos
L
dx C m cos 2
L

2
0 n 0 0

The integral for the terms B and D vanishes as shown below.

L L
mx mx
B D cos L
dx B D sin


L 0
0 [82]
0

Combining equations [81] and [82] gives the following expression for C m when m 1.
L
2 A 2 mx
Cm
L f ( x)
2
x Cx cos
L
dx m 1 [83]
0

If we substitute equations [80] and [83] into equation [77] we obtain the following
equation for u(x,t).
L
1 A L2 L
u ( x, t )
L f ( x )dx B
2 3
C D
2
0
2
n
2
L
A 2 mx t
nx [84]

L
f ( x)
2
x Cx cos
L
dx e

L
cos
L

n 0 0
A 2
x Cx D At B
2

Jacaranda (Engineering) 3333 Mail Code Phone: 818.677.6448

E-mail: lcaretto@csun.edu 8348 Fax: 818.677.7062
The B and D terms in this equation cancel so these terms do not enter the solution. 8
The following steps can be used to modify equation [84]: (1) use the equations in [74] to
1L
replace A and C; (2) use the definition of the average initial condition, U 0 f ( x)dx ,
L0
to replace the integral in equation [84]; and (3) use the symbol, C m, in place of the full
expression used in [84]. With these steps the following result obtains after some
algebra.
2
n
g L g0 g0 t
nx
u ( x, t ) U 0 L Cn e L
cos
6 2 n1 L [85]
g L g0 2 g g0
x g0 x L t
2L L

If we consider a problem with such a large value of t that the terms in the summation
are negligible, we can integrate equation [85] to find the (spatial) average value of u(x,t)
at any time. In this integration there is a large amount of cancellation leading to the
result shown below.

1L g g0
u (t )
L0
u ( x, t ) U 0 L
L
t [86]

In this case of large , the average value of u is increasing or decreasing with time
depending on the sign of gL g0.

For large times, the exponential term in the summation of equation [85] will vanish;
removing this summation from equation [85] gives the large time limit for the solution
shown in equation [87].

g L g0 g0 g L g0 2 g g0
Lim u ( x, t ) U 0 L x g0 x L t [87]
t 6 2 2L L

8
In many approaches to this problem the functions and w are defined so that (0) = w(0) = 0.
This drops the constants B and D early in the derivation. However, there is no physical reason for
setting these conditions; they are only a convenience to rewrite the derivation in a simpler manner
once we know the solution.

Jacaranda (Engineering) 3333 Mail Code Phone: 818.677.6448

E-mail: lcaretto@csun.edu 8348 Fax: 818.677.7062
 u g g0
At this point the change in u with time is L . This derivative is not a function
t L
of x or t.9 At any value of x then, the change in u(x,t) at any x location, u(x,t 2) u(x,t1),
will be the same for any given time difference, t 2 t1.

g L g0
Lim u ( x, t2 ) u ( x, t1 ) (t2 t1 ) [88]
t1 , t 2 L

At large times then, the spatial profile of u(x,t) will not change with x. At any value of x
the change in u from one time point to another will be the same. This change in u will
be an increase if gL > g0; it will be a decrease if gL < g0, and there will be no change in u
if gL = g0. In the case where g0 = gL, equation [87] gives

L
Lim u ( x, t ) U 0 g 0 x g L g0 [89]
t 2

In this case where the flux is the same constant at both boundaries, the long-time
potential is the average of the initial condition plus a linear gradient whose slope is the
constant flux term, g0 = gL.

Mixed boundary condition in a cylinder

The problem of diffusion in a cylindrical coordinate system, 0 r R, for a fixed

boundary condition at the outer radius was treated above, starting with equations [35]
and [36]. If there is a mixed boundary condition at the outer radius of the cylinder, the
initial and boundary conditions for this problem become.

u(r,0) = u0(r) u/r|r=0,t = 0 ku/r|r=R,t = h(u u) [90]

In order to solve this problem, we will first define a new variable, v(r,t) as we did in
equation [37], to give a zero boundary condition at r = R.

v(r,t) = u(r,t) u [91]

With this definition, the boundary condition at r = R becomes

kv/r|r=R,t = hv(R,t) or (k/h)v/r|r=R,t + v(R,t) = 0 [92]

The solution of the differential equation [35] proceeds in the same manner as done
previously, leading to the following equation, similar to equation [42].
9
A comparison of equations [86] and [87] shows that two related quantities have the same value:
u (t ) g g0 u ( x, t )
L Lim
t L t t
The derivative on the left is the average (over x) at any time; the derivative on the right is the
value at any x for large times. The two equations are consistent. At large times, the derivative at
any x has a constant value. Thus the average (over x) will have the same value. Equation [86]
tells us that the average (over x) will have this same value at any (large) time.

Jacaranda (Engineering) 3333 Mail Code Phone: 818.677.6448

E-mail: lcaretto@csun.edu 8348 Fax: 818.677.7062
 2
v( r , t ) Ce t
J 0 ( r ) [93]

We apply this solution to the boundary condition at r = R, using the result that dJ 0(x)/dx
= J1(x) to give the following result.

k v k 2 2
v ( R, t ) 0 Ce t J1 (R ) Ce t J 0 (R ) [94]
h r rR h

Dividing by the constant, C, and the exponential term and defining = R gives the
following results.

k hR
J1 () J 0 () 0 J 0 () J1 () 0 [95]
hR k

The roots of this equation give the eigenvalues = /R. These eigenvalues will depend
on the value of hR/k. The first ten eigenvalues for various values of hR/k are shown in
the table below.

Table of Eigenvalues, n = nR, for Various Values of hR/k

Index hR/k = hR/k = hR/k = hR/k = hR/k = hR/k = hR/k =
,n 0.001 0.01 0.1 1 10 100 1000
1 0.04471577 0.141245 0.441682 1.25578 2.17950 2.38090 2.40242
2 3.83197 3.83431 3.85771 4.07948 5.03321 5.46521 5.51456
3 7.01573 7.01701 7.02983 7.15580 7.95688 8.56783 8.64508
4 10.1736 10.1745 10.1833 10.2710 10.9363 11.6747 11.7797
5 13.3238 13.3244 13.3312 13.3984 13.9580 14.7834 14.9160
6 16.4707 16.4712 16.4767 16.5312 17.0099 17.8931 18.0530
7 19.6159 19.6164 19.6210 19.6667 20.0829 21.0036 21.1904
8 22.7601 22.7605 22.7645 22.8040 23.1710 24.1147 24.3281
9 25.9037 25.9041 25.9075 25.9422 26.2698 27.2264 27.4660
10 29.0469 29.0472 29.0503 29.0812 29.3767 30.3387 30.6040

With the eigenvalues, the general solution to the differential equation for v(r,t) is the sum
of all the individual eigenvalue solutions. The solution for u(r,t) is the sum of v(r,t) and
u .

Cm e m R m hR / k
2
u (r , t ) m t J 0 ( m r ) u [96]
m 1

The constants in this equation are found by satisfying the initial condition that u(r,0) =
u0(r).

u ( r ,0) u u0 (r ) u C m J 0 ( m r ) m R m hR / k [97]
m 1

We use the orthogonality condition for the Bessel functions to compute the constants.
Multiplying equation [97] by rJ0(nr) and integrating from r = 0 to r = R gives the following
result, using Matlab to get the integral in the denominator. Note that this integral has an

Jacaranda (Engineering) 3333 Mail Code Phone: 818.677.6448

E-mail: lcaretto@csun.edu 8348 Fax: 818.677.7062
additional term, [J0(mR)]2, that is not present in equation [46]. In that equation, the
eigenvalue equation, J0(mR) = 0, sets this term to zero.
R R

rJ 0 ( m r ) u0 (r ) u dr rJ 0 ( m r ) u0 (r ) u dr
0 0
Cm [98]
J
R 2
R
r J 0 ( m r ) 0 ( m R ) J 1 ( m R )
2 2 2
dr
2
0

As before, we consider the simplest case where u 0(r) is a constant equal to U0; this
gives the following result for Cm.

R rR
rJ1 ( m r )
rJ 0 ( mr )U 0 u dr 2U 0 u

m

0 r 0
Cm
J ( J 0 ( m R ) J 1 ( m R ) 2
2
R R2 2
m R ) J1 ( m R )
2 2
2
0 [99]
RJ1 ( m R)
2U 0 u
m 2U 0 u J1 ( m R) 2U 0 u J1 ( m )

R 2
J (
0 m R) 2
J 1 ( m R ) 2


m R J 0 ( m R ) J 1 ( m R )
2 2


m J 0 ( m ) 2 J1 ( m ) 2
Substituting this result for Cm into equation [96] for u(r,t) gives the solution for the
constant initial condition shown below after some algebra.

t

2U 0 u J1 ( m ) 2m
r
u (r , t ) J ( ) 2 J ( ) 2 e R2 J 0 m u
R
[100]
m 1 m 0 m 1 m

This equation can be rearranged into the following form.

2 t
u (r , t ) u 2 J1 ( m ) m
R 2 J r
e

[101]
m 1 m J 0 ( m ) J1 ( m )
0 m
U 0 u 2 2
R

Equation [101] shows the explicit dependence of a dimensionless potential difference

u (r , t ) u
on a dimensionless radial distance, r/R, and a dimensionless time, t/R2. In
U 0 u
addition, this solution depends on the ratio hR/k since the eigenvalues, m depend on
this parameter. The results of calculations for hR/k = 1 and hR/k = 10 are shown in the
Figures 8 and 9 as plots of dimensionless potential versus dimensionless radius with
dimensionless time, tau = t/R2 as a parameter. The comparison between these two
figures is similar to the comparison between Figures 6 and 7. Comparing Figures 6 and
7 showed that a higher values of the parameter hL/k give higher gradients at x = L and
shorter times for the values of u at any x location to get closer to u . In Figure 9, with a

Jacaranda (Engineering) 3333 Mail Code Phone: 818.677.6448

E-mail: lcaretto@csun.edu 8348 Fax: 818.677.7062
higher value of hR/k than Figure 8, the gradients at r = R are steeper than in Figure 8
and the times required to reach a given dimensionless u value are shorter.

Spherical geometry

In a spherical coordinate system, 0 r the one-dimensional diffusion equation (no

angular dependence) has the following form.

u 1 2 u
r [102]
t r 2 r r

The most general initial and boundary conditions for the radial problem are

u(r,0) = u0(r) u/r|r=0,t = 0 u(R,t) = uR(t) [103]

We will first consider the conditions where u R is a constant. We will later see that we will
have to use an expansion of Sturm-Liouville eigenfunctions to match the initial
condition. Anticipating this, we will obtain a problem with zero boundary conditions by
splitting u(r,t) into two parts, as we did in equation [37] for cylindrical geometry.

u(r,t) = v(r,t) + uR [104]

Jacaranda (Engineering) 3333 Mail Code Phone: 818.677.6448

E-mail: lcaretto@csun.edu 8348 Fax: 818.677.7062
The new function, v(r,t) satisfies the differential equation in [102] and has v/r = 0 at r =
0 and v = 0 at r = R. Since uR is a constant, the proposed solution in [104] will satisfy
the original differential equation and the boundary conditions from equations [102] and
[103]. As before, we seek a solution for v(r,t) using separation of variables. We
postulate a solution that is the product of two functions, T(t) a function of time only and
P(r) a function of the radial coordinate, r, only. With this assumption, our solution
becomes.

v(r,t) = P(r)T(t) [105]

We substitute equation [105] for u in equation [102]. Since P(r) is a function of r only
and T(t) is a function of t only, we obtain the following result.

u P ( r )T (t ) T (t )
P(r )
t t t
[106]
1 2 u 1 2 P (r )T (t ) 1 2 P ( r )
2 r 2 r T (t ) 2 r
r r r r r r r r r

If we divide the final equation through by the product P(r)T(t), we obtain the following
result.

Jacaranda (Engineering) 3333 Mail Code Phone: 818.677.6448

E-mail: lcaretto@csun.edu 8348 Fax: 818.677.7062
 1 1 T (t ) 1 1 2 P ( r )
r 2 [107]
T (t ) t P (r ) r 2 r r

The left hand side of equation [107] is a function of t only; the right hand side is a
function of r only. The only way that this can be correct is if both sides equal a constant.
As before, we choose the constant to be equal to2. This gives us two ordinary
differential equations to solve.

dT (t )
2T (t ) has the general solution T (t ) Ce t .
2
The first equation becomes
dt
The second ordinary differential equation in [107] can be written as follows.

d 2 dP (r )
r 2 r 2 P (r ) 0 . [108]
dr dr

We can show that the following proposed solution satisfies this differential equation.

sin r cos r
P (r ) A B . [109]
r r

Differentiating this proposed solution gives.

dP ( r ) cos r sin r sin r cos r

A B A 2 B 2 . [110]
dr r r r r

Substituting this first derivative of the proposed solution into the first term in the
differential equation in [108] gives.

d 2 dP (r ) d 2 cos r sin r sin r cos r

r r A B A 2 B 2
dr dr dr r r r r
d
Ar cos r Br sin r A sin r B cos r . [111]
dr
A cos r 2 Ar sin r B sin r 2 Br cos r A cos r B sin r
2 Ar sin r 2 Br cos r

Substituting this result and the proposed solution into equation [108] shows that the two
terms in the differential equation sum to zero; thus, the proposed solution satisfies the
differential equation.

d 2 dP(r )
r 2 rA sin r 2 rB cos r
dr dr
sin r cos r . [112]
2 2 2 2 2 2
r P (r ) r A B rA sin r rB cos r
r r

Jacaranda (Engineering) 3333 Mail Code Phone: 818.677.6448

E-mail: lcaretto@csun.edu 8348 Fax: 818.677.7062
At r = 0, the cos(r)/r term in the solution becomes infinite. In order to retain a finite
solution at the center of the sphere we must set B = 0. The sin(r)/r term becomes zero
at the center of the sphere, providing a finite term in the solution. With B = 0, the
product solution, v(r,t) = P(r)T(t) becomes.
2 sin r 2 sin r
v (r , t ) T (t ) P (r ) ACe t
De t . [113]
r r

We now have to satisfy the condition that v(r,t) = 0 at r = R,

2 sin R
v ( R, t ) T (t ) P ( R ) 0 De t . [114]
R

This gives the expected eigenvalue solution that R must be an integral multiple of .

n
R n or n [115]
R

The general solution for v(r,t) is a sum of all eigenvalue solutions, each multiplied by a
different constant, Dn. From equation [104], we know that the solution we want for u(r,t)
= v(r,t) + uR. Thus, our general solution for u(r,t) is written as follows

sin( n r ) n
Dn e t
2
u (r , t ) n uR with n [116]
n 1
r R

We still have to satisfy the initial condition that u(r,0) = u 0(r). Substituting this condition
for t = 0 into equation [116] gives

Dn nr
u 0 ( r ) u ( r ,0) r
sin
R
uR [117]
n 1

As usual, the eigenfunctions (1/r)sin(nr/R) form a complete set of orthogonal

eigenfunctions in the interval 0 r R; however, the Sturm-Liouville problem defined in
equation [108] has a weighting factor of r 2 that must be used in the orthogonality
integral. If we multiply both sides of equation [117] by r 2(1/r)sin(mr/R)dr, where m is
another integer, and integrate from zero to R, we get the following result. The steps
used here are similar to those in equation [14].

mr nr mr
R sin R sin sin
R dr R R dr
u0 (r ) u R r
2 2
Dn r
r n 1
r r
0 0 [118]
R R
mr nr 2 mr
Dn sin R
sin
R
dr Dm sin

R
dr
n 1 0 0

Jacaranda (Engineering) 3333 Mail Code Phone: 818.677.6448

E-mail: lcaretto@csun.edu 8348 Fax: 818.677.7062
Solving for Dm and evaluating the final integral gives the following result, which is similar
to equation [15].
R
mr
r u0 (r ) u R sin R
dr
2
R
mr
Dm 0
R
r u 0 ( r ) u R sin
dr [119]
2 mr R R
sin
R
dr 0

As an example, consider the simplest case where u0(r) = U0, a constant. In this case we
find Dm as follows.

2U 0 u R
R R
2 mr mr
Dm r u0 (r ) u R sin
dr r sin dr
R
0
R R R 0
2U 0 u R R
2 R
mr mr mr
sin cos
R m R R R 0
2 RU 0 u R 2 RU 0 u R cos m
sin m sin( 0) m cos m 0
m 2 m
[120]

Since cos(m) is 1 when m is even and 1 when m is odd, we can write the final result
for Dm as shown below.

2 RU 0
m m odd
2 RU 0 u R 1 m 1

Dm
m 2 RU 0
m m even
[121]

Substituting this result into equation [116] gives the following solution to the diffusion
equation when u0(x) = U0, a constant, and the boundary temperatures are zero.

2 RU 0 u R 1 m 1 2n t sin( n r ) n
u (r , t )
m
e
r
uR with n
R
[122]
n 1

We can rearrange this equation to introduce nr in the summation and create the
following dimensionless form.

n 2 t
u (r , t ) u R
1 n 1 e 2n t sin( r ) 2 1 n 1 R e nr
U 0 uR
2
nr
n
nr R2 sin(
R
) [123]
n 1 n 1

The dimensionless ratio on the left of this equation depends on the dimensionless
radius r/R and the dimensionless time, expressed as a Fourier number, t/R2.

Jacaranda (Engineering) 3333 Mail Code Phone: 818.677.6448

E-mail: lcaretto@csun.edu 8348 Fax: 818.677.7062
Spherical geometry with mixed boundary condition

We still have the differential equation in [102], but we now have the following boundary
conditions:

u(r,0) = u0(r) u/r|r=0,t = 0 ku/r|r=R,t = h(ur=R u) [124]

In order to solve this problem, we will first define a new variable, v(r,t) as we did in
equation [91], for the mixed boundary condition for the cylinder, to give a zero boundary
condition at r = R.

v(r,t) = u(r,t) u [125]

The boundary conditions for this new variable, v(r,t) are shown below.

v(r,0) = u0(r) u v/r|r=0,t = 0 kv/r|r=R,t = hvr=R [126]

The general solution for v(r,t) is will be the same as the solution to the problem with
fixed temperature given in equation [113]. This equation is copied below.
2 sin r 2 sin r
v (r , t ) T (t ) P (r ) ACe t
De t . [113]
r r

In order to satisfy the boundary condition at r = R, kv/r|r=R,t = hvr=R, with equation

[113], we have to satisfy the following equation.

v 2 cos R sin R sin R

kv hvr R 0 De t k h . [127]
r r R R R2 R

Dividing the final equation by (De-2tk/R2)sin R gives.

cos R hR R hR hR
R 1 1 0 R 1 tan R [128]
sin R k tan R k k

The roots of this equation are the eigenvalues of the solution to the present problem.
These eigenvalues will depend on the dimensionless group hR/k. We can compare this
equation with equation [56], which gives the eigenvalues for the rectangular slab with a
convection boundary condition. The roots of equation [56] are the values of that
satisfy = hL/k cot(). If we define = R in equation [128], the roots of that equation
are the values of that satisfy = (1 hR/k) tan(). For hR/k = 1, this equation can be
written as cot() = 0. Thus, the eigenvalues for hR/k = 1 are simply the zeros of the
cotangent, (2n 1)/2, where n is an integer greater than or equal to one. 10 For other
values of hR/k, we can write the eigenvalue equation as /(1 hR/k) = tan(). For hR/k
< 1, the roots of the equation [128] will be the intersection of the tangent curve with a

10
Note that = 0 is not a root of cot() = 0. As 0, cot() 1.

Jacaranda (Engineering) 3333 Mail Code Phone: 818.677.6448

E-mail: lcaretto@csun.edu 8348 Fax: 818.677.7062
line, starting at zero and having a positive slope. This will yield eigenvalues as shown in
Figure 10. The straight line for values of hR/k greater than one will have a negative
slope as shown in Figure 11.

Regardless of the behavior of the eigenvalues, the resulting solution will still have the
same form as equation [116]; however, this solution has u in place of uR, and the
eigenvalues are different.

sin( n r ) nR
Dne t
2
u (r , t ) n u with tan n R
r hR [129]
n 1 1
k

We still have to satisfy the initial condition that u(r,0) = u 0(r). Substituting this condition
for t = 0 into equation [129] gives

Dn
u0 ( r ) u ( r ,0) r
sin n r u [130]
n 1

As in equation [117], the eigenfunctions (1/r)sin(nr) form a complete set of orthogonal

eigenfunctions in the interval 0 r R, with a weighting factor of r 2 that must be used in
the orthogonality integral.

Jacaranda (Engineering) 3333 Mail Code Phone: 818.677.6448

E-mail: lcaretto@csun.edu 8348 Fax: 818.677.7062
Jacaranda (Engineering) 3333 Mail Code Phone: 818.677.6448
E-mail: lcaretto@csun.edu 8348 Fax: 818.677.7062
Repeating the process used in obtaining equation [118] multiplying both sides of
equation [130] by r2(1/r)sin(mr)dr, where m is another eigenvalue, and integrating from
zero to R gives the following result.
R R
2 sin m r sin m r sin n r
u0 ( r ) u r dr Dn r 2 dr
r r r
0 0 n 1
R R [131]
Dn sin m r sin n rdr Dm sin 2
m rdr
n 1 0 0

Solving for Dm and evaluating the final integral gives the following result, which is similar
to equation [119].
R R
2 sin m r sin m r
u0 ( r ) u r u0 ( r ) u r
2
dr dr
r r
Dm 0 0
R R sin m R cos m R [132]

sin
2
m rdr 2 2 m
0

Jacaranda (Engineering) 3333 Mail Code Phone: 818.677.6448

E-mail: lcaretto@csun.edu 8348 Fax: 818.677.7062
The integral of sin2(mr) was found by symbolic integration using Matlab. This integral
has an additional term that was not present in equation [119], where mR = m which
gives sin(mR) = 0.

As usual we will apply this equation to the simplest case where u 0(r) = U0, a constant. In
this case we find the numerator in the integral for Dm as follows.
R R
R sin m R cos m R
Dm r u0 (r ) u sin m r dr U 0 u r sin m r dr

2 2 m 0 0
2
1
U 0 u sin m r m r cos m r 0R
m
U 0 u
2
sin m R sin(0) m R cos m R 0 U 0 2 u sin m R m R cos m R
m m [133]

Solving for Dm and performing some algebraic rearrangement gives the final result for
Dm.

sin m R m R cos m R
Dm 2U 0 u
R2m m sin m R cos m R [134]

Substituting this expression for Dm into equation [129] gives the following solution for
u(r,t) in this case.

sin R n R cos n R sin( n r )
u ( r , t ) 2 U 0 u R2 n
2
e n t u [135]
n 1 n n sin n R cos n R r

We can obtain a dimensionless form by introducing m = nR and performing some

algebraic rearrangement.

n r
2n
t sin
u (r , t ) u sin n n cos n R
U 0 u
2
n
sin n cos n
e R2
n r
[136]
n 1
R

u ( r , t ) u
This equation shows that the dimensionless ratio, , depends explicitly on the
U 0 u
dimensionless radius, r/R and the dimensionless time (Fourier number) t/R2. In
addition, there is an implicit dependence on the dimensionless group hR/k, known as
the Biot number, through the eigenvalue n. A plot of equation [136] for hR/k = 1 is
shown in Figure 12.

Jacaranda (Engineering) 3333 Mail Code Phone: 818.677.6448

E-mail: lcaretto@csun.edu 8348 Fax: 818.677.7062
Jacaranda (Engineering) 3333 Mail Code Phone: 818.677.6448
E-mail: lcaretto@csun.edu 8348 Fax: 818.677.7062