# Basic Solution to the Heat Equation

As a learning example, I am trying to derive the solution to the basic Heat Equation (https://en.wikipedia.org/wiki/Heat_equation) using Fourier Transforms.

As I understand, the Heat Equation can describe the following scenario: Suppose you have a 1 dimensional rod with constant temperatures at both ends of the rod. At time = 0, the initial distribution of heat across the length of the rod is $$g(x)$$. Assuming nothing else changes - how will this initial distribution of heat diffuse across the rod as time goes on?

The heat equation in one dimension is given by:

$$\frac{\partial u}{\partial t} = \alpha \frac{\partial^2 u}{\partial x^2}$$

Where $$u(x,t)$$ is the temperature at position $$x$$ and time $$t$$, and $$\alpha$$ is the thermal diffusivity.

Step 1: Here is the general formula for Fourier transform of some function $$u(x,t)$$ with respect to $$x$$ as (note: after the Fourier Transform, $$x$$ is replaced with $$k$$):

$$\hat{u}(k,t) = \int_{-\infty}^{\infty} u(x,t) e^{-ikx} dx$$

• Left Hand Side: Take the Fourier Transform with respect to $$x$$.

For some function $$U(x,t)$$ , the Fourier Transform with respect to x can be written as:

$$\hat{U}(x,t) = \mathcal{F}_{x} \{ U(x,t) \}$$

We can see that a Fourier Transform with respect to $$x$$ has no implication on $$t$$ ; thus we can take $$t$$ out of the integral:

$$\mathcal{F} \left( \frac{\partial u}{\partial t} \right) = \int_{-\infty}^{\infty} \left( \frac{\partial U(x,t)}{\partial t} \right) e^{-ikx} \, dx$$

$$\mathcal{F} \left( \frac{\partial u}{\partial t} \right) = \frac{\partial}{\partial t} \int_{-\infty}^{\infty} U(x,t) e^{-ikx} \, dx$$

Thus, for the Left Hand Side of the Heat Equation - the Fourier Transform with respect to $$x$$ can be wrriten as: $$\mathcal{F}\left(\frac{\partial u}{\partial t}\right) = \frac{\partial}{\partial t}\mathcal{F}(u) = \frac{\partial \hat{u}}{\partial t}$$

• Right Hand Side (again, Fourier Transform with respect to $$x$$), we use the property that the Fourier transform of the second derivative is $$-k^2$$ times the Fourier transform of the function: $$\mathcal{F}\left(\alpha\frac{\partial^2 u}{\partial x^2}\right) = \alpha(-k^2)\hat{u} = -\alpha k^2 \hat{u}$$

• Now, we put the Left Hand Side and the Right Hand Side Together: $$\frac{\partial \hat{u}}{\partial t} = -\alpha k^2 \hat{u}$$ $$\frac{\partial}{\partial t} \hat{u}(k,t) = -\alpha k^2 \hat{u}(k,t)$$

After applying the Fourier Transform to both sides, we see that the new equation is only a first order differential equation whereas the original Heat Equation was a second order differential equation. Furthermore, in the original Heat Equation, there were partial derivatives with respect to both $$t$$ and $$x$$. After taking the Fourier Transform, there is only a partial derivative with respect to $$t$$. Thus it appears that the Fourier Transform served to simplify the complexity of the original Heat Equation

Step 2: Solving for $$\hat{u}(k,t)$$:

We now have a first-order ordinary differential equation in $$t$$: $$\frac{\partial \hat{u}}{\partial t} = -\alpha k^2 \hat{u}$$

To solve this, we can separate variables: $$\frac{1}{\hat{u}} \frac{\partial \hat{u}}{\partial t} = -\alpha k^2$$

Integrating both sides with respect to $$t$$ (I think $$dt$$ and $$\partial t$$ can cancel each other out): $$\int \frac{1}{\hat{u}} \frac{\partial \hat{u}}{\partial t} dt = \int -\alpha k^2 dt$$ $$\ln(\hat{u}) = -\alpha k^2 t + C$$

Taking the exponential of both sides (let $$A = e^C$$ ): $$\hat{u} (k,t) = e^{-\alpha k^2 t + C} = e^{-\alpha k^2 t} \cdot e^{C} = Ae^{-\alpha k^2 t}$$

At $$t=0$$, we can see that $$\hat{u}(k,0) = A$$:

$$\hat{u} (k,0) = Ae^{-\alpha k^2 \cdot 0} = Ae^{0} = A$$

Therefore : $$\hat{u}(k,t) = \hat{u}(k,0)e^{-\alpha k^2 t} = Ae^{-\alpha k^2 t}$$

Notice that at this stage, the above equation is no longer a differential equation!

Step 3: Apply the initial condition. Let's assume that across the rod at time =0, the initial distribution of heat is $$g(x)$$:

$$U(x,0) = g(x)$$

Now, we take the Fourier Transform (with respect to $$x$$) of $$U(x,0) = g(x)$$.

Previously, we showed that at time=0:

$$\hat{u} (k,0) = Ae^{-\alpha k^2 \cdot 0} = A$$

Now we approach time=0 from the other side - since we are assuming a heat distribution of $$g(x)$$ at time=0, we can relate these two ideas together:

$$\hat{u} (k,0) = F_x[g(x)] = Ae^{-\alpha k^2 \cdot 0} = A$$

This shows us that there is a relationship between the constant of integration $$A$$ and the initial conditions $$g(x)$$. In other words, the Fourier Transform of $$g(x)$$:

$$F_x[g(x)] = g(k) = A$$

Step 5: Inverse Fourier Transform

In the previous steps, the Fourier Transform served to reduce the complexity of the Heat Equation. Now, after applying the Fourier Transform - we need to bring it back in terms of the original variable (i.e. we went from $$t$$ to $$k$$ - now we have to go from $$k$$ to $$t$$)

The inverse Fourier transform is defined as: $$f(x) = \frac{1}{2\pi} \int_{-\infty}^{\infty} \hat{f}(k) e^{ikx} dk$$

Now, we apply the inverse Fourier transform to our problem for the goal of finding out$$u(x,t)$$. The inverse Fourier transform in our problem can be defined as:

$$u(x,t) = \frac{1}{2\pi} \int_{-\infty}^{\infty} \hat{u}(k,t) e^{ikx} \, dk.$$

Recall that: $$\hat{u}(k,t) = \hat{u}(k,0)e^{-\alpha k^2 t} = Ae^{-\alpha k^2 t} = g(k) e^{-\alpha k^2 t}$$

Substitute $$\hat{u}(k,t) = \hat{g}(k) e^{-\alpha k^2 t}$$:

$$u(x,t) = \frac{1}{2\pi} \int_{-\infty}^{\infty} \hat{g}(k) e^{-\alpha k^2 t} e^{ikx} \, dk.$$

Now, we are back in terms of the original variables of the Heat Equation.

Step 5: Solve the integral - I struggled a lot with this.

I first wrote the Fourier Transform of the Initial Condition:

$$\hat{g}(k) = \int_{-\infty}^{\infty} g(y) e^{-iky} dy$$

Then I substituted this into the original integral:

$$u(x,t) = \frac{1}{2\pi} \int_{-\infty}^{\infty} \left(\int_{-\infty}^{\infty} g(y) e^{-iky} dy\right) e^{-\alpha k^2 t} e^{ikx} \, dk$$

I have often seen that the order of integrals can be rearranged - I tried to do this here (Fubini's theorem):

$$u(x,t) = \frac{1}{2\pi} \int_{-\infty}^{\infty} g(y) \left(\int_{-\infty}^{\infty} e^{-\alpha k^2 t} e^{ik(x-y)} \, dk\right) dy$$

I read references online that said the internal integral here is based on a Gaussian Integral (https://en.wikipedia.org/wiki/Gaussian_integral)

$$I = \int_{-\infty}^{\infty} e^{-\alpha k^2 t} e^{ik(x-y)} \, dk$$

I tried to re-arrange the terms in the exponent:

$$= \int_{-\infty}^{\infty} e^{-\alpha t(k^2 - \frac{i(x-y)}{\alpha t}k)} \, dk$$ $$= e^{-\frac{(x-y)^2}{4\alpha t}} \int_{-\infty}^{\infty} e^{-\alpha t(k^2 - \frac{i(x-y)}{\alpha t}k + (\frac{i(x-y)}{2\alpha t})^2)} \, dk$$

$$I = e^{-\frac{(x-y)^2}{4\alpha t}} \int_{-\infty}^{\infty} e^{-\alpha t(k - \frac{i(x-y)}{2\alpha t})^2} \, dk$$

We still need to do some further manipulations to exploit the Gaussian Integral formula.

Suppose we define a new variable: $$u = k - \frac{i(x-y)}{2\alpha t}$$

Then, $$dk = du$$, and when $$k$$ goes from $$-\infty$$ to $$\infty$$, $$u$$ also goes from $$-\infty$$ to $$\infty$$. The integral becomes:

$$I = e^{-\frac{(x-y)^2}{4\alpha t}} \int_{-\infty}^{\infty} e^{-\alpha t u^2} \, du$$

Now we have a standard Gaussian integral. The general form of this integral is:

$$\int_{-\infty}^{\infty} e^{-a x^2} \, dx = \sqrt{\frac{\pi}{a}}$$

In our case, $$a = \alpha t$$, so:

$$\int_{-\infty}^{\infty} e^{-\alpha t u^2} \, du = \sqrt{\frac{\pi}{\alpha t}}$$

Substituting this result back into our equation:

$$I = e^{-\frac{(x-y)^2}{4\alpha t}} \sqrt{\frac{\pi}{\alpha t}}$$

$$I = e^{-\frac{(x-y)^2}{4\alpha t}} \sqrt{\frac{\pi}{\alpha t}}$$

Substituting this back into our original equation:

$$u(x,t) = \frac{1}{2\pi} \int_{-\infty}^{\infty} g(y) \left(e^{-\frac{(x-y)^2}{4\alpha t}} \sqrt{\frac{\pi}{\alpha t}}\right) dy$$

$$u(x,t) = \frac{1}{\sqrt{4\pi\alpha t}} \int_{-\infty}^{\infty} g(y) e^{-\frac{(x-y)^2}{4\alpha t}} dy$$

This somewhat resembles a Gaussian Distribution as to be expected. It seems that this integral will likely need to be numerically solved (depending on the choice of $$g(x)$$). I think the gaussian here is not a probability distribution, but a weighting function that determines how much the initial temperature at point y contributes to the temperature at point x at time t....

• While the final answer is possibly correct - are my steps correct?
• I am also struggling to understand the following point : A gaussian distribution is random but the heat equation is deterministic. Does this mean that this heat has a random component according to this equation?
• Your second question is, essentially, the core concept of the central limit theorem (or, really, Donsker's, in this context). Commented Jul 16 at 19:05
• @XanderHenderson you are right. The answers below do not mention that. I think this question focus on why gaussian distrubution shares same integral transformation kernel with fourier transformation, instead of considering it from a physical perspective. Commented Jul 17 at 3:18

As of your final answer, it is indeed a correct one. In more general terms, you would need to find a so-called fundamental solution $$\Phi(x,t)$$, which is the solution to the problem: $$\left\{ \begin{array}{} \dfrac{\partial u}{\partial t} &=& \alpha\dfrac{\partial^2 u}{\partial t^2} \\ u(x, 0) &=& \delta(x) \end{array} \right.$$ where $$\delta(x)$$ is a Dirac delta function. Then, for some arbitrary initial condition $$u(x, 0)=g(x)$$, the solution could be obtained by using convolution of the fundamental solution and the initial condition: $$u(x, t) = \Phi(x, t) * g(x) = \int_{\mathbb{R}}\Phi(x-y, t)\,g(y)\,dy$$

And, as it turns out, the fundamental solution for one-dimensional heat equation is given by: $$\Phi (x,t)=\dfrac {1}{\sqrt {4\pi \alpha t}}\exp \left(-{\frac {x^{2}}{4 \alpha t}}\right)$$

As for the question:

I am also struggling to understand the following point: a gaussian distribution is random but the heat equation is deterministic. Does this mean that this heat has a random component according to this equation?

You may think of the heat conduction as particles that get hotter move faster and hit other particles, warming them up. That process is indeed stochastic, since you can't know for sure which direction a particular particle should move the next instant and how much energy per se it should lose or gain. However, the total process at a macroscopic level is quite measurable.

As a side note, you do know, probably, that heat conducion and diffusion process are described by the same equation. What do you think it could imply in terms of what I mentioned above?

There is a conceptual bridge between stochastic models of motion of particles by random steps in time and the smooth Gaussian as the probability density of the final distribution the the position variable of a stochastic process $$t\mapsto X_t(\omega)_{\omega \in \Omega}\quad |\quad X_0=0.$$

Physically Brownian motion (first understood by Einstein) has been mathematically modelled as the Wiener process (invented by Wiener, making proof sound by Ito).

Physically, the central point is to understand the diffusion constant $$D = 10^{-5}\ \frac{cm^2}{s}$$ for the diffusion equation

$$\partial_t \rho(t,x)[\frac{1}{\text{cm}^3}] \ = \ D[\frac{cm^2}{s}] \partial_x^2 \ \rho(t,x)[\frac{1}{\text{cm}^3}]$$ of a collection of small particles in water at a given temperature $$k T$$ and mass $$m$$, that are pushed around at random by its neighbourhood of thermally excited water molecules with a Gaussian momentum distribution (aka Maxwell distribution of velocity or exponential distribution of kinetic energy, defining absolute temperature by $$e^{-\frac{m \ v^2}{2 k \ T}}$$) .

Assume a single particle moving according to discrete time steps $$dt$$ in one space dimension $$X(t+dt, \omega_t) = X(t, \omega_t) + \sigma\ \sqrt{dt}\ G(\omega_t) | X(0,\omega_0)$$ where G is a $$N(0,1)$$ distributed gaussian random variable with $$\mathbb E(G(\omega_t) G(\omega_s))= \delta(g-s)$$ along a specific random path generated by the event family $$\omega\in \Omega$$ , the final position is a sum of statistical identically distributed, independent Gaussian variables with normal distribution $$N(0, \sigma^2 dt )$$.

Under a few assumptions any long sum of $$n$$ iid variables with $$\mathbb E X=0, \mathbb E X^2= \sigma^2 dt$$ converges to a Gaussian, distributed as $$N(0, \sigma^2 \ n \ \text{dt})$$

So we don't need to bother, that the solution of the diffusion equation in one dimension as the probability density of that process has to be

$$\rho(t,x) = \frac{1}{\sqrt{2 \ \pi \ \sigma^2 \ t }} \ e^{-\frac{x^2}{2 \sigma^2 \ t}}$$

Four tests of model reliability:

Total probability 1 and independent of time $$t$$ $$\int_{-\infty}^\infty \ \rho(t,x) \ dx \ = \ 1,$$ solves the diffuison equation by $$\partial_t \rho = \ \rho \ * \ \left(\ \frac{x^2}{2\sigma ^2 t^2}\ - \frac{1}{2 t}\right), \quad \partial_x \rho(t,x) \ = \ - \rho(t,x) \ \frac{x}{\sigma^2 \ t},$$ converges to point measure at $$t\to +0,$$ $$\lim_{t\to +0}\ \int_{-\infty}^\infty f(x) \rho(t,x-y) dx =\lim_{t\to +0}\ \int_{-\infty}^\infty f(x+y) \rho(t,x) dx = \int_{-\infty}^\infty f(x+y) \delta(x) dx = f(y),$$ and its convolution with any start point distribution solves the diffusion equation, too: $$(\partial_t - \sigma^2 \partial_x^2) \ \int_{-\infty}^\infty f(x) \rho(t,x-y) dx = \int_{-\infty}^\infty f(x+y) (\partial_t - \sigma^2 \partial_x^2)\rho(t,x) \ = \ 0$$

Bottom line: The diffusion of an (nearly) infinite collection of independently Brownian-moving particles, starting with a given density, results in the distribution of Brwonian paths with the given start distribution.

The fact that the Gaussian is an exponential of the euclidean 2-norm, with the norm being a product, too, represents the probability density as a product over the dimension, nicely fitting with independency of the stochastic movement in $$n$$-dimensional euclidean space in a cartesian frame.

In the diffusion equation it suffices to replace $$(\partial_t - \sigma^2 \ \partial_x^2) \to (\partial_t - \sigma^2 \Delta_{\mathbb x})$$

As always in pure mathematics, the axiomatic representation takes these simple things as a known fact and concentrates on the economy of its proofs from the set of axioms.