% =============================================================================
% ch12_heat_equation.tex
% Chapter 12: The Heat Equation
% =============================================================================

\section{Heat Equation}
\label{ch:heat_equation}

\subsection{Physical Derivation}
\label{sec:ch12_physical_derivation}

The heat equation is the prototypical parabolic partial differential equation (PDE). It describes how temperature diffuses through a material over time, and it serves as a mathematical model for many other diffusion processes --- from the spread of a pollutant in a river to the flow of electrical charge through a semiconductor.

\paragraph{Fourier's law of heat conduction.} In 1822, Joseph Fourier established the fundamental law governing heat conduction. Consider a thin rod aligned along the $x$-axis. Let $u(x,t)$ denote the temperature at position $x$ and time $t$. Fourier observed that heat flows from hot regions to cold regions, and that the \textbf{heat flux} $J$ (amount of heat energy flowing per unit area per unit time) is proportional to the temperature gradient:
\begin{equation}
  \label{eq:fouriers_law_1d}
  J(x,t) = -\kappa\,\pd{u}{x}(x,t).
\end{equation}
Here $\kappa > 0$ is the \textbf{thermal conductivity} of the material (units: $\mathrm{W/(m{\cdot}K)}$). The minus sign is essential: heat flows in the direction of decreasing temperature, i.e., opposite to the temperature gradient.

\paragraph{Energy conservation in a rod element.} Now consider a small segment of the rod from $x$ to $x+\Delta x$. The amount of heat energy $E$ contained in this segment is
\[
  E = \rho\,c\,A\,\Delta x \cdot u(x,t),
\]
where $\rho$ is the mass density of the material, $c$ is the specific heat capacity (energy per unit mass per degree), and $A$ is the cross-sectional area of the rod.

The rate of change of the energy in this segment must equal the net heat flux into the segment (energy conservation):
\[
  \frac{\diff E}{\diff t} = A\bigl[J(x,t) - J(x+\Delta x, t)\bigr].
\]
Substituting the expressions:
\[
  \rho\,c\,A\,\Delta x\,\pd{u}{t}(x,t)
  = A\bigl[J(x,t) - J(x+\Delta x, t)\bigr].
\]
Divide by $A\,\Delta x$:
\[
  \rho\,c\,\pd{u}{t}(x,t) = -\frac{J(x+\Delta x, t) - J(x,t)}{\Delta x}.
\]
Taking the limit $\Delta x \to 0$, the right side becomes $-\pd{J}{x}$:
\[
  \rho\,c\,\pd{u}{t} = -\pd{J}{x}.
\]

\paragraph{Derivation of the heat equation.} Substitute Fourier's law \cref{eq:fouriers_law_1d} into the energy conservation equation:
\[
  \rho\,c\,\pd{u}{t} = -\pd{}{x}\bigl(-\kappa\,\pd{u}{x}\bigr)
  = \kappa\,\frac{\partial^2 u}{\partial x^2}.
\]
Assuming the material is homogeneous ($\kappa$, $\rho$, and $c$ are constant), we divide by $\rho\,c$ to obtain the \textbf{one-dimensional heat equation}:
\begin{equation}
  \label{eq:heat_equation_1d}
  \pd{u}{t} = \alpha\,\frac{\partial^2 u}{\partial x^2}, \qquad 0 < x < L, \;\; t > 0,
\end{equation}
where
\begin{equation}
  \label{eq:thermal_diffusivity}
  \alpha = \frac{\kappa}{\rho\,c}
\end{equation}
is the \textbf{thermal diffusivity} (units: $\mathrm{m^2/s}$).

\begin{keyresult}
  \textbf{Physical interpretation of $\alpha$.} The thermal diffusivity $\alpha$ governs the rate at which temperature disturbances propagate through a material. A large $\alpha$ means heat diffuses quickly (the material is a good conductor relative to its heat capacity). A small $\alpha$ means temperature changes propagate slowly. In the units of the heat equation, the characteristic diffusion time across a distance $L$ is $t_{\mathrm{diff}} \sim L^2/\alpha$.
\end{keyresult}

\paragraph{Higher dimensions.} In three dimensions, Fourier's law becomes $\bm{J} = -\kappa\nabla u$ and energy conservation gives
\[
  \pd{u}{t} = \alpha\,\nabla^2 u = \alpha\left(\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} + \frac{\partial^2 u}{\partial z^2}\right).
\]
Throughout this chapter, we focus on the one-dimensional case \cref{eq:heat_equation_1d}, which captures the essential mathematical structure.

\subsection{Separation of Variables}
\label{sec:ch12_separation_of_variables}

We now solve the heat equation \cref{eq:heat_equation_1d} on a finite rod $[0,L]$ with homogeneous boundary conditions. The primary method is \textbf{separation of variables}.

\paragraph{The method.} We seek solutions of the form
\[
  u(x,t) = X(x)\,T(t),
\]
where $X(x)$ depends only on space and $T(t)$ depends only on time. Substitute this ansatz into \cref{eq:heat_equation_1d}:
\[
  X(x)\,T'(t) = \alpha\,X''(x)\,T(t).
\]
Assuming neither factor vanishes identically, divide by $\alpha\,X(x)\,T(t)$:
\begin{equation}
  \label{eq:separation_step}
  \frac{T'(t)}{\alpha\,T(t)} = \frac{X''(x)}{X(x)}.
\end{equation}
The left side depends only on $t$, while the right side depends only on $x$. For this equality to hold for all $x$ and $t$, both sides must equal the same \textbf{separation constant}, which we denote by $-\lambda$.

\begin{theorem}[Separation of Variables for the Heat Equation]
  \label{thm:separation_heat}
  Assuming $u(x,t) = X(x)T(t)$, the heat equation \cref{eq:heat_equation_1d} separates into two ordinary differential equations:
  \begin{align}
    T'(t) + \alpha\lambda\,T(t) &= 0, \label{eq:time_ode} \\[6pt]
    X''(x) + \lambda\,X(x) &= 0, \label{eq:space_ode}
  \end{align}
  where $\lambda$ is the separation constant. The choice $\lambda > 0$ is required by the homogeneous boundary conditions and the physical requirement of decay.
\end{theorem}

\paragraph{Justification for the negative sign.} Why do we write the separation constant as $-\lambda$ rather than $+\lambda$? There are three complementary reasons:

\begin{enumerate}
  \item \textbf{Physical reasoning:} Temperature disturbances should decay over time, not grow. If we used $+\lambda > 0$, the time equation $T' = \alpha\lambda T$ would yield $T(t) = e^{\alpha\lambda t}$, an exponentially growing solution, which contradicts the second law of thermodynamics.

  \item \textbf{Boundary conditions:} With homogeneous Dirichlet conditions $X(0) = 0$ and $X(L) = 0$, the spatial equation $X'' + \lambda X = 0$ admits nontrivial solutions only for $\lambda > 0$ (as established in the eigenvalue analysis of \cref{ch:boundary_value_problems}). If $\lambda \leq 0$, only the trivial solution $X \equiv 0$ satisfies both boundary conditions.

  \item \textbf{Consistency:} Using $-\lambda$ gives the time equation $T' = -\alpha\lambda T$, yielding $T(t) = e^{-\alpha\lambda t}$, which decays for $\lambda > 0$.
\end{enumerate}

\paragraph{Solving the separated equations.} The time equation \cref{eq:time_ode} is a simple first-order linear ODE:
\[
  T'(t) = -\alpha\lambda\,T(t)
  \quad\Longrightarrow\quad
  T(t) = A\,e^{-\alpha\lambda t},
\]
where $A$ is an arbitrary constant.

The space equation \cref{eq:space_ode} is exactly the eigenvalue problem studied in \cref{sec:ch11_eigenvalue_problems}. The specific eigenvalues and eigenfunctions depend on the boundary conditions, as we develop in the next subsections.

\subsection{Dirichlet Boundary Conditions}
\label{sec:ch12_dirichlet}

Consider a rod of length $L$ whose ends are held at zero temperature:
\[
  u(0,t) = 0, \qquad u(L,t) = 0, \qquad t > 0.
\]
These are \textbf{homogeneous Dirichlet boundary conditions}. Together with an initial temperature distribution
\[
  u(x,0) = f(x), \qquad 0 < x < L,
\]
we have the initial-boundary value problem:
\begin{equation}
  \label{eq:heat_dirichlet_problem}
  \begin{cases}
    \pd{u}{t} = \alpha\,\frac{\partial^2 u}{\partial x^2}, & 0 < x < L, \;\; t > 0, \\[6pt]
    u(0,t) = 0, \;\; u(L,t) = 0, & t > 0, \\[6pt]
    u(x,0) = f(x), & 0 < x < L.
  \end{cases}
\end{equation}

\paragraph{Eigenvalue problem for $X(x)$.} The boundary conditions $u(0,t) = 0$ and $u(L,t) = 0$ imply $X(0) = 0$ and $X(L) = 0$. The spatial ODE is
\[
  X''(x) + \lambda X(x) = 0, \qquad X(0) = 0, \;\; X(L) = 0.
\]
From \cref{sec:ch11_eigenvalue_problems}, the eigenvalues and eigenfunctions are:
\[
  \lambda_n = \left(\frac{n\pi}{L}\right)^{\!2},
  \qquad
  X_n(x) = \sin\!\left(\frac{n\pi x}{L}\right),
  \qquad n = 1, 2, 3, \dots
\]

\paragraph{Product solutions.} For each eigenvalue $\lambda_n$, the corresponding time factor is
\[
  T_n(t) = e^{-\alpha\lambda_n t}
         = \exp\!\left[-\alpha\left(\frac{n\pi}{L}\right)^{\!2} t\right].
\]
The product solutions are
\[
  u_n(x,t) = \sin\!\left(\frac{n\pi x}{L}\right)\,
             \exp\!\left[-\alpha\left(\frac{n\pi}{L}\right)^{\!2} t\right],
  \qquad n = 1, 2, 3, \dots
\]
Each $u_n(x,t)$ satisfies the PDE and the homogeneous boundary conditions.

\paragraph{General solution.} By linearity, any linear combination of product solutions is also a solution. We form the \textbf{infinite series} (Fourier sine series in space):
\begin{equation}
  \label{eq:heat_dirichlet_solution}
  u(x,t) = \sum_{n=1}^{\infty} b_n\,
           \sin\!\left(\frac{n\pi x}{L}\right)\,
           \exp\!\left[-\alpha\left(\frac{n\pi}{L}\right)^{\!2} t\right].
\end{equation}
This series satisfies the PDE and boundary conditions for any choice of coefficients $\{b_n\}$. To determine the coefficients, we apply the initial condition:
\[
  u(x,0) = \sum_{n=1}^{\infty} b_n\,\sin\!\left(\frac{n\pi x}{L}\right)
         = f(x).
\]
This is a Fourier sine series for $f(x)$ on $[0,L]$. Using the orthogonality of the sine functions (\cref{sec:ch10_orthogonality}), the coefficients are
\[
  b_n = \frac{2}{L}\int_0^L f(x)\,\sin\!\left(\frac{n\pi x}{L}\right)\,\diff x.
\]

\begin{keyresult}
  \label{key:heat_dirichlet}
  \textbf{Heat equation with homogeneous Dirichlet BCs.} For the problem \cref{eq:heat_dirichlet_problem}, the solution is
  \[
    u(x,t) = \sum_{n=1}^{\infty} b_n\,
             \sin\!\left(\frac{n\pi x}{L}\right)\,
             e^{-\alpha(n\pi/L)^2 t},
  \]
  where
  \[
    b_n = \frac{2}{L}\int_0^L f(x)\,\sin\!\left(\frac{n\pi x}{L}\right)\,\diff x.
  \]
\end{keyresult}

\paragraph{Physical interpretation of the series solution.} Each term in the series corresponds to a \textbf{mode} of the temperature distribution. The $n=1$ mode (the fundamental mode) has the lowest decay rate and dominates the long-time behavior:
\[
  u(x,t) \sim b_1\,\sin\!\left(\frac{\pi x}{L}\right)\,
              e^{-\alpha(\pi/L)^2 t}
  \quad\text{as } t \to \infty.
\]
Higher modes ($n \geq 2$) decay much faster because their decay rates scale as $n^2$. After sufficient time, the temperature profile approaches the shape of the fundamental mode.

\paragraph{Worked examples.}

\begin{workedexample}
  Solve the heat equation on a rod of length $L = \pi$ with zero-temperature ends and initial temperature $f(x) = \sin(2x)$:
  \[
    \begin{cases}
      \pd{u}{t} = \alpha\,\frac{\partial^2 u}{\partial x^2}, & 0 < x < \pi, \;\; t > 0, \\[6pt]
      u(0,t) = 0, \;\; u(\pi,t) = 0, & t > 0, \\[6pt]
      u(x,0) = \sin(2x), & 0 < x < \pi.
    \end{cases}
  \]

  \textbf{Solution.} Here $L = \pi$, so the eigenvalues are $\lambda_n = n^2$ and the eigenfunctions are $X_n(x) = \sin(nx)$. The general solution is
  \[
    u(x,t) = \sum_{n=1}^{\infty} b_n\,\sin(nx)\,e^{-\alpha n^2 t}.
  \]
  Apply the initial condition:
  \[
    u(x,0) = \sum_{n=1}^{\infty} b_n\,\sin(nx) = \sin(2x).
  \]
  By orthogonality, $b_n = 0$ for $n \neq 2$ and $b_2 = 1$. Alternatively, compute explicitly:
  \[
    b_n = \frac{2}{\pi}\int_0^{\pi} \sin(2x)\,\sin(nx)\,\diff x.
  \]
  For $n \neq 2$, the integral vanishes by the orthogonality of sines (\cref{thm:orthogonality}). For $n = 2$:
  \[
    b_2 = \frac{2}{\pi}\int_0^{\pi} \sin^2(2x)\,\diff x
        = \frac{2}{\pi}\cdot\frac{\pi}{2} = 1.
  \]
  The solution is
  \[
    u(x,t) = e^{-4\alpha t}\,\sin(2x).
  \]
  This is a single-mode solution: the temperature profile retains its shape (a half-wave of a sine) and simply decays exponentially in amplitude. The decay rate $4\alpha$ corresponds to the second mode.
\end{workedexample}

\begin{workedexample}
  Solve the heat equation on a rod of length $L = \pi$ with zero-temperature ends and a triangular initial temperature distribution:
  \[
    \begin{cases}
      \pd{u}{t} = \alpha\,\frac{\partial^2 u}{\partial x^2}, & 0 < x < \pi, \;\; t > 0, \\[6pt]
      u(0,t) = 0, \;\; u(\pi,t) = 0, & t > 0, \\[6pt]
      u(x,0) = x(\pi - x), & 0 < x < \pi.
    \end{cases}
  \]

  \textbf{Solution.} With $L = \pi$, the solution is
  \[
    u(x,t) = \sum_{n=1}^{\infty} b_n\,\sin(nx)\,e^{-\alpha n^2 t},
  \]
  where
  \[
    b_n = \frac{2}{\pi}\int_0^{\pi} x(\pi - x)\,\sin(nx)\,\diff x.
  \]
  Expand the integrand:
  \[
    b_n = \frac{2}{\pi}\left[\pi\int_0^{\pi} x\sin(nx)\,\diff x
          - \int_0^{\pi} x^2\sin(nx)\,\diff x\right].
  \]
  \textit{First integral.} Integrate by parts with $u = x$, $\diff v = \sin(nx)\,\diff x$:
  \[
    \int_0^{\pi} x\sin(nx)\,\diff x
    = \Bigl[-\frac{x}{n}\cos(nx)\Bigr]_0^{\pi}
      + \frac{1}{n}\int_0^{\pi} \cos(nx)\,\diff x
    = -\frac{\pi}{n}(-1)^n + 0
    = \frac{\pi}{n}(-1)^{n+1}.
  \]
  So the first contribution is
  \[
    \pi\cdot\frac{\pi}{n}(-1)^{n+1} = \frac{\pi^2}{n}(-1)^{n+1}.
  \]

  \textit{Second integral.} Integrate by parts with $u = x^2$, $\diff v = \sin(nx)\,\diff x$:
  \[
    \int_0^{\pi} x^2\sin(nx)\,\diff x
    = \Bigl[-\frac{x^2}{n}\cos(nx)\Bigr]_0^{\pi}
      + \frac{2}{n}\int_0^{\pi} x\cos(nx)\,\diff x.
  \]
  The boundary term gives $-\dfrac{\pi^2}{n}(-1)^n = \dfrac{\pi^2}{n}(-1)^{n+1}$.
  For the remaining integral, use integration by parts again with $u = x$, $\diff v = \cos(nx)\,\diff x$:
  \[
    \int_0^{\pi} x\cos(nx)\,\diff x
    = \Bigl[\frac{x}{n}\sin(nx)\Bigr]_0^{\pi}
      - \frac{1}{n}\int_0^{\pi} \sin(nx)\,\diff x
    = 0 - \frac{1}{n}\Bigl[-\frac{1}{n}\cos(nx)\Bigr]_0^{\pi}
    = \frac{1}{n^2}\bigl((-1)^n - 1\bigr).
  \]
  So the second integral is
  \[
    \int_0^{\pi} x^2\sin(nx)\,\diff x
    = \frac{\pi^2}{n}(-1)^{n+1} + \frac{2}{n}\cdot\frac{(-1)^n - 1}{n^2}
    = \frac{\pi^2}{n}(-1)^{n+1} + \frac{2}{n^3}\bigl((-1)^n - 1\bigr).
  \]

  \textit{Combine.}
  \[
    b_n = \frac{2}{\pi}\left[
           \frac{\pi^2}{n}(-1)^{n+1}
           - \frac{\pi^2}{n}(-1)^{n+1}
           - \frac{2}{n^3}\bigl((-1)^n - 1\bigr)
          \right]
       = \frac{2}{\pi}\cdot\frac{-2}{n^3}\bigl((-1)^n - 1\bigr)
       = \frac{4}{\pi n^3}\bigl(1 - (-1)^n\bigr).
  \]
  Since $1 - (-1)^n = \begin{cases} 2 & n \text{ odd} \\ 0 & n \text{ even} \end{cases}$:
  \[
    b_n = \begin{cases}
      \dfrac{8}{\pi n^3}, & n \text{ odd}, \\[8pt]
      0, & n \text{ even}.
    \end{cases}
  \]
  The solution is
  \[
    u(x,t) = \frac{8}{\pi}\sum_{k=0}^{\infty}
             \frac{\sin\bigl((2k+1)x\bigr)}{(2k+1)^3}
             \exp\!\left[-\alpha(2k+1)^2 t\right].
  \]
  Writing out the first few terms:
  \[
    u(x,t) = \frac{8}{\pi}\Bigl[
             \sin(x)\,e^{-\alpha t}
             + \frac{\sin(3x)}{27}\,e^{-9\alpha t}
             + \frac{\sin(5x)}{125}\,e^{-25\alpha t}
             + \cdots\Bigr].
  \]
  Notice the rapid decay of higher modes: the $n=3$ term decays $9\times$ faster than the fundamental, and the $n=5$ term decays $25\times$ faster.
\end{workedexample}

\begin{figure}[htbp]
  \centering
  \begin{tikzpicture}[scale=0.9]
    % Axes
    \draw[->] (0,0) -- (6.5,0) node[right] {$x$};
    \draw[->] (0,-0.3) -- (0,3.5) node[above] {$u(x,t)$};

    % Rod endpoints
    \draw[dotted, gray] (0,0) -- (0,3.5);
    \draw[dotted, gray] (6,0) -- (6,3.5);
    \node[font=\small, gray] at (0,-0.25) {$0$};
    \node[font=\small, gray] at (6,-0.25) {$L$};

    % Time t=0: triangular profile
    \draw[thick, red] (0,0) -- (3,3.2) -- (6,0);
    \node[font=\small, red, anchor=west] at (4.2,3.3) {$t=0$};

    % Time t=t1: smoothed
    \draw[thick, orange] (0,0)
      .. controls (1,2.5) and (2,2.9) .. (3,2.8)
      .. controls (4,2.7) and (5,1.8) .. (6,0);
    \node[font=\small, orange, anchor=west] at (4.2,2.9) {$t=t_1$};

    % Time t=t2: further decay
    \draw[thick, blue!70] (0,0)
      .. controls (1.5,1.5) and (2.5,1.7) .. (3,1.7)
      .. controls (3.5,1.6) and (4.5,1.0) .. (6,0);
    \node[font=\small, blue!70, anchor=west] at (4.2,1.8) {$t=t_2$};

    % Time t=large: nearly zero
    \draw[thick, gray!60] (0,0)
      .. controls (1.5,0.5) and (2.5,0.6) .. (3,0.6)
      .. controls (3.5,0.55) and (4.5,0.3) .. (6,0);
    \node[font=\small, gray!60, anchor=west] at (4.2,0.7) {$t \to \infty$};

    % Arrow showing time direction
    \draw[->, >=stealth, thick, teal] (6.3,2.5) -- (6.3,0.5)
      node[right, font=\footnotesize, teal] {decay};
  \end{tikzpicture}
  \caption{Temperature profile evolution for the Dirichlet heat equation. The initial triangular distribution (red) smooths out over time as higher modes decay faster. Eventually all temperature dissipates through the zero-temperature boundaries.}
  \label{fig:temperature_evolution}
\end{figure}

\subsection{Neumann Boundary Conditions}
\label{sec:ch12_neumann}

Now consider a rod whose ends are \textbf{insulated}, meaning no heat can flow through the endpoints. The boundary conditions are
\[
  \pd{u}{x}(0,t) = 0, \qquad \pd{u}{x}(L,t) = 0, \qquad t > 0.
\]
These are \textbf{homogeneous Neumann boundary conditions}. Physically, $\pd{u}{x} = 0$ at an endpoint means the temperature gradient vanishes there, so there is no heat flux ($J = -\kappa\,\pd{u}{x} = 0$).

\paragraph{Eigenvalue analysis.} The spatial ODE with Neumann conditions is
\[
  X''(x) + \lambda X(x) = 0, \qquad X'(0) = 0, \;\; X'(L) = 0.
\]
As derived in \cref{sec:ch11_eigenvalue_problems}, the eigenvalues and eigenfunctions are:
\[
  \lambda_0 = 0, \;\; X_0(x) = 1,
  \qquad\text{and}\qquad
  \lambda_n = \left(\frac{n\pi}{L}\right)^{\!2}, \;\;
  X_n(x) = \cos\!\left(\frac{n\pi x}{L}\right),
  \qquad n = 1, 2, 3, \dots
\]

\paragraph{The $\lambda = 0$ mode.} The eigenvalue $\lambda_0 = 0$ deserves special attention. The time equation for this mode is
\[
  T_0'(t) + \alpha\cdot 0\cdot T_0(t) = 0
  \quad\Longrightarrow\quad
  T_0'(t) = 0
  \quad\Longrightarrow\quad
  T_0(t) = \text{constant}.
\]
This means the $n=0$ mode is a \textbf{steady-state component} that does not decay. Physically, it represents the average temperature of the rod, which is conserved because no heat can escape through the insulated ends.

\paragraph{General solution.}
\begin{equation}
  \label{eq:heat_neumann_solution}
  u(x,t) = b_0 + \sum_{n=1}^{\infty} b_n\,
           \cos\!\left(\frac{n\pi x}{L}\right)\,
           \exp\!\left[-\alpha\left(\frac{n\pi}{L}\right)^{\!2} t\right].
\end{equation}

\begin{keyresult}
  \textbf{Heat equation with homogeneous Neumann BCs.} For the problem
  \[
    \pd{u}{t} = \alpha\,\frac{\partial^2 u}{\partial x^2}, \qquad
    u_x(0,t) = 0,\; u_x(L,t) = 0, \qquad
    u(x,0) = f(x),
  \]
  the solution is \cref{eq:heat_neumann_solution} with
  \[
    b_0 = \frac{1}{L}\int_0^L f(x)\,\diff x,
    \qquad
    b_n = \frac{2}{L}\int_0^L f(x)\,\cos\!\left(\frac{n\pi x}{L}\right)\,\diff x,
    \quad n \geq 1.
  \]
  As $t \to \infty$, the exponentially decaying terms vanish and
  \[
    \lim_{t\to\infty} u(x,t) = b_0 = \frac{1}{L}\int_0^L f(x)\,\diff x,
  \]
  i.e., the temperature approaches the \textbf{average} of the initial distribution.
\end{keyresult}

\paragraph{Physical interpretation.} With insulated ends, the total heat energy in the rod is conserved. The temperature distribution smooths out as higher-frequency modes decay, but the overall average temperature remains fixed. In the long run, the rod reaches a uniform temperature equal to the initial average.

\paragraph{Worked example.}

\begin{workedexample}
  Solve the heat equation on a rod of length $L = \pi$ with insulated ends and initial temperature $f(x) = \cos(2x) + 3$:
  \[
    \begin{cases}
      \pd{u}{t} = \alpha\,\frac{\partial^2 u}{\partial x^2}, & 0 < x < \pi, \;\; t > 0, \\[6pt]
      \pd{u}{x}(0,t) = 0, \;\; \pd{u}{x}(\pi,t) = 0, & t > 0, \\[6pt]
      u(x,0) = \cos(2x) + 3, & 0 < x < \pi.
    \end{cases}
  \]

  \textbf{Solution.} With $L = \pi$, the solution has the form
  \[
    u(x,t) = b_0 + \sum_{n=1}^{\infty} b_n\,\cos(nx)\,e^{-\alpha n^2 t}.
  \]
  Compute the coefficients. For $n = 0$:
  \[
    b_0 = \frac{1}{\pi}\int_0^{\pi} \bigl(\cos(2x) + 3\bigr)\,\diff x
        = \frac{1}{\pi}\left[\frac{\sin(2x)}{2} + 3x\right]_0^{\pi}
        = \frac{1}{\pi}\bigl(0 + 3\pi\bigr) = 3.
  \]
  For $n \geq 1$:
  \[
    b_n = \frac{2}{\pi}\int_0^{\pi} \bigl(\cos(2x) + 3\bigr)\,\cos(nx)\,\diff x.
  \]
  Split the integral:
  \[
    b_n = \frac{2}{\pi}\int_0^{\pi} \cos(2x)\cos(nx)\,\diff x
        + \frac{2}{\pi}\int_0^{\pi} 3\cos(nx)\,\diff x.
  \]
  The second integral vanishes: $\int_0^{\pi} \cos(nx)\,\diff x = 0$ for $n \geq 1$.
  For the first integral, orthogonality of cosines gives zero for $n \neq 2$ and:
  \[
    b_2 = \frac{2}{\pi}\int_0^{\pi} \cos^2(2x)\,\diff x
        = \frac{2}{\pi}\cdot\frac{\pi}{2} = 1.
  \]
  The solution is
  \[
    u(x,t) = 3 + e^{-4\alpha t}\cos(2x).
  \]
  \textbf{Check.} At $t = 0$: $u(x,0) = 3 + \cos(2x)$. $\checkmark$

  As $t \to \infty$: $u(x,t) \to 3$, which equals $b_0 = \dfrac{1}{\pi}\int_0^{\pi} (\cos(2x) + 3)\,\diff x = 3$. The constant background temperature $3$ is preserved, while the spatially varying $\cos(2x)$ component decays.
\end{workedexample}

\subsection{Steady-State Solution}
\label{sec:ch12_steady_state}

A \textbf{steady-state solution} is a solution that does not change in time: $\pd{u}{t} = 0$. In the context of the heat equation, this represents the temperature distribution that the system approaches after an infinite amount of time (if such a limit exists).

\paragraph{Derivation.} Setting $\pd{u}{t} = 0$ in \cref{eq:heat_equation_1d} gives
\[
  0 = \alpha\,\frac{\partial^2 u}{\partial x^2} \quad\Longrightarrow\quad \frac{\partial^2 u}{\partial x^2} = 0.
\]
Integrating twice:
\begin{equation}
  \label{eq:steady_state_general}
  u_{\mathrm{ss}}(x) = A\,x + B,
\end{equation}
where $A$ and $B$ are constants determined by the boundary conditions.

\paragraph{Dirichlet BCs with nonzero temperatures.} Suppose the ends of the rod are held at fixed (possibly nonzero) temperatures:
\[
  u(0,t) = T_1, \qquad u(L,t) = T_2.
\]
The steady-state solution satisfies these boundary conditions:
\[
  u_{\mathrm{ss}}(0) = B = T_1,
  \qquad
  u_{\mathrm{ss}}(L) = A\,L + T_1 = T_2.
\]
Solving for $A$:
\[
  A = \frac{T_2 - T_1}{L}.
\]
Therefore:
\begin{equation}
  \label{eq:steady_state_dirichlet}
  u_{\mathrm{ss}}(x) = T_1 + \frac{T_2 - T_1}{L}\,x.
\end{equation}

\begin{keyresult}
  \textbf{Linear steady-state profile.} For the heat equation with fixed end temperatures $u(0) = T_1$ and $u(L) = T_2$, the steady-state temperature distribution is a \textbf{linear function}:
  \[
    u_{\mathrm{ss}}(x) = T_1 + \frac{T_2 - T_1}{L}\,x.
  \]
  This represents a uniform temperature gradient from one end to the other. Heat flows from the hotter end to the colder end at a constant rate.
\end{keyresult}

\paragraph{Physical interpretation.} The linear steady state reflects a balance between the fixed boundary temperatures and the diffusion process. Once the gradient is established, heat flows at a constant rate through the rod (constant flux $J = -\kappa(T_2-T_1)/L$), and the temperature profile no longer changes.

\paragraph{Relation to Neumann BCs.} For insulated ends ($u_x(0) = 0$, $u_x(L) = 0$), the steady state satisfies $A = 0$, giving $u_{\mathrm{ss}}(x) = B$. The constant $B$ equals the average of the initial temperature (as discussed in \cref{sec:ch12_neumann}).

\subsection{Nonhomogeneous Boundary Conditions}
\label{sec:ch12_nonhomogeneous_bc}

When the boundary conditions are nonhomogeneous --- for example, one end held at a nonzero constant temperature --- we cannot directly apply the separation of variables method, which requires homogeneous BCs. The standard approach is a \textbf{shifting technique}: we decompose the solution into a steady-state part that satisfies the nonhomogeneous boundary conditions and a transient part that satisfies homogeneous boundary conditions.

\paragraph{The shifting technique.} Consider the problem
\[
  \begin{cases}
    \pd{u}{t} = \alpha\,\frac{\partial^2 u}{\partial x^2}, & 0 < x < L, \;\; t > 0, \\[6pt]
    u(0,t) = T_1, \;\; u(L,t) = T_2, & t > 0, \\[6pt]
    u(x,0) = f(x), & 0 < x < L.
  \end{cases}
\]
We seek a solution of the form
\[
  u(x,t) = v(x,t) + \phi(x),
\]
where $\phi(x)$ is a time-independent function chosen to satisfy the boundary conditions:
\[
  \phi(0) = T_1, \qquad \phi(L) = T_2.
\]
The simplest choice is the linear steady-state profile \cref{eq:steady_state_dirichlet}:
\[
  \phi(x) = T_1 + \frac{T_2 - T_1}{L}\,x.
\]
Substituting $u = v + \phi$ into the heat equation:
\[
  \pd{v}{t} = \alpha\,\frac{\partial^2 v}{\partial x^2} + \alpha\,\phi''(x).
\]
Since $\phi(x)$ is linear, $\phi''(x) = 0$, so $v$ satisfies the same homogeneous heat equation:
\[
  \pd{v}{t} = \alpha\,\frac{\partial^2 v}{\partial x^2}.
\]
The boundary conditions for $v$ are homogeneous:
\[
  v(0,t) = u(0,t) - \phi(0) = T_1 - T_1 = 0,
  \qquad
  v(L,t) = u(L,t) - \phi(L) = T_2 - T_2 = 0.
\]
The initial condition for $v$ is
\[
  v(x,0) = u(x,0) - \phi(x) = f(x) - \phi(x).
\]

\paragraph{Complete method.}
\begin{enumerate}
  \item Find $\phi(x)$ satisfying the nonhomogeneous BCs (usually the linear steady state).
  \item Set $u(x,t) = v(x,t) + \phi(x)$.
  \item Solve the heat equation for $v(x,t)$ with homogeneous Dirichlet BCs and initial condition $v(x,0) = f(x) - \phi(x)$, using the method of \cref{sec:ch12_dirichlet}.
  \item Recover $u(x,t) = v(x,t) + \phi(x)$.
\end{enumerate}

\paragraph{Worked example.}

\begin{workedexample}
  A metal rod of length $L = \pi$ has one end held at $0^\circ$C and the other at $100^\circ$C. Initially the rod is at a uniform temperature of $50^\circ$C. Find the temperature $u(x,t)$ for $t > 0$.

  \textbf{Solution.} The problem is
  \[
    \begin{cases}
      \pd{u}{t} = \alpha\,\frac{\partial^2 u}{\partial x^2}, & 0 < x < \pi, \;\; t > 0, \\[6pt]
      u(0,t) = 0, \;\; u(\pi,t) = 100, & t > 0, \\[6pt]
      u(x,0) = 50, & 0 < x < \pi.
    \end{cases}
  \]
  \textit{Step 1: Find $\phi(x)$.} The steady-state solution satisfying $u(0) = 0$ and $u(\pi) = 100$ is
  \[
    \phi(x) = 0 + \frac{100 - 0}{\pi}\,x = \frac{100}{\pi}\,x.
  \]

  \textit{Step 2: Define $v(x,t) = u(x,t) - \phi(x)$.} Then $v$ satisfies
  \[
    \begin{cases}
      \pd{v}{t} = \alpha\,\frac{\partial^2 v}{\partial x^2}, & 0 < x < \pi, \\[6pt]
      v(0,t) = 0, \;\; v(\pi,t) = 0, \\[6pt]
      v(x,0) = 50 - \frac{100}{\pi}\,x.
    \end{cases}
  \]

  \textit{Step 3: Solve for $v(x,t)$.} Using the Dirichlet solution from \cref{sec:ch12_dirichlet}:
  \[
    v(x,t) = \sum_{n=1}^{\infty} b_n\,\sin(nx)\,e^{-\alpha n^2 t},
  \]
  where
  \[
    b_n = \frac{2}{\pi}\int_0^{\pi} \left(50 - \frac{100}{\pi}x\right)\sin(nx)\,\diff x.
  \]
  Split the integral:
  \[
    b_n = \frac{2}{\pi}\cdot 50\int_0^{\pi}\sin(nx)\,\diff x
        - \frac{200}{\pi^2}\int_0^{\pi}x\sin(nx)\,\diff x.
  \]
  The first integral:
  \[
    \int_0^{\pi}\sin(nx)\,\diff x = \Bigl[-\frac{\cos(nx)}{n}\Bigr]_0^{\pi}
    = \frac{1 - (-1)^n}{n}.
  \]
  The second integral (from the triangular example above):
  \[
    \int_0^{\pi}x\sin(nx)\,\diff x = \frac{\pi}{n}(-1)^{n+1}.
  \]
  Therefore:
  \[
    b_n = \frac{100}{\pi}\cdot\frac{1 - (-1)^n}{n}
        - \frac{200}{\pi^2}\cdot\frac{\pi}{n}(-1)^{n+1}
    = \frac{100}{\pi n}\Bigl[1 - (-1)^n - 2(-1)^{n+1}\Bigr].
  \]
  Simplify the bracket: $1 - (-1)^n + 2(-1)^n = 1 + (-1)^n$. So:
  \[
    b_n = \frac{100}{\pi n}\bigl(1 + (-1)^n\bigr).
  \]
  This is nonzero only for \textbf{even} $n$: if $n$ is odd, $1 + (-1)^n = 0$; if $n$ is even, $1 + (-1)^n = 2$. Let $n = 2k$:
  \[
    b_{2k} = \frac{200}{\pi(2k)} = \frac{100}{\pi k},
    \qquad k = 1, 2, 3, \dots
  \]

  \textit{Step 4: Recover $u(x,t)$.}
  \[
    u(x,t) = \frac{100}{\pi}\,x
             + \frac{100}{\pi}\sum_{k=1}^{\infty}
               \frac{\sin(2kx)}{k}\,e^{-4\alpha k^2 t}.
  \]
  As $t \to \infty$, the transient part vanishes and
  \[
    u_{\mathrm{ss}}(x) = \frac{100}{\pi}\,x,
  \]
  a linear gradient from $0^\circ$C at $x=0$ to $100^\circ$C at $x=\pi$.

  \textbf{Check at $t=0$:}
  \[
    u(x,0) = \frac{100}{\pi}\,x + \frac{100}{\pi}\sum_{k=1}^{\infty}\frac{\sin(2kx)}{k}.
  \]
  The series $\sum_{k=1}^{\infty}\dfrac{\sin(2kx)}{k}$ is a known Fourier series that equals $\dfrac{\pi - 2x}{2}$ on $(0, \pi)$. Substituting:
  \[
    u(x,0) = \frac{100}{\pi}\,x + \frac{100}{\pi}\cdot\frac{\pi - 2x}{2}
           = \frac{100x}{\pi} + 50 - \frac{100x}{\pi} = 50.
  \]
  $\checkmark$ The initial condition is satisfied.
\end{workedexample}

\subsection{Source Terms}
\label{sec:ch12_source_terms}

The basic heat equation \cref{eq:heat_equation_1d} assumes no internal heat generation. In many practical situations, there is a \textbf{source term} $Q(x,t)$ representing internal heat production (e.g., electrical heating, chemical reactions, nuclear decay):
\begin{equation}
  \label{eq:heat_with_source}
  \pd{u}{t} = \alpha\,\frac{\partial^2 u}{\partial x^2} + Q(x,t), \qquad 0 < x < L, \;\; t > 0.
\end{equation}

\paragraph{Eigenfunction expansion method.} We assume the solution can be expanded in the eigenfunctions of the spatial problem:
\[
  u(x,t) = \sum_{n=1}^{\infty} T_n(t)\,\sin\!\left(\frac{n\pi x}{L}\right),
\]
and similarly expand the source term:
\[
  Q(x,t) = \sum_{n=1}^{\infty} Q_n(t)\,\sin\!\left(\frac{n\pi x}{L}\right),
  \qquad
  Q_n(t) = \frac{2}{L}\int_0^L Q(x,t)\,\sin\!\left(\frac{n\pi x}{L}\right)\,\diff x.
\]
Substituting into \cref{eq:heat_with_source}:
\[
  \sum_{n=1}^{\infty} T_n'(t)\,\sin\!\left(\frac{n\pi x}{L}\right)
  = \alpha\sum_{n=1}^{\infty} T_n(t)\left(-\frac{n^2\pi^2}{L^2}\right)\,
    \sin\!\left(\frac{n\pi x}{L}\right)
  + \sum_{n=1}^{\infty} Q_n(t)\,\sin\!\left(\frac{n\pi x}{L}\right).
\]
By orthogonality, equate coefficients of each eigenfunction:
\begin{equation}
  \label{eq:mode_ode_source}
  T_n'(t) + \alpha\left(\frac{n\pi}{L}\right)^{\!2} T_n(t) = Q_n(t).
\end{equation}
This is a first-order linear ODE for each mode $T_n(t)$, solvable by the integrating factor method.

\begin{keyresult}
  \textbf{Solution of the modal ODE.} The solution to \cref{eq:mode_ode_source} is
  \[
    T_n(t) = e^{-\alpha\lambda_n t}\left[
      T_n(0) + \int_0^t Q_n(\tau)\,e^{\alpha\lambda_n\tau}\,\diff\tau
    \right],
  \]
  where $\lambda_n = (n\pi/L)^2$ and $T_n(0)$ is determined by the initial condition.
\end{keyresult}

\paragraph{Worked example.}

\begin{workedexample}
  Solve the heat equation with a spatially uniform source on a rod of length $L = \pi$ with zero-temperature ends and initially zero temperature:
  \[
    \begin{cases}
      \pd{u}{t} = \alpha\,\frac{\partial^2 u}{\partial x^2} + Q_0, & 0 < x < \pi, \;\; t > 0, \\[6pt]
      u(0,t) = 0, \;\; u(\pi,t) = 0, & t > 0, \\[6pt]
      u(x,0) = 0, & 0 < x < \pi,
    \end{cases}
  \]
  where $Q_0$ is a positive constant representing uniform internal heat generation.

  \textbf{Solution.} With $L = \pi$, the eigenfunctions are $\sin(nx)$ and $\lambda_n = n^2$.

  \textit{Expand the source.} For a constant source $Q(x,t) = Q_0$:
  \[
    Q_n(t) = \frac{2}{\pi}\int_0^{\pi} Q_0\,\sin(nx)\,\diff x
           = \frac{2Q_0}{\pi}\Bigl[-\frac{\cos(nx)}{n}\Bigr]_0^{\pi}
           = \frac{2Q_0}{n\pi}\bigl(1 - (-1)^n\bigr).
  \]
  This is nonzero only for odd $n$:
  \[
    Q_n = \begin{cases}
      \dfrac{4Q_0}{n\pi}, & n \text{ odd}, \\[8pt]
      0, & n \text{ even}.
    \end{cases}
  \]

  \textit{Solve the modal ODE.} With zero initial temperature, $T_n(0) = 0$. For odd $n$:
  \[
    T_n'(t) + \alpha n^2 T_n(t) = \frac{4Q_0}{n\pi}.
  \]
  Using the integrating factor $e^{\alpha n^2 t}$:
  \[
    \frac{\diff}{\diff t}\Bigl[e^{\alpha n^2 t} T_n(t)\Bigr]
    = \frac{4Q_0}{n\pi}\,e^{\alpha n^2 t}.
  \]
  Integrate from $0$ to $t$:
  \[
    e^{\alpha n^2 t} T_n(t) - T_n(0)
    = \frac{4Q_0}{n\pi}\int_0^t e^{\alpha n^2\tau}\,\diff\tau
    = \frac{4Q_0}{n\pi}\cdot\frac{e^{\alpha n^2 t} - 1}{\alpha n^2}.
  \]
  Since $T_n(0) = 0$:
  \[
    T_n(t) = \frac{4Q_0}{n^3\pi\alpha}\Bigl(1 - e^{-\alpha n^2 t}\Bigr).
  \]

  \textit{Reconstruct the solution.} Letting $n = 2k+1$ for odd modes:
  \[
    u(x,t) = \frac{4Q_0}{\pi\alpha}
             \sum_{k=0}^{\infty}
             \frac{1 - e^{-\alpha(2k+1)^2 t}}{(2k+1)^3}\,
             \sin\bigl((2k+1)x\bigr).
  \]

  \textit{Long-time behavior.} As $t \to \infty$, the exponential terms vanish:
  \[
    u_{\mathrm{ss}}(x) = \frac{4Q_0}{\pi\alpha}
                         \sum_{k=0}^{\infty}
                         \frac{\sin\bigl((2k+1)x\bigr)}{(2k+1)^3}.
  \]
  This is the Fourier sine series of a quadratic function. In fact, this series equals $\dfrac{Q_0}{2\alpha}\,x(\pi - x)$, which is the parabolic steady-state profile. This makes physical sense: with uniform heat generation and zero-temperature ends, the temperature reaches a parabolic profile --- hottest in the middle, zero at the ends --- with $\dfrac{\diff^2 u_{\mathrm{ss}}}{\diff x^2} = -\dfrac{Q_0}{\alpha}$, consistent with the steady-state equation $\alpha u'' + Q_0 = 0$.

  The full solution is
  \[
    u(x,t) = \frac{Q_0}{2\alpha}\,x(\pi - x)
             - \frac{4Q_0}{\pi\alpha}\sum_{k=0}^{\infty}
               \frac{e^{-\alpha(2k+1)^2 t}}{(2k+1)^3}
               \sin\bigl((2k+1)x\bigr),
  \]
  showing the approach to the parabolic steady state from an initially cold rod.
\end{workedexample}

\subsection{Summary}
\label{sec:ch12_summary}

The heat equation unifies several concepts developed in \cref{ch:fourier_series} and \cref{ch:boundary_value_problems}: separation of variables, eigenvalue problems, orthogonality, and Fourier series expansion. The choice of boundary conditions determines the eigenfunction basis and the qualitative behavior of the solution.

\begin{table}[htbp]
\centering
\caption{Heat equation solution methods by boundary condition type}
\label{tab:ch12_summary}
\begin{tabular}{l l p{4.5cm}}
\toprule
\textbf{BC type} & \textbf{Eigenfunctions} & \textbf{Key features} \\
\midrule
Dirichlet: &
$\sin(n\pi x/L)$, $n = 1, 2, \dots$ &
Temperature forced to zero at both ends; all modes decay; $u \to 0$ as $t \to \infty$ \\[12pt]
Neumann: &
$1$, $\cos(n\pi x/L)$, $n = 0, 1, 2, \dots$ &
Insulated ends; $\lambda_0 = 0$ mode persists; $u \to$ average of $f(x)$ as $t \to \infty$ \\[12pt]
Mixed (one Dirichlet, one Neumann): &
$\sin\!\bigl((n+\frac{1}{2})\pi x/L\bigr)$, $n = 0, 1, 2, \dots$ &
Eigenvalues $\lambda_n = \bigl((n+\frac{1}{2})\pi/L\bigr)^2$; all modes decay \\[12pt]
Nonhomogeneous Dirichlet: &
Shifting: $u = v + \phi$ &
$\phi(x)$ = steady-state profile; $v$ solved with homogeneous BCs \\[12pt]
Source term $Q(x,t)$: &
Eigenfunction expansion &
Mode-by-mode ODEs: $T_n' + \alpha\lambda_n T_n = Q_n(t)$; solved via integrating factor \\[12pt]
Steady state: &
$u_{\mathrm{ss}}(x) = Ax + B$ &
Linear profile for Dirichlet BCs; constant for Neumann BCs \\
\bottomrule
\end{tabular}
\end{table}

\begin{hintbox}
  \textbf{Problem-solving checklist for the heat equation.}
  \begin{enumerate}
    \item Write down the PDE, boundary conditions, and initial condition.
    \item Classify the boundary conditions (homogeneous/nonhomogeneous, Dirichlet/Neumann).
    \item If BCs are nonhomogeneous, use the shifting technique to homogenize them.
    \item Identify the eigenfunctions and eigenvalues for the spatial problem.
    \item Write the general series solution with undetermined coefficients.
    \item Compute the coefficients from the initial condition using orthogonality.
    \item If a source term is present, expand it in eigenfunctions and solve the modal ODEs.
    \item Interpret the long-time behavior of the solution.
  \end{enumerate}
\end{hintbox}
