Skip to content

Heat Equation

HEAT EQUATION
Consider a region U⊂RnU \subset \mathbb{R}^n that is open and a time variable t>0t > 0, then we can define the heat equation:
∂u∂t=αΔu\begin{equation*}\frac{\partial u}{\partial t} = \alpha \Delta u\end{equation*}
Here, the unknown function u:U‾×[0,∞)→Ru: \overline{U} \times [0,\infty) \to \mathbb{R} represents the temperature distribution, where u(x,t)u(x,t) gives the temperature at position xx and time tt. The Laplacian operator Δ\Delta acts on the spatial variables x=(x1,…,xn)x = (x_1,\ldots,x_n) and is defined as:
Δu=∑i=1n∂2u∂xi2=∂2u∂x12+∂2u∂x22+⋯+∂2u∂xn2\begin{equation*}\Delta u = \sum_{i=1}^n \frac{\partial^2 u}{\partial x_i^2} = \frac{\partial^2 u}{\partial x_1^2} + \frac{\partial^2 u}{\partial x_2^2} + \cdots + \frac{\partial^2 u}{\partial x_n^2}\end{equation*}
The constant α\alpha represents the thermal diffusivity of the medium, which measures how quickly heat spreads through the material.

The foolowing animations are taken from 3blue1brown's video "But what is a partial differential equation? | DE2"

Let's explore the heat equation in the one-dimensional case. Consider a function f:[0,L]→Rf: [0, L] \to \mathbb{R} that represents the initial temperature distribution.
We study a thin metal rod of length LL and analyze how heat flows through it. Initially, the rod is heated non-uniformly, and at time t=0t=0, the temperature at position xx is f(x)f(x). The function u:[0,L]×[0,∞)→Ru: [0,L] \times [0,\infty) \to \mathbb{R} describes the temperature of the rod at position xx and time tt, with the initial condition:
u(x,0)=f(x)for0≤x≤L\begin{equation*}u(x,0) = f(x) \quad \text{for} \quad 0 \leq x \leq L\end{equation*}
The changing of T2T_2 depends on T1+T32\frac{T_1 + T_3}{2}
Physical intuition tells us that over time, the rod should reach a uniform temperature distribution. We want to describe this behavior mathematically. For simplicity, let's consider three adjacent points x1<x2<x3x_1 < x_2 < x_3 with their temperatures at some time tt:
T1=u(t,x1),T2=u(t,x2),T3=u(t,x3).\begin{equation*}T_1 = u(t, x_1), \quad T_2 = u(t, x_2), \quad T_3 = u(t, x_3).\end{equation*}
Let's investigate the temperature change dT2\text{d}T_2 over a small time interval dt\text{d}t relative to the temperatures at points T1T_1 and T3T_3. The fundamental principle is that T2T_2 should approach the average of its neighboring points, T1+T32\frac{T_1 + T_3}{2}, as the system moves toward thermal equilibrium. This can be approximated as:
dT2dt=α(T1+T32−T2)=α2((T3−T2)−(T2−T1))\begin{equation*}\frac{\text{d}T_2}{\text{d}t} = \alpha \left(\frac{T_1 + T_3}{2} - T_2\right) = \frac{\alpha}{2}\left((T_3 - T_2) - (T_2 - T_1)\right)\end{equation*}
The last term represents a finite difference approximation of the Laplacian operator Δ\Delta.
While this explanation is not mathematically rigorous, it provides an intuitive understanding of the heat equation's behavior.
The 3D surface of the solution u(x,t)u(x,t) to the one-dimensional heat equation with initial distribution u(x,0)=f(x)u(x, 0) = f(x). The solution was found using finite difference method.
To visualize the evolution of temperature in our one-dimensional rod, we need to consider the heat equation with the given initial function f(x)f(x). The temperature function u(x,t)u(x,t) requires a three-dimensional coordinate system (x,t,u(x,t))(x, t, u(x,t)) for complete visualization, where:
  • xx represents the position along the rod
  • tt represents time
  • u(x,t)u(x,t) gives the temperature at position xx and time tt

The foolowing material is taken from [EvansPDE2010] p. 46-47

Now, as we understand the idea behind the heat equation we want to find the solution to it. The key to solving the heat equation lies in understanding the fundamental solution Φ(x,t)\Phi(x,t). This remarkable function not only solves the heat equation itself but also enables us to construct solutions for any given initial distribution function f(x)f(x).
FUNDAMENTAL SOLUTION
The function
Φ(x,t):={1(4πt)n/2e−∣x∣24t(x∈Rn,  t>0)0(x∈Rn,  t<0)\begin{equation*}\Phi(x,t) := \begin{cases}\frac{1}{(4\pi t)^{n/2}} e^{-\frac{|x|^2}{4t}} & (x \in \mathbb{R}^n, \; t > 0) \\0 & (x \in \mathbb{R}^n, \; t < 0)\end{cases}\end{equation*}
is called the fundamental solution of the heat equation. Note that Φ\Phi has a singularity at the point (0,0)(0,0).
Using Φ\Phi, we can construct solutions to the initial-value (or Cauchy) problem:
ut=Δuin Rn×(0,∞)u=gon Rn×{t=0}.\begin{align*}u_t &= \Delta u \quad \text{in } \mathbb{R}^n \times (0,\infty) \\u &= g \quad \text{on } \mathbb{R}^n \times \{t=0\}.\end{align*}
A key observation is that for any fixed y∈Rny \in \mathbb{R}^n and t>0t > 0, both (x,t)↦Φ(x,t)(x,t) \mapsto \Phi(x,t) and (x,t)↦Φ(x−y,t)(x,t) \mapsto \Phi(x-y,t) solve the heat equation away from their respective singularities. This property suggests that we can construct a solution through convolution:
u(x,t)=∫RnΦ(x−y,t)g(y) dy=1(4πt)n/2∫Rne−∣x−y∣24tg(y) dy(x∈Rn,  t>0)\begin{equation*}u(x,t) = \int_{\mathbb{R}^n} \Phi(x-y,t)g(y)\,\text{d}y = \frac{1}{(4\pi t)^{n/2}} \int_{\mathbb{R}^n} e^{-\frac{|x-y|^2}{4t}}g(y)\,\text{d}y \quad (x \in \mathbb{R}^n, \; t > 0)\end{equation*}
The fact that this convolution indeed provides a solution to our initial-value problem will be established in the subsequent theorem.

The foolowing material is taken from [EvansPDE2010] p. 47-48

Using the definition of the fundamental solution we can establish the following facts:
1. Regularity of the solution u(x,t)u(x,t), i.e. the amount of derivatives.
2. Next, we verify that the convolution with Φ(x,t)\Phi(x,t) indeed yields a solution to the heat equation.
3. Finally, despite the singularity of Φ(x,t)\Phi(x,t) at (0,0)(0,0), we will prove that our solution satisfies the initial condition by showing
lim⁡t→0+u(x,t)=g(x)\begin{equation*}\lim_{t \to 0^+} u(x,t) = g(x)\end{equation*}
for all x∈Rnx \in \mathbb{R}^n.
SOLUTION OF INITIAL-VALUE PROBLEM
Assume g∈C(Rn)∩L∞(Rn)g \in C(\mathbb{R}^n) \cap L^{\infty}(\mathbb{R}^n)(continuous and bounded), and define uu as above. Then
  • u∈C∞(Rn×(0,∞))u \in C^{\infty}(\mathbb{R}^n \times (0,\infty)),
  • ut(x,t)−Δu(x,t)=0(x∈Rn,t>0)u_t(x,t) - \Delta u(x,t) = 0 \quad (x \in \mathbb{R}^n, t > 0),
  • lim⁡(x,t)→(x0,0)x∈Rn,t>0u(x,t)=g(x0)for each point x0∈Rn\displaystyle\lim_{\substack{(x,t)\to(x^0,0) \\ x\in\mathbb{R}^n, \, t>0}} u(x,t) = g(x^0) \quad \text{for each point } x^0 \in \mathbb{R}^n.
\quad Proof. 1. Since the function 1tn/2e−∣x∣24t\frac{1}{t^{n/2}}e^{-\frac{|x|^2}{4t}} is infinitely differentiable, with uniformly bounded derivatives of all orders, on Rn×[δ,∞)\mathbb{R}^n \times [\delta,\infty) for each δ>0\delta > 0. As a consequence, our solution u(x,t)u(x,t), obtained through convolution with this function, inherits these smoothness properties. Specifically:
u∈C∞(Rn×(0,∞)).\begin{equation*}u \in C^\infty(\mathbb{R}^n \times (0,\infty)).\end{equation*}
This means our solution possesses continuous derivatives of all orders in both space and time variables, making it an extremely well-behaved function from an analytical perspective.
2. Let us now verify that our solution satisfies the heat equation. For each fixed y∈Rny \in \mathbb{R}^n, we know that Φ(x−y,t)\Phi(x-y,t) solves the heat equation, meaning:
Φt(x−y,t)=ΔxΦ(x−y,t)\begin{equation*}\Phi_t(x-y,t) = \Delta_x\Phi(x-y,t)\end{equation*}
Using this property and the linearity of integration, we can show that u(x,t)u(x,t) solves the heat equation:
ut(x,t)−Δu(x,t)=∫Rn[Φt−ΔxΦ](x−y,t)g(y) dy=∫Rn[0]g(y) dy=0(x∈Rn,t>0)\begin{equation*}u_t(x,t) - \Delta u(x,t) = \int_{\mathbb{R}^n} [\Phi_t - \Delta_x\Phi](x-y,t)g(y)\,\text{d}y = \int_{\mathbb{R}^n} [0]g(y)\,\text{d}y = 0 \quad (x \in \mathbb{R}^n, t > 0)\end{equation*}
Thus, we have confirmed that our constructed solution u(x,t)u(x,t) is indeed a solution to the heat equation.
3. Let us now prove that our solution satisfies the initial condition by showing that u(x,t)→g(x)u(x,t) \to g(x) as t→0+t \to 0^+. Fix any x0∈Rnx^0 \in \mathbb{R}^n and ε>0\varepsilon > 0. Since gg is continuous, there exists δ>0\delta > 0 such that
∣g(y)−g(x0)∣<εwhenever ∣y−x0∣<δ,  y∈Rn.\begin{equation*}|g(y) - g(x^0)| < \varepsilon \quad \text{whenever } |y-x^0| < \delta, \; y \in \mathbb{R}^n.\end{equation*}
For ∣x−x0∣<δ2|x-x^0| < \frac{\delta}{2}, we can decompose the difference between u(x,t)u(x,t) and g(x0)g(x^0) as follows:
∣u(x,t)−g(x0)∣=∣∫RnΦ(x−y,t)[g(y)−g(x0)] dy∣≤∫B(x0,δ)Φ(x−y,t)∣g(y)−g(x0)∣ dy⏟I+∫Rn∖B(x0,δ)Φ(x−y,t)∣g(y)−g(x0)∣ dy⏟J.\begin{equation*}|u(x,t) - g(x^0)| = \left|\int_{\mathbb{R}^n} \Phi(x-y,t)[g(y) - g(x^0)]\,\text{d}y\right| \leq \underbrace{\int_{B(x^0,\delta)} \Phi(x-y,t)|g(y) - g(x^0)|\,\text{d}y}_{I} + \underbrace{\int_{\mathbb{R}^n \setminus B(x^0,\delta)} \Phi(x-y,t)|g(y) - g(x^0)|\,\text{d}y}_{J}.\end{equation*}
For the first term II, since Φ(x−y,t)\Phi(x-y,t) is the density of a normal distribution with mean xx and variance tt:
∣g(y)−g(x0)∣<ε⟹I≤ε∫RnΦ(x−y,t) dy=ε\begin{equation*}|g(y) - g(x^0)| < \varepsilon \quad \Longrightarrow \quad I \leq \varepsilon \int_{\mathbb{R}^n} \Phi(x-y,t)\,\text{d}y = \varepsilon\end{equation*}
For the second term JJ, observe that when ∣x−x0∣≤δ2|x-x^0| \leq \frac{\delta}{2} and ∣y−x0∣≥δ|y-x^0| \geq \delta:
∣y−x0∣≤∣y−x∣+δ2≤∣y−x∣+12∣y−x0∣⟹∣y−x∣≥12∣y−x0∣\begin{equation*}|y-x^0| \leq |y-x| + \frac{\delta}{2} \leq |y-x| + \frac{1}{2}|y-x^0| \quad \Longrightarrow \quad |y-x| \geq \frac{1}{2}|y-x^0|\end{equation*}
This leads to
J≤2∥g∥L∞∫Rn∖B(x0,δ)Φ(x−y,t) dy≤Ctn/2∫Rn∖B(x0,δ)e−∣x−y∣24t dy≤Ctn/2∫Rn∖B(x0,δ)e−∣y−x0∣216t dy=C∫Rn∖B(0,δ/t)e−∣z∣216 dz\begin{equation*}J \leq 2\|g\|_{L^\infty} \int_{\mathbb{R}^n \setminus B(x^0,\delta)} \Phi(x-y,t)\,\text{d}y \leq \frac{C}{t^{n/2}} \int_{\mathbb{R}^n \setminus B(x^0,\delta)} e^{-\frac{|x-y|^2}{4t}}\,\text{d}y \leq \frac{C}{t^{n/2}} \int_{\mathbb{R}^n \setminus B(x^0,\delta)} e^{-\frac{|y-x^0|^2}{16t}}\,\text{d}y = C\int_{\mathbb{R}^n \setminus B(0,\delta/\sqrt{t})} e^{-\frac{|z|^2}{16}}\,\text{d}z\end{equation*}
The final integral approaches zero as t→0+t \to 0^+. Therefore, for sufficiently small t>0t > 0 and ∣x−x0∣<δ2|x-x^0| < \frac{\delta}{2}, we have:
∣u(x,t)−g(x0)∣<2ε\begin{align*}|u(x,t) - g(x^0)| < 2\varepsilon\end{align*}
So, the theorem is proved □\Box.

The foolowing material is taken from [EvansPDE2010] p. 51-54

Illustration of ΓT\Gamma_T and UU
DEFINITION
Assume U⊂RnU \subset \mathbb{R^n} is open and bounded region, and fix a time T>0T > 0.
  • We define the parabolic cylinder UT:=U×[0,T].U_T := U \times [0,T].
  • The parabolic boundary of UTU_T is ΓT:=(U‾×{0})∪(∂U×[0,T]).\Gamma_T := (\overline{U}\times\{0\}) \cup (\partial U\times[0,T]).
We want next to derive a kind of analogue to the mean-value property for harmonic functions. There is no such simple formula. However let us observe that for fixed xx the spheres ∂B(x,r)\partial B(x,r) are level sets of the fundamental solution Φ(x−y)\Phi(x-y) for Laplace's equation. This suggests that perhaps for fixed (x,t)(x,t) the level sets of fundamental solution Φ(x−y,t−s)\Phi(x - y, t - s) for the heat equation may be relevant.
We write u∈C2,1(UT)u \in C^{2,1}(U_T) if uu has two continuous spatial derivatives and one continuous time derivative on UTU_T.
MEAN-VALUE PROPERTY OF HEAT EQUATION
Let u∈C2,1(UT)u \in C^{2,1}(U_T) solve the heat equation. For fixed x∈Rnx \in \mathbb{R}^n, t∈Rt \in \mathbb{R}, and r>0r > 0, define
E(x,t;r):={(y,s)∈Rn+1 ∣ s≤t, Φ(x−y,t−s)≥1rn}.\begin{equation*}E(x,t;r) := \left\{(y,s) \in \mathbb{R}^{n+1} \,|\, s \leq t, \,\Phi(x-y,t-s) \geq \frac{1}{r^n} \right\}.\end{equation*}
Then
u(x,t)=14rn∬E(x,t;r)u(y,s)∣x−y∣2(t−s)2 dy ds\begin{equation*}u(x,t) = \frac{1}{4r^n} \iint_{E(x,t;r)} u(y,s) \frac{|x-y|^2}{(t-s)^2} \,dy\,ds\end{equation*}
for each E(x,t;r)⊂UTE(x,t;r) \subset U_T.
Formula above is a sort of analogue for the heat equation of the mean-value formulas for Laplace's equation. Observe that the right-hand side involves only u(y,s)u(y,s) for times s≤ts \leq t. This is reasonable, as the value u(x,t)u(x,t) should not depend upon future times.
\quad Proof sketch. Shift the space and time coordinates so that x=0x = 0 and t=0t = 0. Upon mollifying if necessary, we may assume uu is smooth. Write E(r)=E(0,0;r)E(r) = E(0,0;r) and set
ϕ(r):=1rn∬E(r)u(y,s)∣y∣2s2 dy ds=∬E(1)u(ry,r2s)∣y∣2s2 dy ds.\begin{equation*}\phi(r) := \frac{1}{r^n} \iint_{E(r)} u(y,s) \frac{|y|^2}{s^2} \,dy\,ds= \iint_{E(1)} u(ry,r^2s) \frac{|y|^2}{s^2} \,dy\,ds.\end{equation*}
We compute
ϕ′(r)=∬E(1)∑i=1nuyiyi∣y∣2s2+2rus∣y∣2s dy ds=1rn+1∬E(r)∑i=1nuyiyi∣y∣2s2⏟A+2us∣y∣2s⏟B dy ds.\begin{equation*}\phi'(r) = \iint_{E(1)} \sum_{i=1}^n u_{y_i}y_i\frac{|y|^2}{s^2} + 2ru_s\frac{|y|^2}{s} \,dy\,ds = \frac{1}{r^{n+1}} \iint_{E(r)} \underbrace{\sum_{i=1}^n u_{y_i}y_i\frac{|y|^2}{s^2}}_{A} + \underbrace{2u_s\frac{|y|^2}{s}}_{B} \,dy\,ds.\end{equation*}
To handle the term BB, we introduce a function that vanishes on ∂E(r)\partial E(r):
ψ:=−n2log⁡(−4πs)+∣y∣24s+nlog⁡r\begin{equation*}\psi := -\frac{n}{2}\log(-4\pi s) + \frac{|y|^2}{4s} + n\log r\end{equation*}
and observe ψ=0\psi = 0 on ∂E(r)\partial E(r), since Φ(y,−s)=r−n\Phi(y,-s) = r^{-n} on ∂E(r)\partial E(r). We utilize ψ\psi to write
B=1rn+1∬E(r)4us∑i=1nyiψyi dy ds=−1rn+1∬E(r)4nusψ+4∑i=1nusyiψ dy ds;\begin{equation*}B = \frac{1}{r^{n+1}} \iint_{E(r)} 4u_s \sum_{i=1}^n y_i\psi_{y_i} \,dy\,ds = -\frac{1}{r^{n+1}} \iint_{E(r)} 4nu_s\psi + 4\sum_{i=1}^n u_{sy_i}\psi \,dy\,ds;\end{equation*}
there is no boundary term since ψ=0\psi = 0 on ∂E(r)\partial E(r). Integrating by parts with respect to ss, we discover
B=1rn+1∬E(r)−4nusψ+4∑i=1nuyiψs dy ds=1rn+1∬E(r)−4nusψ+4∑i=1nuyiyi(−n2s−∣y∣24s2) dy ds=1rn+1∬E(r)−4nusψ−2ns∑i=1nuyiyi dy ds−A.\begin{equation*}B = \frac{1}{r^{n+1}} \iint_{E(r)} -4n u_s\psi + 4\sum_{i=1}^n u_{y_i}\psi_s \,dy\,ds = \frac{1}{r^{n+1}} \iint_{E(r)} -4n u_s\psi + 4\sum_{i=1}^n u_{y_i}y_i\left(-\frac{n}{2s} - \frac{|y|^2}{4s^2}\right) \,dy\,ds = \frac{1}{r^{n+1}} \iint_{E(r)} -4n u_s\psi - \frac{2n}{s}\sum_{i=1}^n u_{y_i}y_i \,dy\,ds - A.\end{equation*}
Consequently, since uu solves the heat equation,
ϕ′(r)=A+B=1rn+1∬E(r)−4nΔuψ−2ns∑i=1nuyiyi dy ds=∑i=1n1rn+1∬E(r)4nuyiψyi−2nsuyiyi dy ds=0.\begin{equation*}\phi'(r) = A + B = \frac{1}{r^{n+1}} \iint_{E(r)} -4n\Delta u\psi - \frac{2n}{s}\sum_{i=1}^n u_{y_i}y_i \,dy\,ds = \sum_{i=1}^n \frac{1}{r^{n+1}} \iint_{E(r)} 4nu_{y_i}\psi_{y_i} - \frac{2n}{s}u_{y_i}y_i \,dy\,ds = 0.\end{equation*}
Thus ϕ\phi is constant, and therefore
ϕ(r)=lim⁡t→0ϕ(t)=u(0,0)(lim⁡t→01tn∬E(t)∣y∣2s2 dy ds)=4u(0,0),\begin{equation*}\phi(r) = \lim_{t\to 0}\phi(t) = u(0,0)\left(\lim_{t\to 0}\frac{1}{t^n}\iint_{E(t)} \frac{|y|^2}{s^2}\,dy\,ds\right) = 4u(0,0),\end{equation*}
where the last constant comes from the scaling identity
1tn∬E(t)∣y∣2s2 dy ds=∬E(1)∣y∣2s2 dy ds=4.\begin{equation*}\frac{1}{t^n}\iint_{E(t)} \frac{|y|^2}{s^2}\,dy\,ds = \iint_{E(1)} \frac{|y|^2}{s^2}\,dy\,ds = 4.\end{equation*}
The direct computation of the integral over E(1)E(1) is routine but somewhat long, so we omit it here. □\Box

References