fdCoef {pnd}R Documentation

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 (Fornberg 1988). 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 (Taylor 2016), but instead of matrix inversion, the (Björck and Pereyra 1970) 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 (Higham 1987) show that the relative error is low even for ill-conditioned systems. (Kostyrka 2025) 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

Björck Å, Pereyra V (1970). “Solution of Vandermonde systems of equations.” Mathematics of computation, 24(112), 893–903.

Fornberg B (1988). “Generation of Finite Difference Formulas on Arbitrarily Spaced Grids.” Mathematics of Computation, 51(184), 699–706. doi:10.1090/S0025-5718-1988-0935077-0.

Higham NJ (1987). “Error analysis of the Björck-Pereyra algorithms for solving Vandermonde systems.” Numerische Mathematik, 50(5), 613–632.

Kostyrka AV (2025). “What are you doing, step size: fast computation of accurate numerical derivatives with finite precision.” Working paper.

Taylor CR (2016). “Finite Difference Coefficients Calculator.” https://web.media.mit.edu/~crtaylor/calculator.html.

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)

[Package pnd version 0.1.0 Index]