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 ( |
side |
Integer that determines the type of finite-difference scheme:
|
acc.order |
Order of accuracy: defines how the approximation error scales
with the step size |
stencil |
Optional custom vector of points for function evaluation.
Must include at least |
zero.action |
Character string specifying how to handle near-zero weights:
|
zero.tol |
Non-negative scalar defining the threshold: weights below
|
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)