Today, we will use Python to analytically solve one of the most important partial differential equations out there, the diffusion equation. It is a fundamental equation that arises in many areas of physics, chemistry, biology, and engineering, and has an enormous range of applications. For instance, in physics, the diffusion equation is used to describe the movement of heat through a solid or fluid, as well as the transport of particles in gases and liquids. Also, the Schrödinger equation is (in principle) a kind of diffusion equation. In biology, the diffusion equation is used to study the movement of nutrients and other substances through living tissues, as well as the spread of diseases. In engineering, the diffusion equation is used to design and optimize processes such as heat exchangers, catalytic converters, and solar cells. It is even applied in economics, where the diffusion equation can be used to model the spread of ideas or innovations through a population. Due to its wide range of applications, it is not surprising that it was found independently by several scientists of the 19th century working on completely different topics. The famous George Gabriel Stokes derived the equation in 1845 to describe the motion of fluids, while the even more famous James Clerk Maxwell derived the equation in 1860 to describe the motion of heat.

The equation is given by:

where t is time, 𝛼 a constant, and ∇² is the Laplacian operator. For concreteness, we will solve the following example application. We have a 10cm long bar with insulated sides initially at 100° everywhere. Starting at 𝑡=0, the ends are held at temperature 0°. Our task is to find the temperature distribution at some later time 𝑡.

The method we will use is the separation of variables, i.e. we use the ansatz

where 𝑇 and 𝑋 are functions of a single variable 𝑡 and 𝑥, respectively. Start a new Jupyter notebook and define the symbols for our problem:

Insert the separation ansatz:

Divide by 𝑇(𝑡)𝑋(𝑥):

The key observation here is that the left side depends only on t, whereas the right side depends only on x. So, when you change 𝑡, only the left side changes but not the right side. Conversely, if you change 𝑥, the right side changes, but not the left side. The only way how this can be true is that both sides are equal to a constant. Let’s call that constant −𝑘²:

This is just an ordinary differential equation that we can easily solve:

Now focus on the side of the separated equation:

Again, a simple ODE:

Hmmm… sympy uses the same names of the constants of integration, so we get a name clash. Let’s rename the constants:

It’s more convenient to have that in terms of sines and cosines, so we can rewrite that:

The coefficients are independent, so we can simply give new names to them:

And that’s the form we will work with.

Now consider the first boundary conditions at 𝑥=0. This gives a value for the constant 𝐵:

The other boundary condition gives a constraint on 𝑘:

Now we can assemble a solution from our base solutions.

where we have absorbed the coefficient 𝐶 into the new coefficients 𝑎_𝑛.

But this must still satisfy the initial condition, which we haven’t evaluated yet. So do that now:

We can solve for 𝑎_𝑛 by using the fact that the sines are orthogonal. So multiply both sides by a suitable sine and integrate:

The infinite sum on the left side has only one non-zero term, for which 𝑚=𝑛:

Now insert the values for the coefficients in our solution:

And that is the exact solution to our problem.

Finally, we can plot it. But of course, we have to cut off the series at some point. We will approximate the solution by the first 30 terms:

Now we can plot the solution with the sympy plotting backends spb :