Solving the Heat Equation

Portland State University

1. Introduction

For a metal rod of unit length, placed (for mathematical convenience) on the real interval $[0,\pi]$ , let $u(x,t)$ denote the temperature of the rod at the cross-section $x$ units from the origin at time $t$ (we assume the temperature to be constant on each cross-section of the rod. Then the function $u$ is known to obey the Heat Equation: $u_t = ku_{xx}$ , where $k$ is a positive constant that is a property of the material used to form the rod.

Suppose, for simplicity, that $k=1$ , and we have these initial-boundary conditions:

(1) the ends of the rod are held at constant temperature zero for all time: $u(0,t)=u(1,t)=0$ for all $t \ge 0$ , and

(2) the initial temperature of the rest of the rod is 1, i.e., $u(x,0) \equiv 1$ on $(0,\pi)$ .

Then the solution to our our initial-boundary problem, consisting of the heat equation $u_t = u_{xx}$ , along with boundary condition (1) and initial condition (2), is given by the series below: $(*)~~~~~ {\displaystyle u(x,t) = \sum_{k=0}^\infty e^{-(2k+1)^2 t}~\frac{\sin((2k+1)x)}{2k+1}}$

To see that this series really does solve our initial-boundary problem for the heat equation, note that formal term-by-term differentiation (justified by the exponential decay of the time-dependent terms) shows that $u$ satisfies the heat equation, while the substitution $x=0$ and $x=\pi$ show that (1) is satisfied. As for (2): the series for $t=0$ is the Fourier sine series for the function $f\equiv 1$ over the interval $[0,\pi]$ . See [1, Sec. 1.8.3-1.8.4, pp. 63--70], [2, Sec. 18.10, especially pp. 999--1004], or [3] for more details, and to learn how to discover this formula for $u(x,t)$ .

2. Plotting the solution

We'll define here a function that returns the sum of the first $n$ terms of the series (*), and another one, that plots this sum as a function of $x$ , using the time $t$ as a parameter.

SAGE commands featured here:

def u(x,t,n): 
    return (4/pi)*sum([exp(-(2*k+1)^2*t)*sin((2*k+1)*x)/(2*k+1) for k in range(n)])

Here is the sum of the first five terms of the series (*)

var('t')
u(x,t,5).show() #Sum of first five terms.

\newcommand{\Bold}[1]{\mathbf{#1}}\frac{4 \, {\left(35 \, e^{\left(-81 \, t\right)} \sin\left(9 \, x\right) + 45 \, e^{\left(-49 \, t\right)} \sin\left(7 \, x\right) + 63 \, e^{\left(-25 \, t\right)} \sin\left(5 \, x\right) + 105 \, e^{\left(-9 \, t\right)} \sin\left(3 \, x\right) + 315 \, e^{\left(-t\right)} \sin\left(x\right)\right)}}{315 \, \pi}

For $t = 0$ the infinite series (*) is the Fourier sine series of the function identically 1 on the interval $(0,\pi)$ . We superimpose the graph of this function (in red) on a plot of the sum of the first 25 terms of that sine series.

fs = u(x,0,25) #Plot of fs approx. to initial function
init_temp =  x^0  #The init. temp distn u(x,0) = 1

Pfs = fs.plot((0,pi+.1))                #Plot fs partial sum in blue
one = init_temp.plot((0,pi),rgbcolor=[1,0,0]) #Plot init temp in red

(Pfs+one).show(xmin=-.1, xmax=pi+.1, ymin=0,figsize=[3,3]) #Show both plots on same set of axes

Now we define a function $P$ that maps the $n$ -th partial sum of (*) to a graphics object having $t$ as parameter (i.e., $P:[0,\infty)\times\{{\rm Natural~numbers}\}\rightarrow \{\rm Graphics~objects\}$ ). The command tests what we've done.

def P(t,n):
    fs = u(x,t,n)
    return fs.plot((-.1,pi+.1),rgbcolor=[sin(t),cos(t),exp(-t)]);  #colors change with t

Next we plot, for t=.01, .5, 1,2,and 3, the 25th partial sum of series for $u(x,t)$ . Compare the plot for $t=.01$ with the one above for $t=0$ , noting how effectively the exponential "multiplier" smooths out the partial sum.

show(P(.01,25)+P(.5,25)+P(1,25)+P(2,25)+P(3,25),ymin=0,ymax=1)

3. Visualizing the solution with

#"continuous" parameter t set automatically for 50 equally spaced values in [0,5]
@interact
def _(t=(0,5)):
     show(P(t,25), ymax=1)

4. Animating the solution via `@interact`

Here's an animation inspired by a Maple program written by Paul Bourdon. We represent the rod by 100 large dots on the $x$ -axis, and the temperature $u(x,t)$ of (the center of) each dot by a color, ranging from red (hot) to white(cold). The distribution of colors at time t represents $u(x,t)$ . The animation shows the temperature at each point cooling from the initial temperature 1 toward zero as $t\rightarrow\infty$ .

The first step is to define the function that colors the 100 consecutive points on the $x$ -axis according to the value $u(x,t)$ .

def cool(t,n):
    g = Graphics()
    for k in range(100):
        v = u(pi*k/100,t,n)
        g += point((pi*k/100,0),pointsize=100, rgbcolor=[1,1-sqrt(v),1-sqrt(v)])
    return g

cool(.05,10).show(axes=false) #Testing our definition

#"continuous" parameter t set automatically for 50 equally spaced values in interval [.02,5]
@interact
def _(t=(0.02,5)):
     show(cool(t,25), axes=false)

References

[1] H. Dym and H.P. McKean, Fourier Series and Integrals, Academic Press 1972.

[2] A. Jeffrey, Advanced Engineering Mathematics, Harcourt/Academic Press 2002.

[3] Weisstein, Eric W. Heat Conduction Equation. From MathWorld--A Wolfram Web Resource. <a href="http://mathworld.wolfram.com/HeatConductionEquation.html" target="_blank"http://mathworld.wolfram.com/HeatConductionEquation.html