fdCoef: Finite-difference coefficients for arbitrary grids

In pnd: Parallel Numerical Derivatives, Gradients, Jacobians, and Hessians of Arbitrary Accuracy Order

Finite-difference coefficients for arbitrary grids

Description

This function computes the coefficients for a numerical approximation to derivatives of any specified order. It provides the minimally sufficient stencil for the chosen derivative order and desired accuracy order. It can also use any user-supplied stencil (uniform or non-uniform). For that stencil \{b_i\}_{i=1}^n, it computes the optimal weights \{w_i\} that yield the numerical approximation of the derivative:

\frac{d^m f}{dx^m} \approx h^{-m} \sum_{i=1}^n w_i f(x + b_i\cdot h)

Usage

fdCoef(
  deriv.order = 1L,
  side = c(0L, 1L, -1L),
  acc.order = 2L,
  stencil = NULL,
  zero.action = c("drop", "round", "none"),
  zero.tol = NULL
)

Arguments

deriv.order	Order of the derivative (m in \frac{d^m f}{dx^m}) for which a numerical approximation is needed.
side	Integer that determines the type of finite-difference scheme: 0 for central (AKA symmetrical or two-sided; the default), 1 for forward, and -1 for backward. Using 2 (for 'two-sided') triggers a warning and is treated as 0. with a warning. Unless the function is computationally prohibitively, central differences are strongly recommended for their accuracy.
acc.order	Order of accuracy: defines how the approximation error scales with the step size h, specifically O(h^{a+1}), where a is the accuracy order and depends on the higher-order derivatives of the function.
stencil	Optional custom vector of points for function evaluation. Must include at least m+1 points for the m-th order derivative.
zero.action	Character string specifying how to handle near-zero weights: "drop" to omit small (less in absolute value than zero.tol times the median weight) weights and corresponding stencil points, "round" to round small weights to zero, and "none" to leave all weights as calculated. E.g. the stencil for f'(x) is (-1, 0, 1) with weights (-0.5, 0, 0.5); using "drop" eliminates the zero weight, and the redundant f(x) is not computed.
zero.tol	Non-negative scalar defining the threshold: weights below zero.tol times the median weight are considered near-zero.

Details

This function relies on the approach of approximating numerical derivarives by weghted sums of function values described in \insertCitefornberg1988generationpnd. It reproduces all tables from this paper exactly; see the example below to create Table 1.

The finite-difference coefficients for any given stencil are given as a solution of a linear system. The capabilities of this function are similar to those of \insertCitetaylor2016finitepnd, but instead of matrix inversion, the \insertCitebjorck1970solutionpnd algorithm is used because the left-hand-side matrix is a Vandermonde matrix, and its inverse may be very inaccurate, especially for long one-sided stencils.

The weights computed for the stencil via this algorithm are very reliable; numerical simulations in \insertCitehigham1987errorpnd show that the relative error is low even for ill-conditioned systems. \insertCitekostyrka2025whatpnd computes the exact relative error of the weights on the stencils returned by this function; the zero tolerance is based on these calculations.

Value

A list containing the stencil used and the corresponding weights for each point.

References

\insertAllCited

Examples

fdCoef()  # Simple two-sided derivative
fdCoef(2) # Simple two-sided second derivative
fdCoef(acc.order = 4)$weights * 12  # Should be (1, -8, 8, -1)

# Using an custom stencil for the first derivative: x-2h and x+h
fdCoef(stencil = c(-2, 1), acc.order = 1)

# Reproducing Table 1 from Fornberg (1988) (cited above)
pad9 <- function(x) {l <- length(x); c(a <- rep(0, (9-l)/2), x, a)}
f <- function(d, a) pad9(fdCoef(deriv.order = d, acc.order = a,
                                zero.action = "round")$weights)
t11 <- t(sapply((1:4)*2, function(a) f(d = 1, a)))
t12 <- t(sapply((1:4)*2, function(a) f(d = 2, a)))
t13 <- t(sapply((1:3)*2, function(a) f(d = 3, a)))
t14 <- t(sapply((1:3)*2, function(a) f(d = 4, a)))
t11 <- cbind(t11[, 1:4], 0, t11[, 5:8])
t13 <- cbind(t13[, 1:4], 0, t13[, 5:8])
t1 <- data.frame(OrdDer = rep(1:4, times = c(4, 4, 3, 3)),
                 OrdAcc = c((1:4)*2, (1:4)*2, (1:3)*2, (1:3)*2),
                 rbind(t11, t12, t13, t14))
colnames(t1)[3:11] <- as.character(-4:4)
print(t1, digits = 4)

pnd documentation built on Sept. 9, 2025, 5:44 p.m.