Initial Boundary Value Problems for Heat Equations

Let us consider a heat conduction problem for a straight bar of uniform cross section and homogeneous material. With appropriate Cartesian coordinate system, the bar of length ℓ is oriented in the horizontal x-direction from x = 0 to x = ℓ. Suppose further that the sides of the bar are perfectly insulated so that no heat passes through them. We also assume that the cross-sectional dimensions are so small that the temperature u can be considered constant on any given cross section. Then u is a function of the axial coordinate x and the time t.

Graphics3D[Cylinder[{{0,0,0},{3,0,0}}, 1/3]]

The variation of temperature in the bar is governed by the partial differential equation, called the heat equation or diffusion equation:

\begin{equation} \label{EqBheat.1} \frac{\partial u}{\partial t} = \alpha\, \frac{\partial^2 u}{\partial x^2} \qquad \mbox{or for short} \qquad u_{t} = \alpha\, u_{xx} . \end{equation}
In general, a positive coefficient α>0, known as the thermal diffusivity, may depend on spatial variables, temperature, and pressure. However, we assume the thermal diffusivity to be constant in our problems, for simplicity. The parameter α depends only on the material from which the rod is made and is defined by
\[ \alpha = \frac{\kappa}{\rho\cdot c_p} , \]
where Typical values of α are given in the following table.
Material Thermal diffusivity [m2/s] Thermal diffusivity [mm2/s]
Gold 1.27 × 10-4 127
Silver, pure (99,9%) 1.6563 × 10-4 165.63
Sandstone 1.12--1.19 × 10-6 1.15

In order to determine a unique solution to the heat transfer equation, we need to specify its initial temperature at t=0 and boundary conditions at end points. It is assumed that the initial temperature distribution in the bar is known

\begin{equation} \label{EqBheat.2} u(x,0) = f(x) , \end{equation}
where f is some given function. The physics of the problem dictates that this function must also satisfy some regularity conditions to match the temperature at the boundary. However, sometimes the corresponding initial boundary value problem (IBVP for short) can be considered without these regularity conditions with the expense of not uniform convergence of solution.

Initial Boundary Value Problems for heat equations

The general solution of a partial differential equation is very sensitive to the boundary conditions. As such, the same equation may posses different properties under different sets of boundary conditions. In this section, we illustrate this fact by examining one dimensional heat conduction problems with different sets of boundary conditions. First, we consider a long metal bar of length ℓ which is much longer than it is thick. This means that we can treat the distribution of heat as a function of just one spatial variable x ∈ [0, ℓ] and time variable t. The bar is assumed to be insulated along its length so that it can only gain and lose heat through its ends. This means that the temperature distribution u(x, t) depends only on these two variables.

Keep in mind that, throughout this section, we will be solving the same one-dimensional homogeneous partial differential equation, Eq.\eqref{EqBheat.1}, which is called the diffusion equation (also known as the heat transfer equation). The heat conduction equation can now be paired up with a set of boundary conditions, of which we consider three most common types:

 

Diffusion equation with the Dirichlet boundary condition

Consider a heat transfer problem for a thin straight bar (or wire) of uniform cross section and homogeneous material. Let the x-axis be chosen along the axis of the bar, and let x=0 and x=ℓ denote the ends of the bar. Suppose further that the lateral surface of the rod are perfectly insulated so that no heat transfers through them. We also assume that the cross-sectional dimensions are so small that the temperature u can be considered constant on any cross section. Then u is a function only of the axial coordinate x and the time t.

Assuming that no heat is being generated within the rod, the variation of temperature in it is governed by the partial differential equation (of parabolic type):

\[ \frac{\partial u}{\partial t} = \alpha \frac{\partial^2 u}{\partial x^2} \qquad\mbox{or for short} \qquad u_t = \alpha^2 u_{xx} , \qquad 0< x < \ell , \quad 0< t < t^{\ast} < \infty , \]
where α > 0 is a constant known as the thermal diffusivity. In addition, we assume that the ends x=0 and x=ℓ of the rod are held at fixed temperatures. More over, these temperatures are assumed to be zeroes, u(0,t) = 0℃ and u(ℓ,t) = 0℃ during all heat conduction process (for all t). This requirement is essential for the method of separation of variables that will be used to obtain the solution of Eq.\eqref{EqBheat.1}. In the following section, it will be shown how to bypass the homogeneity of the boundary conditions. For now, we consider the following boundary conditions:
\[ u(0,t) =0 , \qquad u(\ell, t) =0 , \quad t> 0. \]
Besides the boundary conditions, we assume that we know the initial temperature distribution in the bar:
\[ u(x,0^{+}) = \lim_{t\downarrow 0} u(x,t) = f(x) , \qquad 0< x < \ell . \]
Therefore, we get the initial boundary value problem (IBVP for short):
\begin{align*} &\mbox{Heat equation:} \qquad &u_t &= \alpha\,u_{xx} \qquad\mbox{for } 0 < x < \ell \mbox{ and } t > 0, \\ &\mbox{Boundary condition:} \qquad &u (0,t) &= 0 , \quad u (\ell , t) = 0, \qquad t > 0 , \\ &\mbox{Initial condition:} \qquad &u(x,0^{+}) &= f(x) , \qquad 0 < x < \ell , \\ &\mbox{Regularity conditions:} \qquad &f(0) &= u (0,0) = 0 , \quad f(\ell ) = u (\ell , 0) = 0, \qquad \mbox{and} \qquad \lim_{t\to +\infty} u(x,t) = 0. \end{align*}
As previously stated, the regularity conditions are dictated by the physics of the problem. They consist of two parts: the first part requires the matching for the initial temperature with the homogeneous boundary conditions to maintain the continuity of the temperature distribution inside the rod. This condition can be disregarded because it does not affect the solvability of the IBVP. However, violation of this condition makes the given problem ill-posed and its numerical approximation may suffer instability or the Gibbs phenomena. The second condition is almost always satisfied because temperature approaches zero with time due to dissipation of energy.

The following animation gives the temperature distribution in a finite bar subject to homogeneous Dirichlet boundary conditions:

The fundamental problem of heat conduction is to find u(x, t) that satisfies the heat equation and subject to the boundary and initial conditions. Under some light conditions on the initial function f, the formulated initial boundary value problem has a unique solution. This problem can be considered as a boundary value problem in xt-plane (see figure below). The solution u(x, t) of the heat equation is sought in the semi-infinite strip 0 < x < ℓ,    t > 0, subject to the requirement that u(x, t) must assume a prescribed values at each point on the boundary of this strip.
Boundary value problem for the heat conduction equation.
     
wave=Plot[3+0.1*Sin[0.5*t*Pi], {t,0,3}, Filling ->Axis, PlotRange-> {0,4}];
left = Graphics[ Text[Style["u(0,t) = 0", FontSize -> 16], {-0.5, 1.5}]];
right = Graphics[ Text[Style["u(\[ScriptL],0) = 0", FontSize -> 16], {3.5, 1.5}]];
ax = {Arrowheads[Large], Arrow[{{-0.3, 0}, {3.6, 0}}]} ;
ax1 = Graphics[{Thick, Black, ax}] ;
ay = {Arrowheads[Large], Arrow[{{0, -0.3}, {0, 4.2}}]} ;
ay1 = Graphics[{Thick, Black, ay}];
txt1 = Graphics[Text[Style["x", FontSize -> 14], {3.4, 0.35}]];
txt2 = Graphics[ Text[Style["t", FontSize -> 14], {0.25, 4.1}]];
txt3 = Graphics[ Text[Style["x = \[ScriptL]", FontSize -> 14], {3, 3.7}]];
delta = Graphics[ Text[Style[ "\!\(\*SubscriptBox[\(u\), \(t\)]\) = \[Alpha] \ \!\(\*SubscriptBox[\(u\), \(xx\)]\) ", FontSize -> 16], {1.5, 1.5}]];
down = Graphics[ Text[Style[ "u(x,0) = f(x)"], {1.5, -0.3}]];
line = Graphics[ {{Black,Thick},Line[{{3,0},{3,3.4}}]}];
Show[ax1, ay1, wave, txt1, txt2, txt3, delta, left, right, down, line, AspectRatio -> 2/3]

As it is seen from the figure, the boundary for the xt-domain is not smooth because it has two corner points (0,0) and (ℓ,0). In this case, we have to impose the regularity conditions to make the corresponding IBVP well-posed.

Note that the equation and boundary conditions are homogeneous---only under these assumptions the separation of variables method is applicable. This gives us a way to seek its solution as superposition of some nontrivial solutions of the auxiliary problem:

\[ u_t = \alpha\, u_{xx} \qquad (0< x < \ell) \quad u(0,t) = u(\ell ,t ) =0. \]
So we seek nontrivial (which means not identically zero) solutions of the above auxiliary problem that are represented as a product of two functions:
\[ u(x,t) = X(x)\, T(t) , \]
where functions X = X(x) depends only on spatial variable x and T = T(t) depends on time variable t. Substituting this product into the heat equation, we get
\[ X(x)\,\dot{T}(t) = \alpha\, X''(x)\,T(t) \qquad \Longrightarrow \qquad \frac{\dot{T}}{\alpha^2 T(t)} = \frac{X'' (x)}{X(x)} , \]
in which the variables are separated. Since the left hand side depends only on time variable t, the right side depends on spatial variable x, and these variables are independent, these functions could be equal only when they are constant. Denoting this constant of separation as −λ, we arrive at two ordinary differential equations:
\[ \begin{split} \dot{T} + \alpha\, \lambda\,T(t) &= 0 , \\ X'' (x) + \lambda\,X(x) &= 0 . \end{split} \]
The boundary condition at x=0 yields
\[ u(0,t) = X(0)\,T(t) =0 \qquad \Longrightarrow \qquad X(0) =0 , \]
otherwise we would get a trivial solution. From another boundary condition at x=ℓ, we have
\[ u(\ell ,t) = X(\ell )\,T(t) =0 \qquad \Longrightarrow \qquad X(\ell ) =0 , \]
Considering together the differential equation for X(x) and the homogeneous boundary conditions, we obtain the familiar Sturm--Liouville problem that we discussed previously:
\[ X'' (x) + \lambda\,X(x) = 0 , \qquad X(0) =0, \quad X(\ell ) =0 . \]
This problem has nontrivial (not identically zero) solutions only for positive values of parameter λ of which we have a discrete number
\[ \lambda_n = \left( \frac{n\pi}{\ell} \right)^2 \qquad \Longrightarrow \qquad X_n = \sin \frac{n\pi x}{\ell} , \qquad n=1,2,\ldots . \]
Therefore, Xn(x) are associated eigenfunctions that could be multiplied by any nonzero constant.

Now turning to differential equation for T and substituting \( n^2 \pi^2 /\ell^2 \) for λ, we have

\[ T' (t) + \alpha \left( \frac{n\pi}{\ell} \right)^2 T(t) =0 . \]
Hence, T(t) = Tn(t) is proportional to the exponent \( \exp \left\{ - \frac{\alpha\, n^2 \pi^2 t}{\ell^2} \right\} , \ n=1,2,\ldots . \) Multiplying functions Tn and Xn, we obtain for each n partial nontrivial solution of the heat equation satisfying the homogeneous boundary conditions of Dirichlet type:
\[ u_n (x,t) = X_n (x)\,T_n (t) = e^{- n^2 \pi^2 \alpha\, t / \ell^2} \,\sin \frac{n\pi x}{\ell} , \qquad n=1,2,\ldots . \]
Each partial nontrivial solution un(x, t) satisfies the heat conduction equation \eqref{EqBheat.1} and the homogeneous boundary conditions. Obviously, an arbitrary linear combination of these functions will also satisfy the homogeneous heat equation and boundary conditions (in our case, they are of Dirichlet type). Now we apply the principle of superposition by seeking a solution to the given IBVP that is a linear combination of all partial nontrivial solutions
\begin{equation} \label{EqBheat.3} u(x,t) = \sum_{n\ge 1} c_n u_n (x,t) = \sum_{n\ge 1} c_n e^{- n^2 \pi^2 \alpha t / \ell^2} \,\sin \frac{n\pi x}{\ell} . \end{equation}
Assuming that the series \eqref{EqBheat.3} converges uniformely (this condition is required to interchange the order of limit and summation), the the sum function u(x, t) satisfies the partial differential equation \( u_t = \alpha\, u_{xx} \) and homogeneous boundary conditions u(0,t) = 0 and u(ℓ,t) = 0 because the individual terms un have the same property. It remains only to satisfy the initial condition
\begin{equation} \label{EqBheat.4} u(x,0^{+}) = \lim_{t\downarrow 0} \sum_{n\ge 1} c_n u_n (x,t) = \sum_{n\ge 1} c_n u_n (x,0^{+}) = \sum_{n\ge 1} c_n \,\sin \frac{n\pi x}{\ell} = f(x) , \qquad 0< x < \ell . \end{equation}
So if we can move the limit under the sign of summation (this is possible when the series converges uniformely), then what is remained to do is to determine the coefficients cn so that the identity
\[ \sum_{n\ge 1} c_n \,\sin \frac{n\pi x}{\ell} = f(x) , \qquad 0< x < \ell , \]
holds In other words, we have to determine coefficients cn so that the Fourier sine series converges to the initial temperature distribution f(x) for 0<x<ℓ. According to Euler--Fourier formula, we have
\begin{equation} \label{EqBheat.5} c_n = \frac{2}{\ell} \int_0^{\ell} f(x)\,\sin \frac{n\pi x}{\ell}\,{\text d}x , \qquad n=1,2,\ldots . \end{equation}
The expression \eqref{EqBheat.3} for the temperature is not easy to analyze, but the negative exponential factor in each term of the series causes the series to converge if the coefficients \eqref{EqBheat.5} do not grow exponentially. If the initial temperature is represented by a function f(x) of bounded variation, then these coefficients decrease and the total series \eqref{EqBheat.3} converges quite rapidly for t ≥ t0 > 0. For small t, the series solution \eqref{EqBheat.3} may converge, but not uniformly, and proving the validity of the initial condition can be challengeable.

This means that the truncated series

\[ v_N (x,t) = \sum_{n=1}^N X_n (x) \,T_n (t) = \sum_{n=1}^N u_n (x,t) \]
provides a good approximation to the true solution in the domain t ≥ t0 > 0 because every term un(x, t) satisfies exactly the heat conduction equation and the boundary conditions, and the the exponential tail of the remainder can be ignored. Therefore, the truncated series can be used for actual approximation to the solution for the moderate values of t, but it loses its fast convergence for small t, which requires additional attention.

Example 1A: Consider a uniform rod of length ℓ = 3 meters long, insulated on sides, with thermal diffusivity 9. If initially the rod has a constant temperature of 20℃, but homogeneous Dirichlet boundary conditions are maintained, we come to the following IBVP:

\begin{equation*} \begin{split} &\text{(PDE)} \quad u_{t} = 9\, u_{xx} , \quad 0<x<3, \ 0< t\le t^{\ast} < \infty , \\ &\text{(BC)} \hspace{6.2mm} u(0,t) =0 , \quad u (3 ,t) =0 , \quad 0< t\le t^{\ast} < \infty , \\ &\text{(IC)} \qquad u(x,0) = 20 . \end{split} \end{equation*}
The temperature in the rod satisfies the heat conduction equation \eqref{EqBheat.1} with ℓ = 3 and initial temperature f(x) = 20 for 0 < x < 3. Thus, from Eq.\eqref{EqBheat.3}, the solution is
\[ u(x,t) = \sum_{n\ge 1} c_n e^{-n^2 \pi^2 t} \sin \frac{n\pi x}{3} , \]
where
\[ c_n = \frac{2}{3} \int_0^3 20\, \sin \frac{n\pi x}{3}\,{\text d} x = \frac{40}{n\pi} \left( 1 - \cos n\pi \right) = \frac{40}{n\pi} \times \begin{cases} 2, & \ \mbox{ for odd indices, } n= 2k+1 , \\ 0, & \ \mbox{ for even indices, } n= 2k. \end{cases} \]
Mathematica helps:
40/3*Integrate[Sin[n*Pi*x/3], {x, 0, 3}]
(40 (3 - 3 Cos[n \[Pi]]))/(3 n \[Pi])
Finally, we obtain
\[ u(x,t) = \frac{80}{\pi} \sum_{k \ge 0} \frac{1}{2k+1} \,e^{-(2k+1)^2 \pi^2 t} \,\sin \frac{\pi \left( 2k+1 \right) x}{3} . \]
The series solution converges rapidly in the domain t ≥ t0 > 0, but the initial conditions are difficult to verify because the series does not converge uniformly for small t and we cannot set t = 0 under the sign of summation because it converges conditionally.

Now we demonstrate how this problem can be solved with Mathematica. Firs, we define the heat partial differential equation that represents the temperature in the bar:

heatEq = D[u[x, t], t] == 9 * D[u[x, t], {x, 2}]
\( u^{(0,1)}[x, t] = 9\,u^{(2,0)}[x, t] \)
Then we solve the heat equation with the appropriate boundary conditions
heatEqSoln = u[x, t] /. DSolve[{heatEq, u[0, t] == 0, u[3, t] == 0, u[x, 0] == 20}, u[x, t], {x, t}][[1, 1]] /. K[1] -> n
which gives exactly the same formula as before.
      We plot the temperature distribution inside the bar.
Plot3D[heatEqSoln30[x, t], {x, 0, 3}, {t, 0, 1} , ColorFunction -> "TemperatureMap" , PlotRange -> All , AxesLabel -> {"Length", "Time", "Temperature"}]
       Temperature distribution.            Mathematica code

As it is seen from the graph, there are ripples near the edge t = 0, which indicates poor convergence of the series. This is a confirmation that the given initial boundary value problem is ill-posed; so the corresponding numerical calculations must be done with care. The Gibbs phenomenon becomes clear when plotting the graph from a different perspective.
      The graph reveals a Gibbs phenomena near the edge t = 0. We take the truncated sum with 31 terms and plot it.
u30[x_, t_] = 80/Pi*Sum[ Exp[-(2*k + 1)^2*Pi^2*t]*Sin[Pi*(2*k + 1)*x/3]/(2*k + 1), {k, 0, 30}];
Plot3D[u30[x, t], {x, 0, 3}, {t, 0, 1} , ColorFunction -> "TemperatureMap" , PlotRange -> All , AxesLabel -> {"Length", "Time", "Temperature"}]
       Temperature distribution.            Mathematica code

In order to eliminate the Gibbs phenomena, we use the Cesàro summation:
\[ C_{30}(x,t) = \frac{80}{\pi} \sum_{k = 0}^{30} \left( 1 - \frac{k}{31} \right) \frac{1}{2k+1} \,e^{-(2k+1)^2 \pi^2 t} \,\sin \frac{\pi \left( 2k+1 \right) x}{3} . \]
      The graph reveals a Gibbs phenomena near the edge t = 0. We take the truncated sum with 31 terms and plot it.
C30[x_, t_] = 80/Pi*Sum[(1 - k/31)* Exp[-(2*k + 1)^2*Pi^2*t]*Sin[Pi*(2*k + 1)*x/3]/(2*k + 1), {k, 0, 30}];
Plot3D[C30[x, t], {x, 0, 3}, {t, 0, 1} , ColorFunction -> "TemperatureMap" , PlotRange -> All , AxesLabel -> {"Length", "Time", "Temperature"}]
       Cesàro summation for temperature approximation.            Mathematica code

We plot temperature as a function of the position across the bar at different time values. First, we compare truncated versions of regular Fourier approximation and its Cesàro counterpart for tiny t = 0.001.
      We put two temperature profiles at t = 0.001 together for comparison,
Plot[{Callout[u30[x, 0.001], "Fourier", LabelStyle -> {Bold, Orange}], Callout[C30[x, 0.001], "Cesaro", LabelStyle -> {Bold, Blue}]}, {x, 0, 3}, PlotStyle -> Thickness[0.01]]
       Fourier and Cesàro approximations.            Mathematica code

Starting with t = 0.01, there is no visible difference in regular Fourier summation and corresponding Cesàro approximation.
Plot[u30[x, 0.1], {x, 0, 3}, PlotStyle -> Thickness[0.01]]
Plot[u30[x, 0.2], {x, 0, 3}, PlotStyle -> Thickness[0.01]]
Plot[u30[x, 2.1], {x, 0, 3}, PlotStyle -> Thickness[0.01]]
Of course, Mathematia complains that the exponential terms are too small for actual calculations. Therefore, 30-term approximation is not needed and first few terms is sufficient for accurate approximations.
           
       Temperature distribution at t = 0.1.            Temperature distribution at t = 0.2.            Temperature distribution at t = 2.1.

In order plot the temperature as a function of the time at different position values along the bar, we compare its distribution with just first term
\[ s_0 (x,t) = \frac{80}{\pi}\, e^{-\pi^2 t} \sin \left( \frac{\pi x}{3} \right) . \]
      We choose a typical point x = 2 and plot the truncated series with 30 terms along with its first term, s0(2, t).
s0[x_, t_] = 80/Pi*Exp[-Pi^2*t]*Sin[Pi*x/3];
Plot[{Callout[u30[2, t], "Fourier", LabelStyle -> {Bold, Orange}], Callout[s0[2, t], "0 term", LabelStyle -> {Bold, Blue}]}, {t, 0, 1}, PlotStyle -> Thickness[0.008], PlotRange -> {0, 20}]
       Temperature distribution with time at x = 2.            Mathematica code

If you looking at 1% accuracy (which is sufficient for plotting), then the figure shows that the single term s0(x, t) gives a good approximation.
      Finally, we demonstrate annimation for temperature distribution with time.
animation1 = Table[Plot[u30[x, t], {x, 0, 3}, PlotRange -> {0, 20}, PlotStyle -> Thickness[0.01] ], {t, 0, 0.5, 0.01}];
Export["1AAnimation.gif", animation1, "AnimationRepetitions" -> \[Infinity]]

    ■