The Heat Equation, explained

Part 1: Derivation and examples

A partial differential equation is an equation that relates a function of more than one variable to its partial derivatives. To introduce PDEs, we’re going to solve a simple problem: modelling the temperature in a thin metal bar as a function of position and time. Along the way, we will derive the one-dimensional heat equation from physical principles and solve it for some simple conditions:

In this equation, the temperature T is a function of position x and time t, and k, ρ, and c are, respectively, the thermal conductivity, density, and specific heat capacity of the metal, and k/ρc is called the diffusivity.

Physical motivation

We would like to study how heat will distribute itself over time through a long metal bar of length L. One end of the bar is assumed to be at x=0 and the other is at x=L. The bar is much longer than it is thick so we can treat the distribution of heat as a function of just x and t. Assuming that the bar’s specific heat capacity is known, we will know how heat is distributing itself if we can find a function for the temperature T(x,t).

The bar is insulated along its length so that it can gain and lose heat only through its ends. This means that the temperature distribution depends only on these three things:

• The initial temperature distribution T(x,0). This is called the initial condition.
• The temperature at the ends of the bar, T(0,t) and T(L,t). These are called the boundary conditions.
• A law governing the rules for the transfer of heat from point to another within the body. The heat equation is a mathematical representation of such a physical law.

A problem that proposes to solve a partial differential equation for a particular set of initial and boundary conditions is called, fittingly enough, an initial boundary value problem, or IBVP.

To keep things simple so that we can focus on the big picture, in this article we will solve the IBVP for the heat equation with T(0,t)=T(L,t)=0°C. These are called homogeneous boundary conditions.

Deriving the heat equation

The heat equation can be derived from conservation of energy: the time rate of change of the heat stored at a point on the bar is equal to the net flow of heat into that point. This process clearly obeys the continuity equation. If Q is the heat at each point and V is the vector field giving the flow of the heat, then:

According to the Second Law of Thermodynamics, if two identical bodies are brought into thermal contact and one is hotter than the other, then heat must flow from the hotter body to the colder one at a rate proportional to the temperature difference. Therefore V is proportional to the negative gradient of the temperature, so V=-k∇T where k is the thermal conductivity of the metal. In one dimension this reduces to V=(-k∂T/∂x)x where x is the unit vector in the +x-direction. Furthermore, Q=ρcT, so we get the heat equation by plugging in these expressions for V and Q:

Solving the heat equation

Before we can go any further, we need to show that there must exist a unique solution to the heat equation for any physically meaningful initial and boundary conditions. The formal proof of this is beyond the scope of this article, so instead we’ll use an empirical argument. The laws of thermodynamics tell us that no matter what the temperature distribution of the bar is initially, the system mustundergo a process that brings the bar to thermal equilibrium, and we showed in the previous section that this process must obey the heat equation, so solutions to the heat equation exist for physically meaningful initial and boundary conditions. Furthermore, one of the basic assumptions of classical physics is that identical experimental conditions must lead to identical results, so the particular way in which the bar progresses to thermal equilibrium is uniquely specified by the initial and boundary conditions.

This means that if f(x,t) and g(x,t) are two different functions that satisfy the same IBVP for the heat equation, then f and g have the same form. Furthermore the heat equation is linear so if f and g are solutions and α and β are any real numbers, then αf+βg is also a solution. So we can conclude that the solution is going to be a linear combination of functions of the same form.

Consider the following function, which we could have guessed by inspection or by trial and error:

Where n is a positive integer greater than zero. This function satisfies the heat equation:

This function also satisfies the boundary conditions since sin(0)=sin(nπ)=0. Therefore the general solution is:

The problem will be solved if we can find the coefficients Aₙ such that this general solution satisfies the initial condition. That is, we need to find Aₙ such that:

This is called a Fourier sine series expansion for the initial conditions. The coefficients Aₙ called the Fourier coefficients.

Computing the Fourier coefficients.

The initial condition T(x,0) is a piecewise continuous function on the interval [0,L] that is zero at the boundaries. It turns out that the set of functions with these properties is a vector space under addition and scalar multiplication. We will call this vector space 𝕊ᴸ.

This vector space be equipped with an inner product. For f,g∈𝕊ᴸ, a possible inner product is:

An inner product generalizes the idea of the dot product. We can find the components of a geometric vector by projecting it onto the axis using the dot product with the unit vectors, which form a basis for ℝⁿ. In the same way, if we can find a basis for 𝕊ᴸ then we can project any f∈𝕊ᴸ onto the basis functions in order to represent f as a linear combination of basis functions.

For integers m,n>0, functions of the form sin(nπx/L) are orthonormal:

So we can represent any function f∈𝕊ᴸ as a linear combination of functions from the basis set:

And the coefficients of that linear combination will be given by the Euler integral:

To demonstrate, let’s find the Fourier coefficients of a unit sawtooth pulse:

Clearly, saw(x)∈𝕊¹ so:

And the Fourier coefficients are given by:

So the Fourier series expansion for the sawtooth wave is:

This animation shows how the Fourier series approaches the sawtooth as the number of sine terms in the sum increases.

The error near x=1 is called the Gibbs phenomenon. The Gibbs phenomenon is an unavoidable error that causes the Fourier series of a discontinuous function to overestimate the function’s value at a discontinuity by about 9%. The Gibbs phenomenon can never be completely eliminated, but as the number of terms in the Fourier series approaches infinity, the error converges to being localized entirely at the point of discontinuity. So for example if we included an infinite number of terms in the Fourier series expansion for the sawtooth, we would find that the series would be exactly equal to x for 0≤x<1 and at x=1 the series would have a value of about 1.09.

What this tells us is that solving the homogeneous IBVP for the heat equation amounts to using the Euler integral to find the Fourier coefficients:

Example: Bar initially at uniform temperature

An insulated, meter-long bar with k/ρc=0.1m²/s (unrealistic but easier for graphing purposes) is initially 100°C throughout its length before cooling elements set to 0°C are clamped to its ends at t=0. The initial and boundary conditions are:

The Fourier coefficients are:

Let’s verify that this gives a Fourier series that matches the initial and boundary conditions:

So the solution is:

Now let’s graph the solution:

And as a three-dimensional plot:

Example: A temperature “spike”

Now suppose that the bar is initially 0°C everywhere along its length except for the middle 10 centimeters where the temperature is 100°C. This time the diffusivity is 0.0075m²/s.

The easiest way to compute the Fourier coefficients is to put them into a more general form:

So the solution is:

Although the bar will eventually reach an equilibrium state with uniform temperature 0°C, after five seconds the temperature decreases very slowly so I decided to end the animation there. Here’s the 3D plot:

Example: Random heat distribution

This time, the bar, which we will assume to be made of copper with diffusivity 1.11×10⁻⁴m²/s, has a random heat distribution that’s been measured at 10 centimeter intervals according to the following table:

In situations like this one, it’s best to find the Fourier coefficients by numerically integrating rather than trying to find a closed-form expression. I used MATLAB for this, and if you want to see the code then you can get it here on Pastebin.

Here’s an animation of the solution:

After 240 seconds, the temperature changes only very slowly so I decided to end the simulation there. The most interesting behavior happens in the first 60 seconds, so I only went that far for making the 3D plot:

Concluding remarks

So that about wraps this first part up. Now that you know how to solve the heat equation for the simplest cases, you’re ready to use the heat equation to analyze more interesting problems. I will cover some of these in an article that will be out in the coming days.

Copyright

The images, animations, and code in this article are my own original work. I don’t care if you want to use them elsewhere as long as you attribute me. The MATLAB code that I linked should run without issue on any reasonably modern machine but obviously I cannot promise that it will work flawlessly for everyone.