Finite Difference Method

What the finite difference method does

It converts a differential equation or derivative calculation into a system of numerical equations that a computer can evaluate.

This page develops forward, backward, central, and second-derivative formulas from Taylor expansions, then connects them to boundary-value and time-dependent problems. It also explains consistency, stability, convergence, grid refinement, round-off, and health-economic uses such as numerical sensitivity and continuous-time decision models.

A grid replaces a continuous domain

Finite differences begin by selecting discrete points within the continuous independent-variable domain. On a uniform one-dimensional grid, neighbouring points are separated by constant spacing $h$.

$$x_i=x_0+ih, \qquad i=0,1,\ldots,N$$

The numerical solution stores an approximation $u_i\approx u(x_i)$ at each point. Smaller spacing can improve resolution, but it increases computational work and does not guarantee accuracy when the numerical scheme is unstable or dominated by round-off.

Forward and backward first differences

One-sided formulas use the current point and a neighbour on only one side. They are useful at boundaries or when a time-stepping method is defined in a particular direction.

The forward difference is:

$$f'(x) \approx \frac{f(x+h)-f(x)}{h}$$

The backward difference is:

$$f'(x) \approx \frac{f(x)-f(x-h)}{h}$$

For a sufficiently smooth function, both formulas have truncation error of order $O(h)$. Halving $h$ should therefore reduce the leading truncation error by approximately a factor of two before floating-point effects become important.

The central first difference

A central difference uses values symmetrically located around the evaluation point. Symmetry cancels the leading first-order truncation term and usually provides greater accuracy than a one-sided difference at an interior point.

$$f'(x) \approx \frac{f(x+h)-f(x-h)}{2h}$$

For a sufficiently smooth function, its leading truncation error is $O(h^2)$. Halving $h$ should reduce that leading error by approximately a factor of four within the asymptotic convergence region.

The central second difference

Second derivatives describe curvature and appear in diffusion, smoothing, optimisation, and many boundary-value problems. A symmetric three-point stencil gives the standard central approximation.

$$f''(x) \approx \frac{f(x+h)-2f(x)+f(x-h)}{h^2}$$

This formula also has truncation error $O(h^2)$ when the required higher derivatives exist. The division by $h^2$ makes it especially sensitive to cancellation and noisy input values when $h$ is extremely small.

Taylor series explain the approximation orders

Finite difference formulas can be derived by expanding neighbouring values around $x$. The retained terms produce the derivative approximation, while the omitted terms determine the local truncation error.

$$f(x+h)=f(x)+hf'(x)+\frac{h^2}{2}f''(x) +\frac{h^3}{6}f'''(x)+O(h^4)$$ $$f(x-h)=f(x)-hf'(x)+\frac{h^2}{2}f''(x) -\frac{h^3}{6}f'''(x)+O(h^4)$$

Subtracting the second expansion from the first and dividing by $2h$ produces the central first difference. Adding the expansions, subtracting $2f(x)$, and dividing by $h^2$ produces the central second difference.

A worked derivative example

Let $f(x)=x^2$, estimate derivatives at $x=2$, and use $h=0.1$. The exact values are $f'(2)=4$ and $f''(2)=2$.

The forward estimate is:

$$\frac{f(2.1)-f(2)}{0.1} =\frac{4.41-4}{0.1} =4.1$$

The backward estimate is:

$$\frac{f(2)-f(1.9)}{0.1} =\frac{4-3.61}{0.1} =3.9$$

The central estimate is:

$$\frac{f(2.1)-f(1.9)}{0.2} =\frac{4.41-3.61}{0.2} =4.0$$

The central second-difference estimate is:

$$\frac{4.41-2(4)+3.61}{0.1^2} =\frac{0.02}{0.01} =2.0$$

The central formulas are exact here because all omitted derivatives of the quadratic are zero. Exactness for this example does not mean they are exact for general nonlinear functions.

A stencil shows which points are used

A finite difference stencil describes the neighbouring grid values and weights used to approximate a derivative. It makes data dependencies, boundary requirements, and computational structure visible.

For the central second derivative, the weights are:

$$\frac{1}{h^2} \begin{bmatrix} 1&-2&1 \end{bmatrix}$$

applied to $[f_{i-1},f_i,f_{i+1}]^{\mathsf T}$. Wider stencils can achieve higher formal accuracy but need more neighbouring values, complicate boundaries, and may behave poorly around discontinuities or under-resolved features.

Turning a boundary-value problem into linear equations

Consider the second-order boundary-value problem $-u''(x)=g(x)$ on $0<x<L$, with specified endpoint values. Replacing the second derivative at each interior grid point produces one algebraic equation per unknown interior value.

$$-\frac{u_{i-1}-2u_i+u_{i+1}}{h^2} =g_i$$

The interior equations form a sparse tridiagonal linear system. Boundary values enter the right-hand side for equations whose stencil touches an endpoint, so an incorrect boundary implementation can contaminate the entire numerical solution.

Boundary conditions change the discrete equations

Dirichlet conditions specify the solution value, Neumann conditions specify a derivative or flux, and Robin conditions combine values and derivatives. Each type requires a discrete equation consistent with the chosen accuracy and physical meaning.

A first-order one-sided Neumann condition at $x_0$ can be written as:

$$u'(x_0) \approx \frac{u_1-u_0}{h} =q$$

A higher-order boundary stencil or ghost-point construction may be needed to preserve the interior scheme's global order. Boundary formulas should be derived and tested rather than improvised from visual symmetry.

Time-dependent problems need a space-time grid

A partial differential equation with time and space uses separate increments $\Delta t$ and $\Delta x$. The discrete value $u_i^n$ approximates $u(x_i,t_n)$.

$$x_i=x_0+i\Delta x, \qquad t_n=t_0+n\Delta t$$

The numerical scheme specifies how values at one or more time levels and neighbouring spatial points determine the next solution. Accuracy and stability depend on the combined choice of temporal and spatial approximations.

An explicit scheme for the diffusion equation

For the one-dimensional diffusion equation $u_t=\alpha u_{xx}$, forward differencing in time and central differencing in space produce an explicit update. Every new interior value is calculated directly from values at the previous time level.

$$\frac{u_i^{n+1}-u_i^n}{\Delta t} =\alpha \frac{u_{i+1}^n-2u_i^n+u_{i-1}^n}{(\Delta x)^2}$$

Defining $r=\alpha\Delta t/(\Delta x)^2$ gives:

$$u_i^{n+1} =u_i^n+r \left(u_{i+1}^n-2u_i^n+u_{i-1}^n\right)$$

The method is first-order accurate in time and second-order accurate in space under smoothness assumptions. Its simplicity comes with a time-step stability restriction.

A worked diffusion update

Suppose three neighbouring values at time level $n$ are $u_{i-1}^n=0$, $u_i^n=1$, and $u_{i+1}^n=0$, with $r=0.25$. The update represents diffusion away from the central peak.

$$u_i^{n+1} =1+0.25(0-2(1)+0) =0.5$$

If adjacent grid points initially have zero with the central value as their only non-zero neighbour, their next values become 0.25 under the same stencil. The total is conserved locally in this idealised interior pattern, but global conservation also depends on boundary conditions and the complete grid.

Stability restricts the explicit diffusion step

For the standard one-dimensional explicit diffusion scheme on a uniform grid, von Neumann analysis gives a stability condition. The restriction links the time step to the square of the spatial step and becomes more demanding as the grid is refined. With non-negative diffusivity, the usual requirement is:

$$r=\frac{\alpha\Delta t}{(\Delta x)^2} \leq\frac{1}{2}$$

Refining the spatial grid by halving $\Delta x$ therefore requires reducing $\Delta t$ by approximately a factor of four to preserve the same $r$. Violating the condition can produce growing oscillations and numerical blow-up even when the differential equation's true solution is smooth and stable.

Implicit schemes trade direct updates for linear solves

Backward Euler evaluates the spatial difference at the new time level. The resulting method requires solving a linear system at each step but is unconditionally stable for the standard linear diffusion problem in the relevant norm.

$$\frac{u_i^{n+1}-u_i^n}{\Delta t} =\alpha \frac{u_{i+1}^{n+1}-2u_i^{n+1}+u_{i-1}^{n+1}} {(\Delta x)^2}$$

Unconditional stability does not mean unconditional accuracy. A very large time step can suppress important dynamics and introduce substantial temporal error even though the solution does not diverge.

Crank–Nicolson averages two time levels

Crank–Nicolson averages the centred spatial operator at time levels $n$ and $n+1$. For the linear diffusion equation, it is second-order accurate in time and requires a linear solve at each step.

$$\frac{u_i^{n+1}-u_i^n}{\Delta t} =\frac{\alpha}{2} \left( \delta_{xx}u_i^{n+1} +\delta_{xx}u_i^n \right)$$

Although stable for the standard problem, Crank–Nicolson can display non-physical oscillations when time steps are large relative to sharp initial gradients. Stability, monotonicity, and accuracy are different properties and should be checked separately.

Consistency, stability, and convergence

Consistency means the discrete equations approach the differential equation as the grid is refined. Stability means numerical perturbations do not grow without appropriate bound, while convergence means the numerical solution approaches the exact solution.

For suitable well-posed linear initial-value problems, consistency plus stability leads to convergence. A method with high local accuracy can still fail if it is unstable, and a stable method can converge too slowly to be useful at a practical grid size.

Local and global truncation error

Local truncation error is the residual obtained when the exact smooth solution is inserted into one discrete equation. Global error is the accumulated difference between the numerical and exact solutions across all steps and grid points.

The two orders need not be the same because local errors propagate through the numerical recurrence. Reporting a stencil's formal local order without a grid-convergence test can overstate the accuracy of the completed model.

Grid refinement reveals the observed order

If a method has error $E(h)\approx Ch^p$, reducing the grid spacing by a constant factor should reduce error predictably. Solutions on successively finer grids can estimate the observed convergence order without knowing the exact solution.

For refinement ratio 2 and comparable scalar results $Q_h$, $Q_{h/2}$, and $Q_{h/4}$:

$$p_{obs} \approx \log_2 \left| \frac{Q_h-Q_{h/2}} {Q_{h/2}-Q_{h/4}} \right|$$

The estimate is meaningful only when the solutions are in the asymptotic convergence region and the denominator is not dominated by round-off or cancellation. Irregular grids, nonlinear solvers, changing boundaries, and discontinuities require more careful comparisons.

Richardson extrapolation can improve an estimate

When the leading error behaves as $Ch^p$, two grid results can be combined to cancel that term. The extrapolated estimate is often more accurate than either original result.

$$Q_{ext} =Q_{h/2} +\frac{Q_{h/2}-Q_h}{2^p-1}$$

Extrapolation relies on a stable, known-order regime. It should not be used to disguise non-convergence, inconsistent grids, or a solution that changes qualitatively with refinement.

Round-off competes with truncation error

Decreasing $h$ reduces truncation error initially, but subtraction of nearly equal floating-point values increases cancellation error. Numerical differentiation therefore has an intermediate step size that is more accurate than either a very large or extremely small step.

For a forward difference, a simplified total-error model is:

$$E(h) \approx C_1h+C_2\frac{\varepsilon_{mach}}{h}$$

For a central difference, the leading truncation term is typically $C_1h^2$ while the round-off term still grows roughly as $\varepsilon_{mach}/h$. Step-size testing should look for a stable plateau rather than assuming that the smallest representable $h$ is best.

Non-uniform grids need different coefficients

Uniform-grid formulas cannot be copied unchanged onto irregularly spaced points. Coefficients must account for the left and right distances, and the approximation order depends on the resulting stencil.

Non-uniform grids can concentrate resolution where the solution changes rapidly and reduce computation elsewhere. Abrupt changes in grid spacing can degrade conditioning and accuracy, so grid design and coefficient generation should be validated together.

Advection requires direction-aware discretisation

For transport-dominated equations, centred spatial differences can generate oscillations because information travels in a preferred direction. Upwind schemes use values from the direction of incoming characteristics and can improve stability at the cost of numerical diffusion.

Higher-order upwind, flux-limited, or finite-volume methods may be needed when conservation and sharp fronts matter. Selecting a stencil solely by formal order ignores the equation's propagation structure.

Finite differences for model sensitivities

Health-economic models often use finite differences to estimate how an outcome changes with a parameter. These numerical derivatives support local sensitivity analysis, calibration, optimisation, uncertainty propagation, and validation of analytic gradients.

For model output $M(\boldsymbol{\theta})$, a central sensitivity to parameter $\theta_j$ is:

$$\frac{\partial M}{\partial\theta_j} \approx \frac{ M(\boldsymbol{\theta}+h\mathbf{e}_j) -M(\boldsymbol{\theta}-h\mathbf{e}_j) }{2h}$$

The perturbation should respect parameter bounds and any dependent parameters. A derivative near a hard threshold or discrete model branch may be unstable or not exist, in which case scenario analysis can be more informative than a local slope.

Health-economic and decision-model applications

Finite differences appear in health economics whenever continuous change, optimisation, or unavailable analytic derivatives must be represented numerically. Their usefulness depends on whether the grid and differential equation reflect the decision process rather than merely providing a convenient computation.

• Continuous-time decision models can discretise differential equations governing disease progression, treatment, or resource stocks.
• Optimal-control and dynamic-programming models can solve Hamilton–Jacobi–Bellman or related equations on state and time grids.
• Pharmacoeconomic models can approximate time-varying concentrations, hazards, or treatment effects when a differential-equation structure is justified.
• Calibration and optimisation can use finite-difference gradients when analytic or automatic derivatives are unavailable.
• Error propagation can use finite-difference sensitivities to build a local gradient for covariance-based uncertainty calculations.

Many ordinary cohort models use direct matrix recursion rather than finite-difference PDE methods. The method should be named only when derivatives or differential equations are actually discretised, not for every calculation performed at time steps.

Implementing derivative formulas in Excel

Excel can transparently evaluate finite differences for a function or model output. The worksheet should show the base parameter, perturbation $h$, plus and minus inputs, corresponding outputs, derivative estimate, and repeated estimates at alternative step sizes.

If $x$ is in B2, $h$ is in B3, and the formula for $f(x)$ is evaluated in dedicated cells, the central derivative is =(f_plus-f_minus)/(2*$B$3). The second derivative is =(f_plus-2*f_base+f_minus)/($B$3^2); named cells are preferable to repeating a long model formula.

For a full model, plus and minus scenarios must change only the intended parameter and must recalculate all dependent outputs. Volatile formulas, rounded displayed values, and manual copy-paste scenarios can invalidate a derivative even when the final formula looks correct.

Implementing a time-stepping grid in Excel

A small explicit model can place space across columns and time down rows. Each interior cell links only to the three neighbouring values in the preceding row and to a clearly labelled stability ratio $r$.

The update formula has the structure =current+r*(right-2*current+left). Boundary columns need separate formulas reflecting the declared boundary conditions, and a visible check should flag r>0.5 for the standard one-dimensional explicit diffusion scheme.

Excel is suitable for a teaching model and independent checks but becomes fragile for large, multidimensional, stiff, or nonlinear systems. Production work should use tested numerical libraries, sparse solvers, controlled tolerances, and automated convergence diagnostics.

Validation and numerical safeguards

Validation should test the discretisation, solver, boundaries, conservation properties, and decision outputs. Agreement on one grid can be accidental and does not demonstrate convergence.

• Analytic-solution checks: Problems with known solutions reproduce expected values and convergence orders.
• Polynomial checks: Difference formulas exactly differentiate polynomials up to the degree implied by their construction.
• Boundary checks: Every boundary condition is satisfied to the method's stated order and uses correct units and signs.
• Stability checks: Time steps and grid ratios meet applicable scheme-specific conditions.
• Refinement checks: Successively finer grids approach a stable result at the expected rate.
• Conservation checks: Mass, probability, population, or other conserved quantities remain balanced when the continuous model requires conservation.
• Positivity checks: States that must remain non-negative do not become negative through an unsuitable stencil or time step.
• Solver checks: Residuals, iteration tolerances, and sparse-system solutions meet prespecified thresholds.
• Replication checks: A hand calculation and an independent implementation reproduce selected grid updates and final outputs.

Decision validation should also show that incremental costs, outcomes, or recommended actions remain stable at a finer grid. A numerically converged solution can still be substantively wrong if the differential equation, initial conditions, or parameters are inappropriate.

Common errors and misunderstandings

Finite differences are simple to write but easy to misuse. Most serious failures involve a stencil applied outside its assumptions, an untested grid, or confusion between stability and accuracy.

• Making $h$ extremely small can increase round-off and cancellation instead of improving the derivative.
• Using uniform-grid coefficients on uneven points produces the wrong derivative approximation.
• Ignoring boundary accuracy can lower the global order even when the interior stencil is higher order.
• Violating a stability restriction can create oscillation or divergence unrelated to the true solution.
• Treating implicit as automatically accurate confuses freedom from a stability limit with adequate time resolution.
• Reporting formal order without refinement does not show that the implemented model reaches that order.
• Using central differences across a discontinuity can average incompatible regimes and generate non-physical oscillation.
• Calling every time-step model finite difference obscures whether an actual derivative or differential equation was discretised.

The finite difference method approximates a declared mathematical model. It does not correct biased parameters, an inappropriate causal structure, or an omitted clinical process.

A practical workflow

A defensible finite difference analysis starts with the continuous equation and ends with demonstrated convergence of the decision-relevant result. The workflow keeps the grid, stencil, boundaries, stability, and numerical error open to audit.

1. Define the continuous problem: State the differential equation, domain, units, initial conditions, boundary conditions, and desired output.
2. Choose the grid: Set spatial and temporal coordinates, spacing, refinement strategy, and treatment of irregular regions.
3. Derive the stencils: Use Taylor expansion or interpolation to obtain every interior and boundary formula and its formal order.
4. Select the scheme: Match explicit, implicit, centred, upwind, or other choices to the equation and propagation behaviour.
5. Check stability: Derive or verify all step-size and coefficient restrictions before running the model.
6. Implement transparently: Preserve sparse structure, units, indexing, signs, and independently testable intermediate values.
7. Verify calculations: Reproduce simple derivatives, one-step updates, and known analytic solutions.
8. Demonstrate convergence: Refine the grid, estimate observed order, and distinguish truncation from round-off or solver error.
9. Test decision robustness: Confirm that relevant costs, outcomes, and decisions are stable across justified numerical settings.
10. Report completely: Document equations, grids, stencils, boundaries, solvers, tolerances, errors, checks, and limitations.

What should be reported

A reproducible report states the differential equation, dependent and independent variables, units, domain, initial and boundary conditions, grid structure, finite difference stencils, formal accuracy, and time-stepping scheme. It also reports stability restrictions, solver and stopping tolerances, software precision, grid-refinement results, observed convergence order, conservation or positivity checks, and comparisons with analytic or independent solutions.

For health-economic applications, the report should identify which model quantities are generated by the finite difference solution and how numerical error affects incremental costs, effects, net benefit, or the optimal decision. Numerical settings should be included in sensitivity analysis whenever plausible alternatives can change a decision-relevant result.