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

κ is thermal conductivity [W/(m·K)]
ρ is density [kg/m³]
c_p is specific heat capacity [j/(kg·K)]

Typical values of α are given in the following table.

Material	Thermal diffusivity [m²/s]	Thermal diffusivity [mm²/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:

The first type of boundary conditions or Dirichlet conditions (named after Peter Gustav Lejeune Dirichlet) when the temperature is prescribed at a boundary
\[ \left. u(x,t) \right\vert_{x\in \mbox{ boundary}} = T(t) . \]
The second type of boundary conditions or Neumann conditions ( named after Carl Neumann) when the heat flow (in general, the surface heat flux) is specified
\[ \left. \nabla u(x,t) \right\vert{x\in \mbox{ boundary}} = \left. \frac{\partial u(x,t)}{\partial x} \right\vert_{x\in \mbox{ boundary}} = q(t) . \]
The third type of boundary conditions or mixed boundary conditions usually include the following two cases.
- On one part of the boundary the temperature is specified (the Dirichlet conditions), while on another part the normal derivative is specified (Neumann conditions).
- When we observe convective heat transfer at the boundary
  \[ \left( - \kappa\, \nabla \,u \right)_{x\in \mbox{ boundary}} = h\,A \left( u - T_{\infty} \right)_{x\in \mbox{ boundary}} , \]
  where κ is the thermal conductivity of the material (with units W/(m K)), h is the convection heat transfer coefficient (with unites Watt/m² -K), and A is the surface area (units m²). The convection heat transfer coefficient is sometimes referred to as a film coefficient. It represents the thermal resistance of a relatively stagnant layer of fluid or air between a heat transfer surface and the fluid medium. It can also be defined through the Nusselt number: \( \displaystyle h = \left( Nu \cdot \kappa \right) /L , \) where Nu is the Nusselt number, κ is the thermal conductivity, and L is the characteristic length.

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^{+}) &= \lim_{t\downarrow 0} u(x,t) = 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, X_n(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) = T_n(t) is proportional to the exponent \( \exp \left\{ - \frac{\alpha\, n^2 \pi^2 t}{\ell^2} \right\} , \ n=1,2,\ldots . \) Multiplying functions T_n and X_n, 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 u_n(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 u_n 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 c_n 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 c_n 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 ≥ t₀ > 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 ≥ t₀ > 0 because every term u_n(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.

x=2


syms x n t
y1 = (x.*(120-x).*sin((n.*pi.*x)./60));
I1 = int(y1,x,0,60);
Cn = ((1/30).*I1);
utemp =
subs(Cn.*sin(n.*pi.*x.*(1/60)).*exp(((n.*pi.*(1/60)).^2).*t.*-4),n,1:1000);
U = sum(utemp);
U = subs(U,t,0:5);
U = subs(U,x,20);
figure(1)
plot(0:5,U)

picture:


syms x n t
y1 = (x.*(120-x).*sin((n.*pi.*x)./60));
I1 = int(y1,x,0,60);
Cn = ((1/30).*I1);
utemp =
subs(Cn.*sin(n.*pi.*x.*(1/60)).*exp(((n.*pi.*(1/60)).^2).*t.*-4),n,1:10);
U = sum(utemp);
U = subs(U,t,0:5);
U = subs(U,x,1);
figure(1)
plot(0:5,U)

picture:


syms x n t
y1 = (x.*(120-x).*sin((n.*pi.*x)./60));
I1 = int(y1,x,0,60);
Cn = ((1/30).*I1);
utemp =
subs(Cn.*sin(n.*pi.*x.*(1/60)).*exp(((n.*pi.*(1/60)).^2).*t.*-4),n,1:10);
U = sum(utemp);
U = subs(U,t,0:5);
U = subs(U,x,1);
figure(1)
plot(0:5,U)
figure(1)
plot(0:5,U)

picture:


syms x n t
y1 = (x.*(120-x).*sin((n.*pi.*x)./60));
I1 = int(y1,x,0,60);
Cn = ((1/30).*I1);
utemp =
subs(Cn.*sin(n.*pi.*x.*(1/60)).*exp(((n.*pi.*(1/60)).^2).*t.*-4),n,1:10);
U = sum(utemp);
U = subs(U,t,1);
U = subs(U,x,0:60);
figure(1)
plot(0:60,U)

picture:


syms x n t
y1 = (x.*(120-x).*sin((n.*pi.*x)./60));
I1 = int(y1,x,0,60);
Cn = ((1/30).*I1);
utemp =
subs(Cn.*sin(n.*pi.*x.*(1/60)).*exp(((n.*pi.*(1/60)).^2).*t.*-4),n,1:10);
U = sum(utemp);
U = subs(U,t,2.5);
U = subs(U,x,0:60);
figure(1)
plot(0:60,U)

picture:


syms x n t
y1 = (x.*(120-x).*sin((n.*pi.*x)./60));
I1 = int(y1,x,0,60);
Cn = ((1/30).*I1);
utemp =
subs(Cn.*sin(n.*pi.*x.*(1/60)).*exp(((n.*pi.*(1/60)).^2).*t.*-4),n,1:10);
U = sum(utemp);
U = subs(U,t,4);
U = subs(U,x,0:60);
figure(1)
plot(0:60,U)

picture: