Diffusion equation

So far, we have explored ordinary differential equations where physical quantities depend only on one variable (usually time). However, many physical phenomena involve quantities that change in both space and time simultaneously. These situations require partial differential equations (PDEs). The diffusion equation is one of the most important PDEs in physics, appearing in heat conduction, particle diffusion, and even quantum mechanics (as the time-dependent Schrödinger equation).

Physical Model

Understanding Diffusion

Imagine dropping a drop of ink into a glass of still water. Initially, all the ink particles are concentrated in one small region. Over time, you’ll observe the ink spreading outward, gradually coloring the entire glass of water. This spreading process is called diffusion.

At the microscopic level, diffusion results from the random motion of particles (Brownian motion). Each ink particle undergoes a “random walk” - constantly changing direction due to collisions with water molecules. While individual particle motion is unpredictable, the collective behavior follows predictable statistical laws.

The animation below demonstrates this process: particles start concentrated at the center and gradually spread throughout the domain due to their random motion.

width = 600
height = 600
margin = ({top: 20, right: 30, bottom: 20, left: 40})
plotHeight = 100

viewof simulation = {
  // Create main container
  const container = d3.create("div")
    .style("display", "flex")
    .style("flex-direction", "column");

  // Create SVG for particle simulation
  const svg = container.append("svg")
    .attr("width", width)
    .attr("height", height - plotHeight)
    .attr("viewBox", [0, 0, width, height - plotHeight]);

  // Create SVG for histogram
  const histogramSvg = container.append("svg")
    .attr("width", width)
    .attr("height", plotHeight)
    .attr("viewBox", [0, 0, width, plotHeight]);

  const numParticles = 1000;
  const D = 0.5; // Diffusion coefficient
  const numBins = 80;

  // Create particles at the center
  const particles = Array.from({length: numParticles}, () => ({
    x: width / 2,
    y: (height - plotHeight) / 2,
    vx: 0,
    vy: 0
  }));

  // Setup scales for histogram
  const xScale = d3.scaleLinear()
    .domain([0, width])
    .range([margin.left, width - margin.right]);

  const yScale = d3.scaleLinear()
    .domain([0, numParticles/5])
    .range([plotHeight - margin.bottom, margin.top]);

  // Create histogram generator
  const histogram = d3.bin()
    .domain(xScale.domain())
    .thresholds(xScale.ticks(numBins))
    .value(d => d.x);

  // Create histogram group
  const histogramGroup = histogramSvg.append("g");

  // Add axes
  histogramSvg.append("g")
    .attr("transform", `translate(0,${plotHeight - margin.bottom})`)
    .call(d3.axisBottom(xScale));

  // Animation function
  function animate() {
    particles.forEach(particle => {
      // Random walk implementation based on diffusion equation
      const randomAngle = Math.random() * 2 * Math.PI;
      const displacement = Math.sqrt(2 * D);

      particle.x += displacement * Math.cos(randomAngle);
      particle.y += displacement * Math.sin(randomAngle);

      // Bounce off walls
      if (particle.x < 0) particle.x = 0;
      if (particle.x > width) particle.x = width;
      if (particle.y < 0) particle.y = 0;
      if (particle.y > (height - plotHeight)) particle.y = height - plotHeight;
    });

    // Update particle positions
    circles
      .attr("cx", d => d.x)
      .attr("cy", d => d.y);

    // Update histogram
    const bins = histogram(particles);

    const bars = histogramGroup.selectAll("rect")
      .data(bins);

    bars.enter()
      .append("rect")
      .merge(bars)
      .attr("x", d => xScale(d.x0))
      .attr("y", d => yScale(d.length))
      .attr("width", d => Math.max(0, xScale(d.x1) - xScale(d.x0) - 1))
      .attr("height", d => yScale(0) - yScale(d.length))
      .attr("fill", "steelblue");

    bars.exit().remove();
  }

  // Create circles for particles
  const circles = svg.selectAll("circle")
    .data(particles)
    .join("circle")
    .attr("r", 2)
    .attr("fill", "steelblue")
    .attr("opacity", 0.6);

  // Start animation
  d3.timer(animate);

  return container.node();
}

Figure 1: Simulation showing the diffusion of particles in a 2D box. Particles move randomly based on the diffusion equation. The histogram shows the distribution of particles along the x-axis.

The Diffusion Equation

Rather than tracking individual particles, we work with concentration - the number of particles per unit volume at each point in space. The diffusion equation describes how this concentration changes over time:

\[\begin{equation} \frac{\partial c(\mathbf{r},t)}{\partial t}=D\nabla^2 c (\mathbf{r},t) \end{equation}\]

Where:

\(c(\mathbf{r},t)\) is the concentration at position \(\mathbf{r}\) and time \(t\)
\(D\) is the diffusion coefficient (units: length²/time)
\(\nabla^2\) is the Laplacian operator (measures the “curvature” of the concentration)

The diffusion coefficient \(D\) characterizes how quickly diffusion occurs - larger \(D\) means faster spreading.

Simplifying to One Dimension

To make our numerical solution manageable, we’ll consider diffusion along a single direction (say, the x-axis). This reduces our equation to:

\[\begin{equation} \frac{\partial c(x,t)}{\partial t}=D\frac{\partial^2 c(x,t)}{\partial x^{2}} \end{equation}\]

This equation tells us: The rate of change of concentration at any point equals the diffusion coefficient times the curvature of the concentration profile at that point.

From Continuous to Discrete

To solve this equation numerically, we must discretize both space and time - converting our continuous problem into a discrete one that computers can handle.

Spatial Discretization

We divide our spatial domain into a grid of points separated by distance \(\Delta x\). Instead of knowing \(c(x,t)\) at every point \(x\), we’ll track concentrations at specific grid points: \(c_0, c_1, c_2, \ldots, c_N\).

The second spatial derivative (curvature) at point \(i\) can be approximated using three neighboring points:

\[\begin{equation} \frac{\partial^{2} c(x,t)}{\partial x^2}\bigg|_{x=x_i} \approx\frac{c_{i+1}^{n}-2c_{i}^{n}+c_{i-1}^{n}}{\Delta x^2} \end{equation}\]

This finite difference formula captures the essential physics: if the concentration at point \(i\) is higher than the average of its neighbors, the curvature is negative, and diffusion will tend to reduce the concentration at point \(i\).

Matrix Representation

We can represent this spatial discretization as a matrix operation. For our concentration vector \(\mathbf{C} = [c_0, c_1, c_2, \ldots, c_N]^T\), the spatial derivative becomes:

\[\frac{\partial^2 \mathbf{C}}{\partial x^2} = \frac{1}{\Delta x^2} \mathbf{M} \mathbf{C}\]

where \(\mathbf{M}\) is the second-difference matrix:

\[\mathbf{M} = \begin{bmatrix} -2 & 1 & 0 & 0 & \cdots & 0\\ 1 & -2 & 1 & 0 & \cdots & 0\\ 0 & 1 & -2 & 1 & \cdots & 0\\ \vdots & \vdots & \ddots & \ddots & \ddots & \vdots\\ 0 & \cdots & 0 & 1 & -2 & 1\\ 0 & \cdots & 0 & 0 & 1 & -2\\ \end{bmatrix}\]

Each row of this matrix implements the finite difference stencil for one grid point.

Dirichlet Boundary Conditions (Fixed Values)

The matrix above assumes so-called Dirichlet boundary conditions, where the concentration is fixed at specific values at the boundaries (typically zero). This represents situations where particles are absorbed when they reach the boundaries.

Physical example: A pipe with absorbing walls where particles disappear upon contact.

Matrix structure: The first and last rows have the standard -2, 1 pattern, enforcing \(c_{-1} = 0\) and \(c_{N+1} = 0\).

Neumann Boundary Conditions (Zero Flux)

For open boundary conditions (zero flux: \(\frac{\partial c}{\partial x} = 0\) at boundaries), the matrix would be modified:

\[\mathbf{M}_{\text{Neumann}} = \begin{bmatrix} -1 & 1 & 0 & 0 & \cdots & 0\\ 1 & -2 & 1 & 0 & \cdots & 0\\ 0 & 1 & -2 & 1 & \cdots & 0\\ \vdots & \vdots & \ddots & \ddots & \ddots & \vdots\\ 0 & \cdots & 0 & 1 & -2 & 1\\ 0 & \cdots & 0 & 0 & 1 & -1\\ \end{bmatrix}\]

These boundary conditions are now called Neumann boundary conditions.

Physical example: A system with reflective or impermeable boundaries where no material can cross.

Matrix changes: First and last diagonal entries change from -2 to -1.

Temporal Discretization: The Crank-Nicolson Scheme

For time discretization, we use the Crank-Nicolson scheme, which provides excellent stability and accuracy. The key idea is to evaluate the spatial terms at the average of two consecutive time steps:

\[\begin{equation} \frac{c_i^{n+1} - c_i^n}{\Delta t} = \frac{D}{2}\left[\frac{\partial^2 c}{\partial x^2}\bigg|^{n+1} + \frac{\partial^2 c}{\partial x^2}\bigg|^{n}\right] \end{equation}\]

Rearranging this equation and using our matrix notation:

\[\begin{equation} \left(\mathbf{I}-\frac{\Delta t \cdot D}{2} \mathbf{M} \right )\mathbf{C}^{n+1}=\left(\mathbf{I}+\frac{\Delta t \cdot D}{2} \mathbf{M} \right )\mathbf{C}^{n} \end{equation}\]

This gives us our solution algorithm: \[\begin{equation} \mathbf{A}\mathbf{C}^{n+1}=\mathbf{B}\mathbf{C}^{n} \end{equation}\]

where \(\mathbf{A}=\left(\mathbf{I}-\frac{\Delta t \cdot D}{2} \mathbf{M} \right )\) and \(\mathbf{B}=\left(\mathbf{I}+\frac{\Delta t \cdot D}{2} \mathbf{M} \right )\).

Numerical Solution

Now we implement our numerical scheme in Python. We need to specify both boundary conditions (what happens at the edges) and initial conditions (the starting concentration profile).

Problem Setup

We’ll solve diffusion on a domain of length \(L=1\) with Dirichlet boundary conditions: the concentration is fixed at zero at both ends. This might represent a situation where particles are absorbed when they reach the boundaries.

For our initial condition, we’ll use a Gaussian (bell-shaped) distribution centered at \(x=L/2\):

\[\begin{equation} c(x,0)=\frac{1}{\sigma\sqrt{2\pi }}e^{-\frac{(x-L/2)^2}{2\sigma^2}} \end{equation}\]

Parameter Choices

The parameter selection in this diffusion equation implementation involves several carefully considered choices:

Spatial Parameters:

Domain length (L = 1.0): A unit domain simplifies the mathematics and provides a clear reference scale
Grid points (NX = 500): This provides high spatial resolution with \(dx \approx 0.002\), ensuring accurate capture of concentration gradients
Grid spacing calculation: \(dx = L/(NX+1)\) accounts for the boundary points being excluded from the interior grid

Temporal Parameters: The critical choice here is \(dt = dx^2/4\). This is not arbitrary but based on stability theory:

For explicit schemes, the stability condition is \(D \cdot dt/dx^2 \leq 0.5\)
Here, \(dt = dx^2/4\) ensures \(D \cdot dt/dx^2 = D/4 \approx 0.025\) (with \(D=0.1\)), well within the stable range
This conservative choice prioritizes numerical stability over computational speed

Physical Parameter: The diffusion coefficient \(D = 0.1\) provides a moderate diffusion rate that allows observable spreading over the simulation time \(T = 0.5\).

How to analyze stability

Understanding numerical stability is crucial for reliable simulations. The key insight is that when we discretize differential equations, we must ensure that small errors don’t grow exponentially and destroy our solution.

The Stability Number

The stability number \(r = D \cdot dt/dx^2\) is a dimensionless parameter that tells us whether our numerical scheme will behave well:

Physical meaning: \(r\) represents the ratio of how far information can diffuse in one time step (\(\sqrt{D \cdot dt}\)) compared to our grid spacing (\(dx\))
Mathematical meaning: It measures how much the solution can change between adjacent grid points in one time step

Interpreting Stability Values

Values \(r > 0.5\):

Can cause catastrophic instability in explicit schemes (solution “blows up”)
Even in stable implicit schemes like Crank-Nicolson, accuracy degrades significantly
Tip: If your simulation produces wildly oscillating or growing solutions, check if \(r\) is too large

Values \(r < 0.01\):

Very conservative - the solution will be accurate but computationally expensive
You’re taking much smaller time steps than necessary
Tip: If your simulation runs very slowly, you might be able to increase \(dt\) safely

Values \(0.01 \leq r \leq 0.5\):

The “sweet spot” for most diffusion problems
Good balance of computational efficiency and numerical accuracy
Tip: Aim for this range in your simulations

Why Crank-Nicolson Still Needs Stability Analysis

While the Crank-Nicolson scheme is theoretically “unconditionally stable” (it won’t blow up for any time step), accuracy still matters:

Large time steps can introduce significant numerical errors
The solution might be stable but completely wrong
Rule of thumb: Stability analysis helps ensure both stability AND accuracy

Practical Stability Analysis Steps

Calculate your stability number: \(r = D \cdot dt/dx^2\)
Check the range: Is \(0.01 \leq r \leq 0.5\)?
If \(r\) is too large: Decrease \(dt\) or increase \(dx\) (coarsen the grid)
If \(r\) is too small: You can safely increase \(dt\) to speed up computation
Always verify: Run a test with known analytical solution to validate your parameters

In our example, choosing \(dt = dx^2/4\) gives \(r \approx 0.025\), which ensures both stability and accuracy while keeping computational cost reasonable.

Initial Conditions

We create our initial Gaussian concentration profile:

Let’s visualize our initial condition:

Matrix Construction

We build the second-difference matrix using SciPy’s sparse matrix tools. Let’s start with the main matrix for our simulation, which implements Dirichlet boundary conditions (fixed concentration values of zero at both boundaries).

For comparison, here’s how to construct the matrix with Neumann boundary conditions (zero-flux at boundaries, meaning \(\frac{\partial c}{\partial x} = 0\) at both ends):

Time Evolution

Now we solve the system iteratively. Let’s break this down into manageable steps:

Step 1: Pre-compute the Matrices

Remember our equation \(\mathbf{A}\mathbf{C}^{n+1}=\mathbf{B}\mathbf{C}^{n}\)? Since the matrices A and B contain only our grid spacing, time step, and diffusion coefficient - all of which stay constant - we can calculate them once at the beginning rather than recalculating them thousands of times in our loop. This saves enormous amounts of computation time.

Step 2: Set Up Storage

As our simulation runs through time, we want to save snapshots of how the concentration profile looks at different moments. Think of this like taking photos of the spreading ink drop at regular intervals. We’ll store these “snapshots” in a list so we can analyze and visualize the evolution later.

Step 3: Main Time Stepping Loop

This is where the real physics happens! At each time step, we take our current concentration profile and use our diffusion equation to predict what it will look like one small time step later. We repeat this process thousands of times to watch the diffusion unfold. It’s like predicting the future in tiny increments - each prediction becomes the starting point for the next prediction.

The key insight is that we’re solving \(\mathbf{A}\mathbf{C}^{n+1}=\mathbf{B}\mathbf{C}^{n}\) at each step: we know \(\mathbf{C}^{n}\) (current concentration), so we calculate \(\mathbf{B}\mathbf{C}^{n}\) and then solve for the unknown \(\mathbf{C}^{n+1}\) (future concentration).

Step 4: Quick Validation Check

Even the best numerical methods can sometimes produce unphysical results due to rounding errors or inappropriate parameters. As good computational physicists, we should always check whether our results make sense from a physics perspective. Are there negative concentrations (impossible!)? Do particles disappear from the boundaries as expected? Is mass conserved (particles don’t just vanish)?

Visualization and Analysis

Evolution of Concentration Profiles

Let’s visualize how the concentration profile evolves over time:

Notice how the sharp initial peak gradually spreads out and flattens over time, while the total area under the curve remains constant (conservation of mass).

Space-Time Evolution

A particularly insightful visualization is the space-time plot, which shows the entire evolution as a 2D image:

Quantitative Analysis

Let’s analyze some key properties of our solution:

This analysis confirms that our numerical solution matches the theoretical prediction: the width of a diffusing distribution grows as \(\sqrt{2Dt}\).

Extensions and Explorations

Different Boundary Conditions

Try modifying the boundary conditions to see their effects:

Neumann boundaries (zero flux at boundaries):

# Modify the matrix M to implement zero-flux conditions
M[0, 0] = -1/dx**2    # First row
M[-1, -1] = -1/dx**2  # Last row

Periodic boundaries (what goes out one side comes in the other):

# Add wraparound connections
M[0, -1] = 1/dx**2
M[-1, 0] = 1/dx**2

Parameter Exploration

Investigate how different parameters affect the solution:

Diffusion coefficient D: Try values from 0.01 to 1.0
Initial width σ: Compare narrow vs. wide initial distributions
Grid resolution: See how dx affects accuracy
Time step: Explore the stability limit

Physical Applications

This same mathematical framework applies to:

Heat conduction: Replace concentration with temperature, D becomes thermal diffusivity
Quantum mechanics: The time-dependent Schrödinger equation has the same mathematical structure
Finance: Black-Scholes equation for option pricing
Biology: Population dynamics and chemotaxis

Where to Go From Here

The diffusion equation is your gateway to the vast world of partial differential equations. Building on what you’ve learned here, you could explore:

Two-dimensional diffusion: Extend to diffusion in a plane
Nonlinear diffusion: When D depends on concentration
Reaction-diffusion systems: Combine diffusion with chemical reactions
Wave equations: Replace first time derivative with second
Advanced numerical methods: Finite element methods, spectral methods

The techniques you’ve mastered - discretization, matrix methods, and iterative solution - form the foundation for computational physics and engineering. These tools will serve you well as you tackle increasingly complex physical phenomena.