Finite-Difference Discretization of Diffusion

From Diffusion to a Grid Update

A tissue diffusion model describes concentration as changing continuously in space and time. A numerical simulation instead stores concentrations at selected locations and times. The finite-difference scheme specifies how those stored values are updated. Its spatial grid and time step are therefore part of the numerical model: a poor choice can produce oscillations or misleading results even when the underlying diffusion equation is appropriate.

Consider one-dimensional diffusion with constant diffusivity \(D>0\), no reaction term, and concentration \(c(x,t)\):

\frac{\partial c}{\partial t} = D\frac{\partial^2 c}{\partial x^2}.

Here \(D\) has units of length squared per time. For example, if \(x\) is measured in millimetres and \(t\) in minutes, \(D\) is in \(\mathrm{mm^2/min}\). The equation must be accompanied by an initial concentration profile and boundary conditions. In the examples below, the spatial domain is a uniform grid with specified boundary concentrations. We use \(x_i=i\Delta x\) and \(t_n=n\Delta t\), and write \(c_i^n\) for the numerical approximation to \(c(x_i,t_n)\).

Definition: A finite-difference discretization replaces spatial and time derivatives with expressions formed from values at grid points. The explicit forward-time, centered-space (FTCS) scheme uses known concentrations at time \(t_n\) to compute those at \(t_{n+1}\). The implicit backward-time, centered-space scheme uses unknown concentrations at \(t_{n+1}\) in its spatial-difference expression, so the new values must be solved for together.

The central second difference approximates the spatial curvature at an interior grid point:

\frac{\partial^2 c}{\partial x^2}(x_i,t_n) \approx \frac{c_{i-1}^n-2c_i^n+c_{i+1}^n}{(\Delta x)^2}.

This approximation follows by expanding \(c(x_i-\Delta x,t_n)\) and \(c(x_i+\Delta x,t_n)\) in Taylor series about \(x_i\). Adding the expansions cancels the first-derivative terms; after subtracting \(2c(x_i,t_n)\) and dividing by \((\Delta x)^2\), the leading term is \(c_{xx}(x_i,t_n)\), with an error of order \((\Delta x)^2\). The centered difference therefore approximates curvature to second order in space, provided the concentration profile is sufficiently smooth.

Explicit and Implicit Schemes

For the explicit scheme, approximate the time derivative by a forward difference and the spatial derivative by the centered difference at the known time \(t_n\). Define the dimensionless ratio

r=\frac{D\Delta t}{(\Delta x)^2}.

Substitution into the diffusion equation gives the update for each interior point:

c_i^{n+1} = c_i^n+r\left(c_{i-1}^n-2c_i^n+c_{i+1}^n\right) = (1-2r)c_i^n+rc_{i-1}^n+rc_{i+1}^n.

The boundary values at each time step are supplied by the chosen boundary conditions. The explicit formula is inexpensive because each new interior value is calculated directly from the previous time level. Its time approximation is first order, while its centered spatial approximation is second order.

In the implicit scheme, evaluate the centered second difference at the new time \(t_{n+1}\), and use a backward time difference. The resulting equations are

c_i^{n+1} - r\left(c_{i-1}^{n+1}-2c_i^{n+1}+c_{i+1}^{n+1}\right) = c_i^n,

or, equivalently,

-r c_{i-1}^{n+1}+(1+2r)c_i^{n+1}-r c_{i+1}^{n+1}=c_i^n.

For multiple interior points these equations form a linear system. Known boundary values contribute to the right-hand side when they are nonzero. Backward differencing is first order in time and the centered spatial difference remains second order in space. The need to solve a system at every time step is the main computational cost compared with the explicit update.

Worked Example: One Explicit Update

Suppose a hypothetical tissue model has \(D=0.0100\ \mathrm{mm^2/min}\), \(\Delta x=0.100\ \mathrm{mm}\), and \(\Delta t=0.200\ \mathrm{min}\). There are three interior points, \(i=1,2,3\), between boundary points \(i=0\) and \(i=4\). At time \(t_n\), let the boundary concentrations be \(c_0^n=c_4^n=0\ \mathrm{mg/L}\), and the interior concentrations be \(c_1^n=2.00\), \(c_2^n=1.00\), and \(c_3^n=0.500\ \mathrm{mg/L}\). Calculate one FTCS update.

First calculate the ratio:

r=\frac{(0.0100\ \mathrm{mm^2/min})(0.200\ \mathrm{min})}{(0.100\ \mathrm{mm})^2} =\frac{0.00200}{0.0100}=0.200.

Apply the update separately at each interior point:

\begin{aligned} c_1^{n+1} &=2.00+0.200(0-2(2.00)+1.00) =1.40\ \mathrm{mg/L},\\ c_2^{n+1} &=1.00+0.200(2.00-2(1.00)+0.500) =1.10\ \mathrm{mg/L},\\ c_3^{n+1} &=0.500+0.200(0-2(0.500)+1.00) =0.500\ \mathrm{mg/L}. \end{aligned}

For the last point, the curvature term is \(0-1.00+1.00=0\), so its concentration does not change in this step. All three new concentrations are in \(\mathrm{mg/L}\), and the boundary values remain as specified.

Worked Example: Solving an Implicit Step

Use the same hypothetical \(r=0.200\), but consider two interior points with zero boundary concentrations. Let \(c_1^n=1.00\ \mathrm{mg/L}\) and \(c_2^n=0\ \mathrm{mg/L}\). The implicit equations are

\begin{aligned} 1.4c_1^{n+1}-0.2c_2^{n+1}&=1.00,\\ -0.2c_1^{n+1}+1.4c_2^{n+1}&=0. \end{aligned}

The determinant of the coefficient matrix is \(1.4(1.4)-0.2(0.2)=1.92\), so it is nonzero and the system has a unique solution. Solving gives

c_1^{n+1}=\frac{1.4}{1.92}=0.7292\ \mathrm{mg/L}, \qquad c_2^{n+1}=\frac{0.2}{1.92}=0.1042\ \mathrm{mg/L}.

Both updated concentrations are nonnegative. The second point becomes positive because the implicit system couples the new concentration at that point to the new value at its neighbour.

Why the Explicit Time Step Is Restricted

A stability analysis asks whether small errors in a numerical solution grow as time steps are taken. For the explicit diffusion scheme on a uniform infinite or periodic grid, consider a Fourier mode \(e^{\mathrm{i}i\theta}\), where \(\theta=k\Delta x\) is the dimensionless spatial wave number. If its amplitude changes by a factor \(G\) in one step, substitution into the FTCS update gives

G=1+r(e^{-\mathrm{i}\theta}-2+e^{\mathrm{i}\theta}) =1-4r\sin^2\left(\frac{\theta}{2}\right).

For stability, the magnitude of the amplification factor must not exceed one for the modes under consideration. If the Fourier-symbol analysis considers all \(\theta\) in the continuous interval from \(0\) to \(2\pi\), then \(\sin^2(\theta/2)\) ranges from 0 to 1. For \(r\geq 0\), \(G\) is at most 1, and its smallest value is \(1-4r\). Requiring this value to be at least \(-1\) yields

1-4r\geq -1 \quad\Longrightarrow\quad r\leq\frac{1}{2} \quad\Longrightarrow\quad \Delta t\leq\frac{(\Delta x)^2}{2D}.

This is the standard grid-independent sufficient restriction for the one-dimensional FTCS scheme. It also ensures that the three weights \(1-2r,r,r\) in the update are nonnegative and sum to one, so each updated value is a weighted average of its neighbours and itself. This property helps prevent a new value from falling outside the range of those three old values.

There is an important finite-grid qualification. A finite periodic grid allows only discrete wave numbers. On such a grid, the exact modal condition is \(r\leq 1/(2s_{\max})\), where \(s_{\max}\) is the largest allowed value of \(\sin^2(\theta/2)\). This equals \(1/2\) when an allowed mode has \(\theta=\pi\), as on an even grid; on an odd grid it can be slightly larger. Thus the bound \(r\leq 1/2\) remains a sufficient condition for every such grid, but the continuous-wave-number argument should not be mistaken for an exact finite-grid threshold in every case. For ordinary use, the sufficient bound is simple and robust.

For the implicit scheme, the same Fourier-mode calculation gives

G_{\mathrm{implicit}} = \frac{1}{1+4r\sin^2(\theta/2)}.

For \(r\geq0\), this amplification factor is between 0 and 1 for every mode. Thus the linear backward-time scheme for constant-coefficient diffusion is unconditionally stable in this Fourier sense: stability does not impose an upper bound on \(\Delta t\). This does not mean that any time step gives an accurate solution. A large step can smooth the profile too much or miss important transient behaviour.

Worked Example: Calculating a Stable Explicit Time Step

Suppose a hypothetical model uses \(D=0.0200\ \mathrm{mm^2/min}\) and a spatial spacing of \(\Delta x=0.100\ \mathrm{mm}\). The standard explicit restriction gives

\Delta t\leq\frac{(0.100\ \mathrm{mm})^2}{2(0.0200\ \mathrm{mm^2/min)}} =\frac{0.0100\ \mathrm{mm^2}}{0.0400\ \mathrm{mm^2/min}} =0.250\ \mathrm{min}.

A step of \(0.200\ \mathrm{min}\) has \(r=0.400\) and satisfies the bound. A step of \(0.300\ \mathrm{min}\) has \(r=0.600\) and violates it. On an even periodic grid, the highest-frequency mode has \(\theta=\pi\), giving \(G=1-4r=1-2.4=-1.4\) for \(r=0.600\). Its magnitude exceeds one, so that mode alternates sign and grows from step to step. This is a numerical instability, not evidence that physical diffusion creates alternating concentrations.

Choosing and Checking a Scheme

The explicit and implicit schemes make different trade-offs. Explicit updates are straightforward and inexpensive per step, but the stability restriction can require many small steps. The implicit scheme requires a linear solve but avoids that diffusion-driven stability restriction. Both schemes still have discretization error, and neither removes the need to check whether the spatial and temporal resolution is adequate for the model’s purpose.

Feature	Explicit FTCS	Implicit backward-time
Update	Direct calculation from the old time level	Linear system for the new time level
Time accuracy	First order	First order
Diffusion stability	Requires the usual \(r\leq 1/2\) bound	Unconditionally stable for linear diffusion in Fourier analysis
Practical caution	Small \(\Delta x\) sharply restricts \(\Delta t\)	A stable large step may still be too inaccurate

Several checks help identify problems in an implementation:

Check units and the ratio. \(D\Delta t/(\Delta x)^2\) is dimensionless. A value computed with inconsistent length or time units has no valid stability interpretation.
Check the explicit bound before advancing. If \(r>1/2\), the standard sufficient condition is violated; oscillations or growing grid-scale patterns may result.
Inspect the update weights. When \(r>1/2\), the coefficient \(1-2r\) is negative. A computed concentration can then overshoot even if the old concentrations were nonnegative.
Check the boundary treatment. The update at the first and last interior points must use the specified boundary values or the appropriate boundary-difference rule. Incorrect indices can mimic instability or create artificial sources.
Do not equate stability with accuracy. A bounded solution can still be too smoothed or otherwise differ materially from the continuous solution. Stability is a necessary numerical safeguard, not a validation of model predictions.

A useful alternative to simply reducing \(\Delta t\) is to use an implicit update when the explicit stability limit is computationally burdensome. The choice depends on the desired temporal resolution, the cost of solving the system, and the need to resolve rapid changes. This tutorial treats constant \(D\) and a uniform one-dimensional grid; spatially varying coefficients, nonlinear fluxes, and more complex geometries require additional discretization decisions.

Common Mistake

Pitfall: “The implicit scheme is unconditionally stable, so I can choose a very large time step without affecting the simulation.” Stability only says that numerical errors do not grow without bound under the stated linear diffusion analysis. It does not guarantee that the simulated concentration profile is accurate. A time step that is too large can erase transient structure or skip biologically relevant changes. Choose time resolution based on the dynamics to be represented, and assess numerical adequacy rather than relying on stability alone.

Check Your Understanding

Answer each question, then compare with the solutions below.

For \(D=0.00500\ \mathrm{mm^2/min}\) and \(\Delta x=0.0500\ \mathrm{mm}\), calculate the standard maximum explicit time step.
An explicit calculation uses \(r=0.300\). Write its three update weights and state whether they are nonnegative and sum to one.
For an explicit periodic-grid mode with \(\sin^2(\theta/2)=0.25\) and \(r=0.600\), calculate its amplification factor. Does this particular mode grow in magnitude?
Why does implicit stability not establish that a chosen time step is accurate?

Solutions. (1) \(\Delta t\leq(\Delta x)^2/(2D)=(0.0500\ \mathrm{mm})^2/[2(0.00500\ \mathrm{mm^2/min})]=0.250\ \mathrm{min}\). (2) The weights are \(1-2r=0.400\), \(r=0.300\), and \(r=0.300\). They are nonnegative and sum to \(1.00\). (3) \(G=1-4(0.600)(0.25)=0.400\). This mode does not grow in magnitude, although \(r=0.600\) exceeds the grid-independent sufficient bound; stability must be checked for all allowed modes, not inferred from one mode. (4) Stability limits growth of numerical errors but does not ensure that the numerical solution resolves the physical time course. A large stable step may still smooth or skip relevant changes.