version of 11 November 1997

VIII.  The mathematics of hot rods.†

In many ways the simplest of our three model partial differential equations is the heat equation, also known as the diffusion equation,
	ut = k —2u							(HE-3D)
or, specializing for now to one dimension,

	ut = k uxx.							(HE)

A physical derivation of this equation is given in an appendix; the constant 
	k = s/r ks,
where s is the thermal conductivity, r is the density, and ks is the specific heat capacity of a homogeneous substance.  Ordinary substances have values of k ranging from about 5 to 9000 cm2/hr (see the table at the end of the chapter).

We shall solve the heat equation with the method of separation of variables, as we did for the wave equation.  The solutions will be quite similar, and you may find it instructive to compare with what we do in chapter VI.  We begin with the one-space-dimensional heat equation (HE), which can be written in terms of a heat operator

	Hk [f]  := b( f(1,k)  f(∂,∂t) - f(∂2,∂x2)  )f(t,x),

so that (HE) is the equation for the kernel of this operator.  This equation describes the temperature u(t,x) in a long, thin rod along the x-axis, if the rod is wrapped along its length with a perfect insulator.  The heat equation is a linear, homogeneous equation, so, just as for the wave equation, the superposition principle holds; if we can find some solutions u1, ..., uk, then any linear combination of these solutions is again a solution.

With the method of separation of variables, we try to find solutions of a special form, 
		u = T(t) X(x). 
As before, boundary and initial conditions will be important, and arise from the physical situation we wish to study.  The end points of the rod may be in thermal contact with ice water, which will ensure that we have homogeneous Dirichlet boundary conditions
(DBC)
	u(t,a) = u(t,b) = 0.
(Actually, stirred wet crushed ice would guarantee these boundary conditions in the laboratory better than ice water.)

The fact that these boundary conditions are homogeneous (something equals 0) is going to be useful in our analysis, so you may wonder if something goes wrong if we use a different temperature scale or hold the ends at a different fixed temperature T different from 0.  Fortunately, we only need to redefine u by a constant to have these same boundary conditions.  Since Hk [T] = 0, if 
	v(t,x) := u(t,x) - T, 						(8.1)
then
	Hk[v] = Hk[u] = 0.
By this simple change, the problem with some other fixed temperature at the ends is converted into out standard problem with homogeneous boundary conditions of the Dirichlet type.

Neumann boundary conditions (NBC),
	ux(t,a) = ux(t,b) = 0,
correspond physically to the situation when the ends of the rod are also insulated, since the heat flow through either end of the rod is proportional to the temperature gradient at that end.  If there is no heat flow, the gradient must be 0.

We begin by considering the case of Dirichlet boundary conditions, and suppose that the rod has length 1, with a = 0 and b = 1.  We shall solve this problem with the method of separation of variables, by looking for special solutions analogous to normal modes for the wave equation.  Ignoring the initial conditions for the time being, we begin by seeking product solutions,
	u(t,x) = T(t) X(x)
to the heat equation (HE), which also satisfy the homogeneous boundary conditions.

The efficient way to proceed is to divide the wave equation by the product solution:

	f(1,T(t) X(x)) Hk [T(t) X(x)] = f(T',k T) - f(X'',X) .

If u solves the heat equation, this expression must be zero. Thus 
           T'/(k T)  = X''/X.
Since the left side is independent of x and the right side is independent of t, both sides must in fact a constant.  We don't know its value yet, so let us call it - m.  (It will turn out to be negative.  

Thus   
	-X''(x) =  m X(x)						(8.2)
	
This is precisely the same eigenvalue equation as in chapter VI, for the wave equation with (DBC).  We know that the eigenfunctions and eigenvalues are:

	X(x) := sin(n p x),
	mn = (np)2.

The other part of the partial differential equation which has been separated differs from what we got with the wave equation:

              T'(t) = -  mn k T(t)						(8.3))
              
This is a simple first-order equation, the solutions of which are exponential functions,
	T(t) = C exp(- mn k t).
The eigenvalue problem (8.2) fixed the possibilities for mn, leaving us with normal modes of the form 
	cn exp(- n2p2 k t) sin(npx), n = 1, ...

With DBC at a and b, b-a = L, we would get:

	cn exp(- n2p2 k t) sin(np(x-a)/L), n = 1, ...;

the easiest way to see this is to change variables to x' = (x-a)/L.  In passing, we note that we could have set k=1 throughout the derivation of the normal modes, and then changed variables to t' = k t at the end.  

Since the heat equation is linear, we can make linear combinations of the solutions we have found and thereby produce new solutions. The most general solution to the heat equation that we get in this way is:

	∑n=1 cn exp(- n2p2 t) sin(npx).					(8.4)

At this stage it may not be completely clear that this is the absolutely most general solution, but the completeness theorem for the Fourier sine series will guarantee that this can indeed describe the solution to the heat equation with an arbitrary square-integrable initial condition.

The initial conditions appropriate to the heat equation are somewhat different from what we had for the wave equation, again because the time variable enters differently.  The heat equation is first-order in time, so we need only one initial condition, the initial temperature distribution in the rod, 
	u(t,x) = f(x).                                                             (IC)

There is no need to look at ut at t=0.  Physically, this would be inappropriate, because our experience teaches us that the temperature changes only in response to the initial temperature, not to some unmeasurable effect of the rate of change of the initial temperature.  If this were not the case, cookbooks would be nearly impossible.  Fortunately, the initial conditions to bake a cake are, typically, batter at room temperature, and if we have identical batter in identical pans at room temperature, we need to bake the cake for identical lengths of time.  No other physical initial condition is necessary.

Model Problem VIII.1.  Find the future temperature within a rod of length 1, k=1, insulated along its length, and with ends in thermal contact with ice water, if at time 0, the rod is at a uniform temperature u(0,x) = 50.

Solution (using software for calculations).  At time t=0, the general solution becomes a Fourier sine series,

	∑n=1 cn  sin(npx),

so the coefficients cn are determined as the Fourier sine  series coefficients of the initial function,  f(x) = 50.  We know the formula for the Fourier sine coefficients, (2.14), so we simply calculate the integral to get:

	cn = f(100,n p) (1 - (-1)n ). 

This is 0 when n is even and otherwise 200/(n p).  Our solution to the heat equation is thus:
	u(t,x) = ∑n odd (200/np) exp(- n2p2 t) sin(npx).
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

Suppose now that the rod is insulated not only along its length but also on its ends.  Then u(t,0) and u(t,L) are not specified, but instead we have:

          ∂u[t,0]/∂x = 0,      ∂u[t,L]/∂x = 0                         (NBC)

Model Problem VIII.2.  Set up the general solution of the insulated rod with Neumann BC., with the simplifications k = L = 1.

Solution.
This is very similar to the analysis in chapter VI, and the eigenfunctions are the basis functions of the Fourier cosine  series, 
	{1, cos(n p x)}.
The normal modes are 
	cn exp(- n2p2 t) cos(npx), n = 0, 1, ...
and the general solution is:

	∑n=0 cn exp(- n2p2 t) cos(npx),					(8.5)
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
Notice that after a long time, all of the terms of this series will become negligible except the one corresponding to n=0.  Remembering that c0 is nothing other than the average temperature, we see that if the rod is entirely insulated, then after a long time the temperature has settled down to equilibrium at the one temperature which is consistent with conservation of energy.  

Another way to see the same effect is as follows.  I will do the calculation in three dimensions, since it is not much harder than in one dimension.  In the derivation of the heat equation in the appendix, we reason that the rate of change of heat content is

If the object is insulated, then n◊—u = 0.  (This higher-dimensional Neumann condition is what simplifies to ux = 0 in one space dimension.)  It follows that the total heat content 
	
remains constant.  For a homogeneous material, this is just a multiple of the average temperature.

Something more subtle happens in the case of DBC, where energy is lost (or gained) through the ends of the rod.  In that case, the general solution (8.4) ensures that not only will the temperature settle down to an equilibrium u(t,x) Æ 0 as x Æ ∞, but it usually will do so with a particular temperature profile, for assuming that c1 ≠ 0, we find that 
	exp(p2t) u(t,x) = c1 + r(t,x),
where the remainder terms in 
	r(t,x) = ∑n=2 cn exp(- (n2-1)p2 t) sin(npx)
all vanish rapidly as tÆ∞.  As can be seen from the calculations in the Mathematica notebook, the temperature profile of the rod rapidly approaches that of a simple sine function, times a time-varying factor tending to 0. 

Model Problem VIII.3.  Suppose now that we have non-homogeneous boundary conditions, such as what we would get if the end of the rod at x=0 is held at temperature 100 and the end at x=1 at temperature 0.  Find the general solution to the heat equation with k=1.

Solution.
We can no longer convert this problem to the usual one by the simple trick (8.1).  A cleverer idea is to think about what should happen at equilibrium in this problem.  At equilibrium, ut = 0, so the temperature will depend on space alone according to 
	u(eq)xx = 0.								(8.6)
In this one-dimensional situation, this means that u(eq)(x) = c1 x + c2 - the general solution to (8.6).  With the given boundary conditions, we can easily solve for the constants, and find that 
	u(eq)(x) = 100 - 100 x.

Since the heat equation is linear we can adjust the solution to obtain our familiar boundary conditions by subtracting the equilibrium, letting
	v(t,x) := u(t,x) - u(eq)(x).

It is easy to see that
	Hk[u - u(eq)] = 0,
and the general solution is thus
	u(t,x) = 100 - 100 x + ∑n=1 cn exp(- n2p2 t) sin(npx)
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

The properties of a solution to the heat equation are quite different from those of the wave equation.  Instead of wave motion and a tolerance for discontinuity, we will find
1.  Equilibration - after a long time, the solution becomes independent of time.
2.  Intolerance for extremes and discontinuity.  Even if the solution starts off as quite a wild function f(x), it will instantaneously become a smooth function, which can be differentiated arbitrarily often.  

We investigated equilibration above, but have not yet seen how the heat equation enforces smoothness.  One aspect of this phenomenon can be seen from the general solution, say with DBC, (8.4).  We learned in chapters II-III that a function can behave rather terribly and still be the sum of a Fourier series, where the squares of the coefficients must tend to 0 as n Æ∞, owing to the Parseval formula (3.4).  Because of this, suppose that all we know about the coefficients is that they are bounded by some tremendous number B.  Now, if we could say that we are safe in neglecting the terms in a Fourier series beyond some finite point, n > N, then we would be left with a finite sum of sines or cosines, and that would clearly be differentiable as often as we would like.  We also learned in chapter V that we can safely differentiate Fourier series, provided that the differentiated series makes sense.

With this in mind, let's differentiate our general solution, say, 900 times with respect to x.  It is not hard to see from (8.4) that we get
	∑n=1 cn exp(- n2p2 t) (np)900 sin(npx),                               (8.7)
where all we know about the cn is that they are bounded by B.  (We would have some minuses or cosines if we differentiated a number of times not divisible by 4.)  If t=0, we have amplified all the coefficients, but if t > 0, then for large n, the exponential goes to zero faster than any power of n, and the series (8.7) converges uniformly.  This happens no matter how small t is, and shows that the 900-th x-derivative of u exists as a continuous function (a theorem from advanced calculus states that a uniform limit of continuous functions is continuous).

A careful analysis suggested by this argument proves that even if the initial condition is extremely irregular, a solution of the heat equation instantly becomes smooth.  A precise definition of what we mean by “smooth” is as follows.

Definition VIII.4.  A function f Œ C∞ if it can be differentiated arbitrarily often. 

Theorem VIII.5.  Let u solve the heat equation on an interval a < x < b, and t > 0, with initial condition u(0,x) = f(x) square integrable and either Dirichlet or Neumann boundary conditions.  Then for all t > 0, u(t,x) Œ C∞.  

Later we shall consider the heat equation on higher dimensional regions Ω.  The same theorem holds for all x in the interior of Ω.  Also, other standard kinds of boundary conditions can be allowed in this theorem.

Another important property of the heat equation is the maximum principle.  

Theorem VIII.6.  Let u solve the heat equation.  If we consider any finite space-time region, such as a rectangle a ≤ x ≤ b, c ≤ t ≤ d, then the greatest value of the temperature occurs on the boundary of the region, in this case when x = a or b or when t = c or d.  In other words, a hot spot cannot be spontaneously generated in the interior of the material.  The same is true of the minimum, a potential cold spot - the minimum always occurs on the boundary.

As with the smoothing theorem, the maximum principle also applies to higher-dimensional regions.

To prove this, let us first consider the heat equation with a heat sink, which absorbs energy.  The equation describing this would be 
	ut = uxx - Q(t,x),							(8.8)
where Q ≥ 0.  For now, suppose that Q is strictly positive, so Q ≥ e > 0 throughout the rod.  Suppose that there is a hot spot, i.e., a local maximum, at (t0,x0), in the interior of the space-time region.  By basic calculus, the space-time gradient of u must be 0 there (we know that the solution is differentiable), that is, ut(t0,x0) = ux(t0,x0) = 0.  If the point is really a maximum, then uxx(t0,x0) ≤ 0.  But this is a contradiction, since
      uxx = ut+ Q ≥ e > 0.

This shows that there can be no local maximum in the interior of the region, so the only possibility is for the maximum to be on the boundary.

Now suppose that there is no sink, so we have a solution of the heat equation with no sinks or sources, 
	ut = uxx, and that there is a point (t0,x0) where the temperature exceeds that on the boundary of a finite region.  To be specific, suppose that we replace u with the function 
	v(t,x) := u(t,x) + d exp(x),
where d is a small but positive quantity, so small that the maximum value of v is also attained in the interior - this is possible by the intermediate value theorem of calculus, which applies to continuous functions.  What equation does v solve?
	vt = ut, 
while
	vxx = uxx + d exp(x) = ut + d exp(x) = vt + d exp(x).
Solving for vt, we have
	vt = vxx - Q,
where Q =  d exp(x) is strictly positive.  We already considered this case, and found that it cannot happen.  The transformation shows that it cannot happen for the basic heat equation, either.
										QED
The statement about the minimum follows because -u solves the heat equation whenever u does.

Exercises VIII

VIII.1.  Solve the heat equation with DBC at 0 and 1 and:
                     f(x) = |x-1/2|.

VIII.2.  Solve the heat equation with NBC at 0 and 1 and:
                     f(x) = sin(px)

VIII.3.  Put DBC at x=0 and NBC at x=L.  Separate variables, solve the eigenvalue problem for X[x], and write down a general solution.

VIII.4.  Solve the heat equation with IC:
                     f(x) = 1
and
a)  DBC at x=0 and x=1
b)  NBC at x=0 and x=1
c)  DBC at x=0 and NBC at x=1

Plot the solutions and compare.

VIII.5.  Find some initial conditions for which a unit rod with DBC will not have an asymptotic temperature profile with spatial dependence sin(px), and use the computer to study the approach to equilibrium in this case.

VIII.6.  Solve the initial-boundary value problem
	ut = uxx
	u(t,0) = 50, u(t,2) = 100
	u(0,x) = 20 x

VIII.7.  Solve the initial-boundary value problem
	ut = uxx + 1
	u(t,0) = 100, u(t,2) = 100
	u(0,x) = 20 x
Hint:  subtract something off to eliminate the non-homogeneity in the PDE as well as in the BC.

VIII.8.  Use the ratio test and l'Hôpital's rule of advanced calculus to show that the series (8.7) converges uniformly.

VIII.9.  It is possible for a rod to be in equilibrium at a constant temperature, and in that case it attains its maximum, and only, temperature in the interior.  Explain why this state of affairs is consistent with the maximum principle and with its proof.  Prove the maximum principle with 2 or 3 spatial variables, using the same argument.

VIII.10.  Does the maximum principle hold for ut = uxx + Q(t,x), Q > 0?  Either give a proof or a counterexample.  Interpret physically.

Appendix.  Derivation of the heat equation.

Newton was the first to articulate some principles of heat flow, but it was Fourier who created the correct systematic theory of how heat flows through solids.  Inside a solid there is no convective transfer of heat energy and little radiative transfer, so temperature changes only by conduction, as the energy we now recognize as molecular kinetic energy flows from hotter regions to cooler regions.  The first basic principle of heat is 

1.  the heat energy contained in a material is proportional to the temperature, the density of the material, and a physical characteristic of the material called the specific heat capacity.  In mathematical terms,
							(8.a1)
For the other principles of heat transfer, let us do some experiments with the following materials:  a hot stove, some iron rods of different, relatively short  lengths and various widths, and various ceramic rods of different lengths and widths.  Since these will be thought experiments only, it will be safe to use a finger as the probe.  Putting your finger right on the stove will convince you that the energy transfer is proportional to the difference in temperature between your finger and the stove.  Using, if necessary, a different, undamaged finger, you will also find that the rate of heat transfer is inversely proportional to the length of an iron rod intervening between your finger to the stove  (fixing the cross-sectional area).  In other words, the rate of heat flow from one region to another is proportional to the temperature gradient between the two regions.  You will probably also agree that the rate of heat flow will be proportional to the area of the contact; for example, a short pin with one end on a hot stove and the other touching your hand is preferable to putting the palm of your hand on a frying pan.  Finally, a ceramic material on the stove being usually more pleasant to the touch than hot iron, we see that the rate of heat transfer depends on the material, as measured with a physical constant known as the heat conductivity  s.  The second basic principle is thus:
2.  The heat transfer through the boundary of a region is proportional to the heat conductivity, to the gradient of the temperature across the region, and to the area of contact, or, in mathematical terms, if the boundary of the region is ∂Ω, with outward normal vector n, then

If we differentiate eq. (8.a1) with respect to time and apply  Gauss's divergence theorem to (8.a2), we find that dQ/dt can be expressed in two ways as an integral:

Since the region Ω can be an arbitrary piece of the material under study, the integrands must be equal at almost every point.  If the material under study is a slab of a homogeneous substance, then r, ks, and s are independent of the position x, and we obtain the heat equation

where k = s/r ks.  The one-dimensional heat equation 
	ut = k uxx
would apply, for instance, to the case of a long, thin metal rod wrapped with insulation, since the temperature of any cross-section will be constant, due to the rapid equilibration to be expected over short distances.

Here are some approximate physical data, from sources including the CRC Handbook of Chemistry and Physics, 48th ed., Cleveland:  Chemical Rubber Company, 1967, E.S. Dana and W.E. Ford, Textbook of Mineralogy, New York: Wiley, 1964, and C.J. Giankoplis, Transport Processes and Unit Operations, Englewood Cliffs, NJ: Prentice-Hall, 1993.  These characteristics actually depend on such factors as temperature and pressure.  According to the Dulong-Petit law, the specific heat capacity of an element is roughly inversely proportional to its atomic weight, but there is no simple rule for solid compounds.  Thus these are not necessarily accurate. 
		s			r		ks			k
units		cal/(hr-cm-deg)	gm/cm3	cal/(gm-deg)		cm2/hr
glass	 	80			2.8		.2 (guess)		140
iron		6800			7.9		.1			8600
granite	300			2.7		.19 (orthoclase)	585
pork     	75  			1.1		0.6			113
water		5			1		1 (or .5 for ice)	5
wood 		350  			.7		.2 (guess)		2500
air		.002 			.0011		.25 (N2)		7