pysdic.compute_central_finite_difference_coefficients#
- compute_central_finite_difference_coefficients(order, spacing=1.0, accuracy=2)[source]#
Compute the coefficients for the central finite difference approximation of a derivative.
The function returns the Stencil coefficients \(c_j\) used to approximate the n-th order derivative of a function using central finite differences.
\[\frac{d^n f}{dt^n} \approx \frac{1}{h^n} \sum_{j=-M}^{M} c_j f(t + j h)\]where \(h\) is the time step size, and \(m\) is the accuracy order and \(N\) is the number of stencil points required to achieve the desired accuracy (\(N = n + m - 1\)), such that the approximation error is of order \(O(h^m)\). Here, \(M = \frac{N}{2}\) if \(N\) is even, and \(M = \frac{N-1}{2}\) if \(N\) is odd.
\[\frac{1}{h^n} \sum_{j=-M}^{M} c_j f(t + j h) = \frac{d^n f}{dt^n} + O(h^m)\]Note
The coefficients are computed using Gauss elimination on the Vandermonde matrix constructed from the Taylor series expansion.
- Parameters:
order (Integral) – The order of the derivative to approximate (e.g., 1 for first derivative, 2 for second derivative).
spacing (Number, optional) – The time step size \(h\). Default is
1.0.accuracy (Integral, optional) – The desired accuracy order of the approximation. Must be a even integer. Default is
2.
- Returns:
The coefficients for the central finite difference approximation of the derivative in the form of a 1D array
[c_{-Mh}, ..., c_{-h}, c_{0}, c_{h}, ..., c_{Mh}].- Return type:
- Raises:
ValueError – If the order or accuracy are not positive integers. If the accuracy is not an even integer. If the spacing is not a strictly positive number.
Notes
By developing the function \(f(t + k h)\) in Taylor series around \(t\), we have:
\[f(t + j h) = \sum_{k=0}^{\infty} \frac{(j h)^k}{k!} \frac{d^k f}{dt^k}\]We can add the contributions of the stencil points weighted by the coefficients \(c_k\):
\[\frac{1}{h^n} \sum_{j=-M}^{M} c_j f(t + j h) = \frac{1}{h^n} \sum_{j=-M}^{M} c_j \sum_{k=0}^{\infty} \frac{(j h)^k}{k!} \frac{d^k f}{dt^k}\]\[\frac{1}{h^n} \sum_{j=-M}^{M} c_j f(t + j h) = \sum_{k=0}^{\infty} \frac{h^{k-n}}{k!} \left(\sum_{j=-M}^{M} c_j j^k\right) \frac{d^k f}{dt^k}\]Thus \(\forall k \in [0, N]\), we want to enforce the conditions:
\[\sum_{j=-M}^{M} c_j j^k = n! \delta_{k,n}\]where \(\delta_{k,n}\) is the Kronecker delta.
So we have a linear system of equations to solve for the coefficients \(c_j\). We can demonstrate that the matrix of the system is a Vandermonde matrix, which is invertible as long as the stencil points are distinct if we consider symmetry/antisymmetry of the coefficients depending on the order. And the output coefficients are divided by \(h^n\) to account for the time step size. The result is a finite difference approximation of the n-th order derivative with an error of order \(O(h^m)\).
See also
compute_forward_finite_difference_coefficients()To compute the forward finite difference coefficients.
compute_backward_finite_difference_coefficients()To compute the backward finite difference coefficients.
assemble_central_finite_difference_matrix()To assemble the central finite difference operator matrix.
Examples
>>> compute_central_finite_difference_coefficients(1, 1.0, 2) array([-0.5, 0. , 0.5])
>>> compute_central_finite_difference_coefficients(2, 1.0, 2) array([ 1., -2., 1.])
>>> compute_central_finite_difference_coefficients(1, 1.0, 4) array([ 1/12, -2/3, 0, 2/3, -1/12])
>>> compute_central_finite_difference_coefficients(2, 0.1, 4) array([ -100/12, 400/3, -500/2, 400/3, -100/12])