Heat equation, apparently just for fun

It often happens in research that one writes a considerable amount of code only to find that it won’t turn out to be useful in whatever project one is working on. Such was the case with me today: I wrote a solver for the diffusion equation only to find that it won’t suit the purpose for which I originally wrote it. Oh well! Check out these cool numerical solutions anyway.

Recall that the diffusion equation with constant rate of diffusivity is given by

My application required me to solve it in two dimensions only, so that’s what I’ll do here, setting the domain of solution \(\Omega = [0,1]\times[0,1]\). I considered two basic kinds of boundary conditions:

  1. Neumann, given by \( \nabla_{n}\rho = 0 \) (or, more generally, some function of time), where \(n\) is a vector normal to the boundary \(\partial \Omega\).
  2. Dirichlet, given by \(\rho(\partial \Omega) = 0 \) (again, more generally, some function of time).

First, let’s see what happens when we start out with a multivariate normal distribution centered at \( \mu = (\frac{1}{2}, \frac{1}{2})^T \) with covariance matrix given by \( \Sigma = \frac{1}{4} I \) and use Neumann boundary conditions:

png

The left panel is the initial condition, while the right is the state of the system after \(N_t\) timesteps of size \(\Delta t\). Since I won’t be able to use this code for my research, I just made some fun shapes…

png

png

The above are both using the Neumann BCs. For a interesting example of how boundary conditions really matter in PDEs (even with ones as straightforward as these!), constrast the following two simulations. The first uses the Neumann BCs, the second uses Dirichlet.

png

png

Physically, we can interpret the Dirichlet conditions as saying that there is some energy source (or sink, in this case) that is constantly injecting (or withdrawing) energy from \(\Omega\) in order to keep the boundary a fixed temperature. The Neumann conditions, on the other hand, describe the case of a perfectly insulated boundary.

Update

We’d better just solve the Laplace equation too, while we’re at it. This is the steady state heat equation,

and here I’ll solve it with Dirichlet boundary conditions only. What follows isn’t even related to anything else I’m doing–100% of this is just for fun.

With the BCs \( \rho(x, 0) = \cos(2\pi x) + 1,\ \rho(x, 1) = \sin(2\pi x) + 1,\ \rho(0, y) = \rho(1, y) = 1\):

png

With the tent map on the upper and lower boundaries, and zero boundary conditions on the sides:

png

With the tent map on all sides:

png

And yes, just to reassure you that my solutions are accurate, here’s one that many will instinctually recognize: bottom and left walls held constant at one, top and right walls held constant at zero:

png

The code is here. Be careful when running it for many timesteps (e.g., over 100K); it stores the entire state of the system at every timestep in memory, so that can grow quickly. (I’ll end up changing this behavior eventually.)

N.B.: This was compiled with a new compiler!

Written on August 1, 2017