Computing Heat Flow

Simulation of physics frequently involves fields whose behavior is defined by differential equations. Space is divided into small zones and field values are computed in these zones for a sequence of successive times. The zones are organized into a mesh which determines adjacency.

This practice was originally called difference equations emphasizing the connection to the “underlying” differential equations. More recently this practice is called finite element analysis and sometimes described as an approximation to “integral equations” that themselves are nearly equivalent to differential equations.

The first computing plan usually tried in such attempts is called explicit. In an explicit calculation field values are somehow known for all of the zones for some particular time and the corresponding values for the next time are computed from those of the current time in the current and immediately adjacent zones. This means that cause and effect can, within the simulation, travel thru the zones no faster than mesh velocity—one zone per time step.

I will talk here first about the one dimensional heat diffusion equation for it illustrates an interesting problem and several interesting general techniques for doing such simulations. Higher dimensional meshes are discussed too.

Imagine a rod where we are interested in the temperatures of equally spaced points along the rod and how they change over time. Suppose that heat flowing along the rod is the significant factor in knowing what the temperatures will be. Here is an interactive presentation of the heat equation, and another. Hear this one.

Let ϕkj denote the temperature in zone j at time tk. We can then compute directly from the equations:
k+1j − ϕkj]/Δt = [ejkj−1 − ϕkj) + fjkj+1 − ϕkj)]/Δx2
(This comes directly from ∂ϕ/∂t = e(x)∂2ϕ/∂x2)
The expressions in the parenthesis are temperature difference between adjacent zones and those divided by Δx are the temperature gradients which determine the rate at which heat is flowing to point k. The values ej and fj are conductivity of the bar and may in some situations depend on ϕ. These equations determine how fast each point is heating up, which is just the term on the left of the equation. I have assumed heat capacity to be one to simplify equations. (Alternatively varying heat capacity can be coded in the e’s and f’s.)

These equations can easily be put into the form
ϕk+1j = ajϕkj−1 + bjϕkj + cjϕkj+1

The a’s, b’s and c’s depend on Δt, Δx, the e’s and f’s. Choosing a large Δt would seem to get our answer sooner but our equations, as they stand, say that the heat flows during all of a time interval at the rate determined by the gradients at the beginning of the interval whereas in fact the temperature gradients may significantly decrease during the interval. For Δt too large the temperature changes will overshoot and oscillate wildly and soon lose all fidelity to the original problem and differential equation. The heat thus tries to flow too fast and this is a truncation error which we must deal with.

The following implicit equations are much more accurate.
k+1j − ϕkj]/Δt = [ejk+½j−1 − ϕk+½j) + fjk+½j+1 − ϕk+½j)]/Δx2
where ϕk+½j = (ϕkj + ϕk+1j)/2.
This says that the flow rate is determined by the temperature differences averaged over the time interval. Actually it is not quite that good; it uses the average of the beginning and ending gradients. These latter equations are called second order for the error is proportional to Δt2. The explicit equations are first order.

But it is not now immediately clear how to compute the k+1 values from the known k values. These equations can be viewed as a set of linear equations, one per zone. With only local zone calculations the equations can be put in the form
Ajϕk+1j−1 + Bjϕk+1j + Cjϕk+1j+1 = Dj for 0 < j < 100.
and now we solve for the ϕs. We establish boundary conditions by saying that:
ϕk0 = ϕk100 = 0 for k > 0
and that ϕ0j for 0 < j < 100 are provided as initial data somehow. At each time step we now have 99 unknowns.

This may seem daunting but the attentive first year algebra student who can program and has understood the first lecture on solving simultaneous equations is likely to come up with the best solution. Those who have also studied linear transformations and matrices may do much worse.

We are taught in first year algebra to eliminate unknowns from a set of equations, one at a time. The natural and strategic first candidate to eliminate is ϕ11. We take equation A1ϕ10 + B1ϕ11 + C1ϕ12 = D1 and solve for ϕ11 and get:
ϕ11 = (D1−A1ϕ10−C1ϕ12)/B1.

Note that ϕ10 is known from the boundary conditions. Of the remaining equations only the next, with j=2, involves ϕ11. We solve that equation for ϕ12 and each time there is just one equation into which to substitute the new expression. We have luckily avoided the quadratic explosion that accompanies this elimination procedure in the general case.

When we have solved the last equation for ϕ199 there are no unknowns in the resulting expression (ϕ1100 is known from boundary conditions) from which we have a concrete value for ϕ199. The previous equation which expressed ϕ198 in terms of ϕ199 now immediately provides ϕ198. This pattern continues in a backwards sweep thru the mesh until we have ϕ1j for all j.

This method was old hat in 1955 when I learned of it by word of mouth. It is numerically well behaved for physics problems such as these. It effectively moves the temperature signal several zones in one time step. Here is another more complete description of this calculation.

If our original problem had been for a metal loop instead of a metal rod we would not had a natural initial unknown to eliminate. There would have remained a lingering unknown just to the left of our initially eliminated unknown. This would have been resolved when we came around to it from the other side. It turns out that topology rules in these matters. We were driven to simplicity towards a forwards sweep followed by a backwards sweep. The Naïve solution was indeed strategic! There is a way to understand these problems abstracted from the equations. Here is a way to understand the accuracy and stability of this scheme.


It is constructive to consider eliminating unknowns in other orders. To first approximation this requires doing somewhat more arithmetic, but opportunities arise to deploy extra CPU’s or functional units. Imagine eliminating some interior point first. Without carrying out the algebra here you can probably see that we are left with an identical problem except for one fewer mesh points, around which some of the given numbers have changed. We can eliminate every third point with no interference and this can be done in parallel. Soon we have on or two equations which are quickly solved and we reverse the logic and find values for all the unknowns.