<img src="docs/assets/findiff_logo.png" width="100px"> findiff

Coverage

A Python package for finite difference numerical derivatives and partial differential equations in any number of dimensions.

Main Features

Differentiate arrays of any number of dimensions along any axis with any desired accuracy order
Accurate treatment of grid boundary
Can handle uniform and non-uniform grids
Can handle arbitrary linear combinations of derivatives with constant and variable coefficients
Fully vectorized for speed
GPU / JAX / CuPy support — pass JAX or CuPy arrays directly, combine with jax.jit for acceleration
Standard operators from vector calculus: gradient, divergence, curl, Laplacian
Matrix representations of arbitrary linear differential operators
Solve partial differential equations with Dirichlet, Neumann or Robin boundary conditions
Solve eigenvalue problems (e.g. Schrodinger equation, vibration modes)
Direct and iterative solvers with preconditioner support
Calculate raw finite difference coefficients for any derivative and accuracy order
Generate differential operators for arbitrary stencils
Symbolic representation of finite difference schemes
Estimate truncation error by comparing accuracy orders
Solve time-dependent PDEs via Method of Lines (Forward Euler, RK4, Backward Euler, Crank-Nicolson)
New in version 0.11: More comfortable API (keeping the old API available)
New in version 0.12: Periodic boundary conditions for differential operators and PDEs.
New in version 0.13: Compact (implicit) finite differences with spectral-like resolution.
New in version 0.14: Error estimation via accuracy order comparison.
New in version 0.15: Time-dependent PDE solving via Method of Lines. GPU / JAX / CuPy backend support for operator application. (to be released)

Installation

pip install --upgrade findiff

For GPU / JAX support, install JAX separately (findiff detects it automatically):

pip install jax          # CPU-only
pip install jax[cuda12]  # NVIDIA GPU

For CuPy support:

pip install cupy-cuda12x

Documentation and Examples

You can find the documentation of the code including examples of application at https://maroba.github.io/findiff/.

Taking Derivatives

findiff allows to easily define derivative operators that you can apply to numpy arrays of any dimension. JAX and CuPy arrays work too — see GPU / JAX Support below.

Consider the simple 1D case of a equidistant grid with a first derivative $\displaystyle \frac{\partial}{\partial x}$ along the only axis (0):

import numpy as np
from findiff import Diff

# define the grid:
x = np.linspace(0, 1, 100)

# the array to differentiate:
f = np.sin(x)  # as an example

# Define the derivative:
d_dx = Diff(0, x[1] - x[0])

# Apply it:
df_dx = d_dx(f)

Similarly, you can define partial derivatives along other axes, for example, if $z$ is the 2-axis, we can write $\frac{\partial}{\partial z}$ as:

Diff(2, dz)

Diff always creates a first derivative. For higher derivatives, you simply exponentiate them, for example for $\frac{\partial^2}{\partial_x^2}$

d2_dx2 = Diff(0, dx)**2

and apply it as before.

You can also define more general differential operators intuitively, like

$$ 2x \frac{\partial^3}{\partial x^2 \partial z} + 3 \sin(y)z^2 \frac{\partial^3}{\partial x \partial y^2} $$

which can be written as

# define the operator
diff_op = 2 * X * Diff(0)**2 * Diff(2) + 3 * sin(Y) * Z**2 * Diff(0) * Diff(1)**2

# set the grid you use (equidistant here)
diff_op.set_grid({0: dx, 1: dy, 2: dz})

# apply the operator
result = diff_op(f)

where X, Y, Z are numpy arrays with meshed grid points. Here you see that you can also define your grid lazily.

Of course, standard operators from vector calculus like gradient, divergence and curl are also available as shortcuts.

If one or more axis of your grid are periodic, you can specify that when defining the derivative or later when setting the grid. For example:

d_dx = Diff(0, dx, periodic=True)

# or later
d_dx = Diff(0)
d_dx.set_grid({0: {"h": dx, "periodic": True}})

More examples can be found in the documentation and in this blog.

GPU / JAX Support

All operators (Diff, Gradient, Divergence, Curl, Laplacian) accept JAX and CuPy arrays directly. Just pass the array — findiff detects the backend automatically:

import jax
import jax.numpy as jnp
from findiff import Diff, Laplacian

jax.config.update("jax_enable_x64", True)

x = jnp.linspace(0, 2 * jnp.pi, 1000)
dx = x[1] - x[0]
f = jnp.sin(x)

d_dx = Diff(0, dx)
df_dx = d_dx(f)  # returns a JAX array

For best performance, wrap operators with jax.jit to fuse all slice operations into a single compiled kernel:

d_dx_jit = jax.jit(d_dx)
df_dx = d_dx_jit(f)  # first call compiles, subsequent calls are fast

# Works for any operator, including higher-order and vector calculus:
lap = Laplacian(h=[dx, dy, dz])
lap_jit = jax.jit(lap)
result = lap_jit(f_3d)

Note: The matrix-based path (.matrix(), PDE solving, eigenvalue problems) remains NumPy/SciPy only. GPU support applies to the operator application path (operator(array)).

Accuracy Control

When constructing an instance of Diff, you can request the desired accuracy order by setting the keyword argument acc. For example:

d_dx = Diff(0, dx, acc=4)
df_dx = d_dx(f)

Alternatively, you can also split operator definition and configuration:

d_dx = Diff(0, dx)
d_dx.set_accuracy(4)
df_dx = d_dx(f)

which comes in handy if you have a complicated expression of differential operators, because then you can specify it on the whole expression and it will be passed down to all basic operators.

If not specified, second order accuracy will be taken by default.

Compact Finite Differences

Standard finite differences only use function values to approximate a derivative. Compact (or implicit) finite differences also couple derivative values at neighboring points, which gives you spectral-like resolution from small stencils. The tradeoff is that applying the operator requires solving a banded linear system — but for tridiagonal systems that's $O(N)$ and very fast.

You can define a compact scheme explicitly by specifying the left-hand side coefficients and right-hand side offsets:

from findiff import Diff, CompactScheme

scheme = CompactScheme(
    deriv=1,
    left={-1: 1/3, 0: 1, 1: 1/3},     # tridiagonal LHS
    right=[-3, -2, -1, 0, 1, 2, 3],    # RHS offsets (coefficients computed automatically)
)

d_dx = Diff(0, dx, scheme=scheme, periodic=True)
df_dx = d_dx(np.sin(x))  # 6th-order accurate first derivative

Or let findiff pick a scheme for you:

d_dx = Diff(0, dx, compact=3, acc=6, periodic=True)

Here compact=3 sets the LHS bandwidth (must be odd), and findiff widens the RHS stencil until the requested accuracy is reached. Higher derivatives work the same as usual:

d2_dx2 = d_dx ** 2

Non-periodic grids are handled automatically — findiff uses one-sided compact stencils near the boundaries. The matrix() method is also supported. For more details, see the compact finite differences documentation.

Finite Difference Coefficients

Sometimes you may want to have the raw finite difference coefficients. These can be obtained for any derivative and accuracy order using findiff.coefficients(deriv, acc). For instance,

import findiff
coefs = findiff.coefficients(deriv=3, acc=4, symbolic=True)

gives

{'backward': {'coefficients': [15/8, -13, 307/8, -62, 461/8, -29, 49/8],
              'offsets': [-6, -5, -4, -3, -2, -1, 0]},
 'center': {'coefficients': [1/8, -1, 13/8, 0, -13/8, 1, -1/8],
            'offsets': [-3, -2, -1, 0, 1, 2, 3]},
 'forward': {'coefficients': [-49/8, 29, -461/8, 62, -307/8, 13, -15/8],
             'offsets': [0, 1, 2, 3, 4, 5, 6]}}

If you want to specify the detailed offsets instead of the accuracy order, you can do this by setting the offset keyword argument:

import findiff
coefs = findiff.coefficients(deriv=2, offsets=[-2, 1, 0, 2, 3, 4, 7], symbolic=True)

The resulting accuracy order is computed and part of the output:

{'coefficients': [187/1620, -122/27, 9/7, 103/20, -13/5, 31/54, -19/2835], 
 'offsets': [-2, 1, 0, 2, 3, 4, 7], 
 'accuracy': 5}

Matrix Representation

For a given differential operator, you can get the matrix representation using the matrix(shape) method, e.g. for a small 1D grid of 10 points:

d2_dx2 = Diff(0, dx)**2
mat = d2_dx2.matrix((10,))  # this method returns a scipy sparse matrix
print(mat.toarray())

has the output

[[ 2. -5.  4. -1.  0.  0.  0.]
 [ 1. -2.  1.  0.  0.  0.  0.]
 [ 0.  1. -2.  1.  0.  0.  0.]
 [ 0.  0.  1. -2.  1.  0.  0.]
 [ 0.  0.  0.  1. -2.  1.  0.]
 [ 0.  0.  0.  0.  1. -2.  1.]
 [ 0.  0.  0. -1.  4. -5.  2.]]

If you have periodic boundary conditions, the matrix looks like that:

d2_dx2 = Diff(0, dx, periodic=True)**2
mat = d2_dx2.matrix((10,))  # th

Findiff

Install / Use

README