Derivation of the heat equation

The first type of physical situation we want to consider is how the temperature of an object changes over time. We have dealt with this idea previously for a well-mixed fluid at a single temperature with Newton’s Law of cooling. However, we will consider a physical object here, so the temperature will depend on position as well, giving rise to a PDE instead of the ODE we got with Newton’s Law.

Let’s start by considering an insulated wire of length $L$. With this, we mean that no heat can escape through the lateral sides of the wire, and the only way that heat can enter the system is through the two ends. We will let $u(t,x)$ be the temperature in the wire after $t$ seconds, and at the position $x$ meters from the left end of the wire (the units here are not important, they just need to be consistent). See Normally a reference to a previous figure goes here..

---

---

We want to come up with a way to find a differential equation that models how this temperature will change with time and location. The easiest way to do this is with the accumulation equation from modelfirst:section with thermal energy, since that is something that can be accumulated, while temperature really can’t. How does temperature relate to thermal energy? Physical properties tell us that $$\Delta E = c \rho V \delta T$$ where $\Delta E$ is the change (or input) of energy, $c$ is the specific heat of the material, $\rho$ is the mass density of the material, $V$ is the volume, and $\delta T$ is the change in temperature. This expression mainly comes from the fact that the specific heat is the energy required to raise one kilogram of the material by one degree Celsius.

Now, we need to build up the accumulation equation for energy. To do this, we will look at a small interior segment of this wire of length $\Delta x$ that stretches from a position $x$ to $x + \Delta x$. The accumulation equation tells us that $$\text {rate of change of energy} = \text {rate in from left} + \text {rate in from right} + \text {interior sources}$$

For interior sources, we will assume that there is some function $Q(t,x)$ that describes all of the sources of thermal energy per unit volume at time $t$ and position $x$, and our discussion of the relationship between temperature and thermal energy says that the rate of change of energy is $$\frac {\partial E}{\partial t} = c\rho (A\Delta x)\frac {\partial u}{\partial t}$$ since $u$ represents our temperature and $A\Delta x$ is the volume of this small piece of wire. Since the piece is “small”, we will assume that things like the energy change are constant along this piece, which works out using linearization since we are eventually going to limit $\Delta x$ to zero.

At this point, we have the equation $$c(x)\rho (x) (A\Delta x)\frac {\partial u}{\partial t} = \text {rate in from left} + \text {rate in from right} + Q(t,x)(A \Delta x)$$ where we have also added in the fact that $c$ and $\rho$, these material properties, could depend on $x$ if the wire is not uniform and made of the same material throughout. The last things we have to deal with are these rate terms. The normal terminology used for the rate of energy entering or exiting a region is flux of energy. Fourier’s law (same Fourier) says that if $\phi$ is the flux of energy, then $$\phi = -k \nabla u$$ where $u$ represents temperature and $k$ is the thermal conductivity of the material. The point of this expression for now is the fact that the thermal flux is proportional to the derivative of temperature. This means that in this particular case $$\text {flux in from left} = -k(x) A\frac {\partial u}{\partial x}(t, x)$$ and $$\text {flux in from right} = k(x+\Delta x) A \frac {\partial u}{\partial x}(t, x+\Delta x).$$ The minus sign is missing from the second term because the flux coming in from the right is moving in the negative direction, so it needs to pick up an additional minus sign. The $A$ here represents the cross-sectional area of the wire. Thus, we have the full equation

$$c(x)\rho (x) (A\Delta x)\frac {\partial u}{\partial t} = -k(x) A\frac {\partial u}{\partial x}(t, x) + k(x+\Delta x) A \frac {\partial u}{\partial x}(t, x+\Delta x) + Q(t,x)(A \Delta x).$$

Dividing both sides by $A \Delta x$ gives $$c(x)\rho (x) \frac {\partial u}{\partial t} = \frac {1}{\Delta x}\left ( k(x+\Delta x) \frac {\partial u}{\partial x}(t, x+\Delta x) - k(x) \frac {\partial u}{\partial x}(t, x) \right ) + Q(t,x).$$

Now, we want to take the limit as $\Delta x \rightarrow 0$. The one term that still involves $\Delta x$ looks like $$\frac { k(x+\Delta x) \frac {\partial u}{\partial x}(t, x+\Delta x) - k(x) \frac {\partial u}{\partial x}(t, x) }{\Delta x}$$ which, as $\Delta x \rightarrow 0$ is a partial derivative. Taking this limit, we get the equation $$c(x)\rho (x) \frac {\partial u}{\partial t} = \frac {\partial }{\partial x}\left (k(x) \frac {\partial u}{\partial x}\right ) + Q(t,x).$$

In most cases, we will assume that the material we use is uniform. That means that $c$, $\rho$ and $k$ are all independent of $x$. We can move the $k$ outside the derivative, then divide both sides by $c\rho$ to get the equation $$\frac {\partial u}{\partial t} = \alpha \frac {\partial ^2 u}{\partial x^2} + F(t,x)$$ where $$\alpha = \frac {k}{c\rho }$$ is the thermal diffusivity of the material and $$F(t,x) = \frac {Q(t,x)}{c\rho }$$ represents the thermal sources in the system. This is the heat equation, in particular, this is a non-homogeneous heat equation. The homogeneous version is $$\frac {\partial u}{\partial t} = \alpha \frac {\partial ^2 u}{\partial x^2}.$$

For an interpretation, this tells us that the change in heat at a specific point is proportional to the second derivative of the heat along the wire. This makes sense; if at a fixed $t$ the graph of the heat distribution has a maximum (the graph is concave down), then heat flows away from the maximum. And vice versa.

Boundary conditions for the heat equation

There are several main types of boundary conditions that can be used with the heat equation. In general, we will assume that either the ends of the wire are exposed and fixed at a constant temperature, or they are insulated. In the first case, the endpoint being fixed at a constant temperature means that $$u(t,0) = T_0$$ at the left endpoint and $$u(t,L) = T_L$$ at the right endpoint. For insulated ends, this means that thermal energy can not escape from the endpoint. Since the flux is proportional to the derivative of temperature, this means that $$\frac {\partial u}{\partial x}(t, 0) = 0$$ at the left endpoint and $$\frac {\partial u}{\partial x}(t, L) = 0$$ at the right endpoint. When trying to solve these problems, we want to have conditions where something is set to zero at the end point; that will make our lives much easier once we separate variables. These side conditions are said to be homogeneous (i.e., $u$ or a derivative of $u$ is set to zero).

It is possible to have more involved boundary conditions that involve both $u$ and the derivative. These are sometimes more physically natural, since they are similar to Newton’s law of cooling, but are much harder to solve. These may look something like $$\frac {\partial u}{\partial x}(t,0) - 4 u(t,0) = 0,$$ but we will generally avoid these types of conditions here.

We also need an initial condition—the temperature distribution at time $t=0$. That is,

$$\begin{equation*} u(x,0) = f(x) , \end{equation*}$$ for some known function $f(x)$. This initial condition is not a homogeneous side condition.

Multidimensional heat equation

The same ideas can be used to derive the heat equation in plates (two dimensions) or regions (three dimensions) as opposed to just one dimensional wires. The process is the same: take a small region and relate the change of energy to the flux through each side. In the two dimensional case, there are four sides to consider, two where the flux is related to $\frac {\partial u}{\partial x}$ and two where it is related to $\frac {\partial u}{\partial y}$. The differences end up working out the same way and with the same assumptions as before, this results in the heat equation in two dimensions looking like $$\frac {\partial u}{\partial t} = \alpha \left ( \frac {\partial ^2 u}{\partial x^2} + \frac {\partial ^2 u}{\partial y^2}\right ) + F(t,x, y).$$ With three spatial dimensions, there are six different flux terms, and the equation simplifies to $$\frac {\partial u}{\partial t} = \alpha \left ( \frac {\partial ^2 u}{\partial x^2} + \frac {\partial ^2 u}{\partial y^2} + \frac {\partial ^2 u}{\partial z^2} \right ) + F(t,x, y, z).$$

In each case, we end up with a term that is the sum of all of the different second partial derivatives in the spatial directions. This type of a term shows up very frequently, so we give it a name.

With all of these definitions, the homogeneous heat equation, in any number of dimensions, can be written as $$\frac {\partial u}{\partial t} = \alpha \Delta u.$$

Boundary conditions in the multidimensional case work similarly to the one-dimensional case. The temperature can either be fixed at the boundary, giving a value for the function $u$, or insulated, giving that the flux must be zero, or that $$\nabla u \cdot \vec {n} = 0$$ along the boundary of the domain. For a rectangular domain, this will look like $$\frac {\partial u}{\partial x} = 0 \quad \text { or } \quad \frac {\partial u}{\partial y} = 0.$$

Solution by separation of variables

Now, we want to employ techniques from secondpde:section to try to solve the heat equation for a variety of initial and boundary conditions.

Solution: We use the method of separation of variables and guess $u(x,t) = X(x)T(t)$. We will try to make this guess satisfy the differential equation, $u_t = k u_{xx}$, and the homogeneous side conditions, $u(0,t) = 0$ and $u(L,t) = 0$. Then, as superposition preserves the differential equation and the homogeneous side conditions, we will try to build up a solution from these building blocks to solve the nonhomogeneous initial condition $u(x,0) = f(x)$.

First we plug $u(x,t) = X(x)T(t)$ into the heat equation to obtain

$$\begin{equation*} X(x)T'(t) = k X''(x)T(t) . \end{equation*}$$ We rewrite as $$\begin{equation*} \frac{T'(t)}{k T(t)} = \frac{X''(x)}{X(x)}, \end{equation*}$$ which means that both sides of this equation must equal a constant value $-\lambda$.

We obtain the two equations

$$\begin{equation*} \frac{T'(t)}{k T(t)} = -\lambda = \frac{X''(x)}{X(x)} . \end{equation*}$$ In other words, $$\begin{align*} X''(x) + \lambda X(x) &= 0 , \\ T'(t) + \lambda k T(t) &= 0 . \end{align*}$$

The boundary condition $u(0,t) = 0$ implies $X(0)T(t) = 0$. We are looking for a nontrivial solution and so we can assume that $T(t)$ is not identically zero. Hence $X(0) = 0$. Similarly, $u(L,t) = 0$ implies $X(L) = 0$. We are looking for nontrivial solutions $X$ of the eigenvalue problem $X'' + \lambda X = 0$, $X(0) = 0$, $X(L) = 0$. We have previously found that the only eigenvalues are $\lambda _n = \frac {n^2 \pi ^2}{L^2}$, for integers $n \geq 1$, where eigenfunctions are $\sin \left (\frac {n \pi }{L} x\right )$. Hence, let us pick the solutions

$$\begin{equation*} X_n (x) = \sin \left(\frac{n \pi}{L} x \right) . \end{equation*}$$ The corresponding $T_n$ must satisfy the equation $$\begin{equation*} T_n'(t) + \frac{n^2 \pi^2}{L^2} k T_n(t) = 0 . \end{equation*}$$ This is one of our subsection:fourfundamental, and the solution is just an exponential: $$\begin{equation*} T_n(t) = e^{\frac{-n^2 \pi^2}{L^2} k t} . \end{equation*}$$ It will be useful to note that $T_n(0) = 1$. Our building-block solutions are $$\begin{equation*} u_n(x,t) = X_n(x)T_n(t) = \sin \left( \frac{n \pi}{L} x \right) e^{\frac{-n^2 \pi^2}{L^2} k t} . \end{equation*}$$

We note that $u_n(x,0) = \sin \left ( \frac {n \pi }{L} x \right )$. Let us write $f(x)$ as the sine series

$$\begin{equation*} f(x) = \sum_{n=1}^\infty b_n \sin \left(\frac{n \pi}{L} x \right) . \end{equation*}$$ That is, we find the Fourier series of the odd periodic extension of $f(x)$. We used the sine series as it corresponds to the eigenvalue problem for $X(x)$ above. Finally, we use superposition to write the solution as $$\begin{equation*} %\mybxbg{~~ u(x,t) = \sum_{n=1}^\infty b_n u_n(x,t) = \sum_{n=1}^\infty b_n \sin \left( \frac{n \pi}{L} x \right) e^{\frac{-n^2 \pi^2}{L^2} k t} . %~~} \end{equation*}$$

__

Why does this solution work? First note that it is a solution to the heat equation by superposition. It satisfies $u(0,t) = 0$ and $u(L,t) = 0$, because $x=0$ or $x=L$ makes all the sines vanish. Finally, plugging in $t=0$, we notice that $T_n(0) = 1$ and so

$$\begin{equation*} u(x,0) = \sum_{n=1}^\infty b_n u_n(x,0) = \sum_{n=1}^\infty b_n \sin \left( \frac{n \pi}{L} x \right) = f(x) . \end{equation*}$$

Solution: We are solving the following PDE problem:

$$\begin{align*} & u_t = 0.003 u_{xx} , \\ & u(0,t) = u(1,t) = 0 , \\ & u(x,0) = 50\,x\,(1-x) \qquad \text{for } \; 0 < x < 1 . \end{align*}$$

We write $f(x) = 50\,x\,(1-x)$ for $0 < x < 1$ as a sine series. That is, $f(x) = \sum _{n=1}^\infty b_n \sin (n \pi x) ,$ where

$$\begin{equation*} b_n = 2 \int_0^1 50\,x\,(1-x) \sin (n \pi x)dx = \frac{200}{{\pi }^{3}{n}^{3}}-\frac{200\,{\left( -1\right) }^{n}}{{\pi }^{3}{n}^{3}} = \begin{cases} 0 & \text{if } n \text{ even} , \\ \frac{400}{\pi^3 n^3} & \text{if } n \text{ odd} . \end{cases} \end{equation*}$$

The solution $u(x,t)$, plotted in Normally a reference to a previous figure goes here. for $0 \leq t \leq 100$, is given by the series:

$$\begin{equation*} u(x,t) = \sum_{\substack{n=1 \\ n \text{ odd}}}^\infty \frac{400}{\pi^3 n^3} \sin (n \pi x ) e^{-n^2 \pi^2 0.003 t} . \end{equation*}$$

---

---

Finally, let us answer the question about the maximum temperature. It is relatively easy to see that the maximum temperature at any fixed time is always at $x=0.5$, in the middle of the wire. The plot of $u(x,t)$ confirms this intuition. If we plug in $x=0.5$, we get

$$\begin{equation*} u(0.5,t) = \sum_{\substack{n=1 \\ n \text{ odd}}}^\infty \frac{400}{\pi^3 n^3} \sin (n \pi0.5 ) e^{-n^2 \pi^2 0.003 t} . \end{equation*}$$ For $n=3$ and higher (remember $n$ is only odd), the terms of the series are insignificant compared to the first term. The first term in the series is already a very good approximation of the function. Hence $$\begin{equation*} u(0.5,t) \approx \frac{400}{\pi^3} e^{-\pi^2 0.003 t} . \end{equation*}$$ The approximation gets better and better as $t$ gets larger as the other terms decay much faster. Let us plot the function $u(0.5,t)$, the temperature at the midpoint of the wire at time $t$, in Normally a reference to a previous figure goes here.. The figure also plots the approximation by the first term.

---

---

After $t=5$ or so it would be hard to tell the difference between the first term of the series for $u(x,t)$ and the real solution $u(x,t)$. This behavior is a general feature of solving the heat equation. If you are interested in behavior for large enough $t$, only the first one or two terms may be necessary.

Let us get back to the question of when is the maximum temperature one half of the initial maximum temperature. That is, when is the temperature at the midpoint $\frac {12.5}{2} = 6.25$. We notice on the graph that if we use the approximation by the first term we will be close enough. We solve

$$\begin{equation*} 6.25 = \frac{400}{\pi^3} e^{-\pi^2 0.003 t} . \end{equation*}$$ That is, $$\begin{equation*} t = \frac{\ln \frac{6.25\,\pi^3}{400}}{-\pi^2 0.003} \approx 24.5 . \end{equation*}$$ So the maximum temperature drops to half at about $t=24.5$.

___

We mention an interesting behavior of the solution to the heat equation. The heat equation “smoothes” out the function $f(x)$ as $t$ grows. For a fixed $t$, the solution is a Fourier series with coefficients $b_n e^{\frac {-n^2 \pi ^2}{L^2} k t}$. If $t > 0$, then these coefficients go to zero faster than any $\frac {1}{n^p}$ for any power $p$. In other words, the Fourier series has infinitely many derivatives everywhere. Thus even if the function $f(x)$ has jumps and corners, then for a fixed $t > 0$, the solution $u(x,t)$ as a function of $x$ is as smooth as we want it to be.

Insulated ends

Now suppose the ends of the wire are insulated. In this case, we are solving the equation

$$\begin{equation*} u_t = k u_{xx} \qquad \text{with} \quad u_x(0,t) = 0, \quad u_x(L,t) = 0, \quad \text{and} \quad u(x,0) = f(x) . \end{equation*}$$ Yet again we try a solution of the form $u(x,t) = X(x)T(t)$. By the same procedure as before we plug into the heat equation and arrive at the following two equations $$\begin{align*} X''(x) + \lambda X(x) &= 0 , \\ T'(t) + \lambda k T(t) &= 0 . \end{align*}$$

At this point the story changes slightly. The boundary condition $u_x(0,t) = 0$ implies $X'(0)T(t) = 0$. Hence $X'(0) = 0$. Similarly, $u_x(L,t) = 0$ implies $X'(L) = 0$. We are looking for nontrivial solutions $X$ of the eigenvalue problem $X'' + \lambda X = 0$, $X'(0) = 0$, $X'(L) = 0$. We have previously found that the only eigenvalues are $\lambda _n = \frac {n^2 \pi ^2}{L^2}$, for integers $n \geq 0$, where eigenfunctions are $\cos \left ( \frac {n \pi }{L} x\right )$ we include the constant eigenfunction). Hence, let us pick solutions

$$\begin{equation*} X_n (x) = \cos \left( \frac{n \pi}{L} x \right) \qquad \text{and} \qquad X_0 (x) = 1. \end{equation*}$$ The corresponding $T_n$ must satisfy the equation $$\begin{equation*} T_n'(t) + \frac{n^2 \pi^2}{L^2} k T_n(t) = 0 . \end{equation*}$$ For $n \geq 1$, as before, $$\begin{equation*} T_n(t) = e^{\frac{-n^2 \pi^2}{L^2} k t} . \end{equation*}$$ For $n = 0$, we have $T_0'(t) = 0$ and hence $T_0(t) = 1$. Our building-block solutions are $$\begin{equation*} u_n(x,t) = X_n(x)T_n(t) = \cos \left( \frac{n \pi}{L} x \right) e^{\frac{-n^2 \pi^2}{L^2} k t} , \end{equation*}$$ and $$\begin{equation*} u_0(x,t) = 1 . \end{equation*}$$

We note that $u_n(x,0) = \cos \left ( \frac {n \pi }{L} x \right )$. Let us write $f$ using the cosine series

$$\begin{equation*} f(x) = \frac{a_0}{2} + \sum_{n=1}^\infty a_n \cos \left( \frac{n \pi}{L} x \right) . \end{equation*}$$ That is, we find the Fourier series of the even periodic extension of $f(x)$.

We use superposition to write the solution as

$$\begin{equation*} %\mybxbg{~~ u(x,t) = \frac{a_0}{2} + \sum_{n=1}^\infty a_n u_n(x,t) = \frac{a_0}{2} + \sum_{n=1}^\infty a_n \cos \left( \frac{n \pi}{L} x \right) e^{\frac{-n^2 \pi^2}{L^2} k t} . %~~} \end{equation*}$$

Solution: For this problem, we must find the cosine series of $u(x,0)$. For $0 < x < 1$ we have

$$\begin{equation*} 50x\,(1-x) = \frac{25}{3} + \sum_{\substack{n=2 \\ n \text{ even}}}^\infty \left( \frac{-200}{\pi^2 n^2} \right) \cos (n \pi x) . \end{equation*}$$ The calculation is left to the reader. Hence, the solution to the PDE problem, plotted in Normally a reference to a previous figure goes here., is given by the series $$\begin{equation*} u(x,t) = \frac{25}{3} + \sum_{\substack{n=2 \\ n \text{ even}}}^\infty \left( \frac{-200}{\pi^2 n^2} \right) \cos ( n \pi x) e^{-n^2 \pi^2 0.003 t} . \end{equation*}$$

---

---

Note in the graph that as time goes on, the temperature evens out across the wire. Eventually, all the terms except the constant die out, and you will be left with a uniform temperature of $\frac {25}{3} \approx 8.33$ along the entire length of the wire.

___

Let us expand on the last point. The constant term in the series is

$$\begin{equation*} \frac{a_0}{2} = \frac{1}{L} \int_0^L f(x) dx . \end{equation*}$$ In other words, $\frac {a_0}{2}$ is the average value of $f(x)$, that is, the average of the initial temperature. As the wire is insulated everywhere, no heat can get out, no heat can get in. So the temperature tries to distribute evenly over time, and the average temperature must always be the same, in particular it is always $\frac {a_0}{2}$. As time goes to infinity, the temperature goes to the constant $\frac {a_0}{2}$ everywhere.