findiff Documentation

findiff is a Python package designed for computing numerical derivatives and solving partial differential equations (PDEs) in any number of dimensions.

Key capabilities include:

• Numerical Differentiation: Differentiate arrays of any dimension along any axis with configurable accuracy orders. Supports both uniform and non-uniform grids.
• Vector Calculus: Built-in support for standard operators like gradient, divergence, and curl.
• Advanced Operators: Handle arbitrary linear combinations of derivatives with constant or variable coefficients. Supports arbitrary stencils and compact (implicit) finite differences.
• PDE Solving: Solve PDEs with Dirichlet, Neumann, Robin, or Periodic boundary conditions. Supports time-dependent PDEs via the Method of Lines (e.g., Forward Euler, RK4, Backward Euler, Crank-Nicolson).
• Performance: Fully vectorized operations with support for GPU, JAX, and CuPy backends. Can be combined with jax.jit for acceleration.
• Mathematical Tools: Generate matrix representations of operators, solve eigenvalue problems (e.g., Schrodinger equation), and calculate raw finite difference coefficients.

The BoundaryConditions object allows you to specify different types of boundary constraints on your grid edges.

Dirichlet Boundary Conditions

Set a single value to specify the value of the function at the boundary. bc[index] = value

Neumann Boundary Conditions

Specify the derivative at the boundary by passing a 2-tuple: (diff_operator, value). bc[index] = (Diff(axis, dx), value)

Robin (Mixed) Boundary Conditions

Robin conditions follow the form $\alpha u + \beta \frac{\partial u}{\partial n} = g$. You can specify these in two ways:

1. As a 4-tuple: (alpha, diff_op, beta, g)
2. As a 2-tuple: Construct a custom operator using Identity() and Diff, then pass (operator, g).

2D Boundary Conditions

For 2D grids, use slicing to apply conditions to entire edges (e.g., bc[0, :] for the top edge).

# Neumann Example
bc[0] = 0                         # u(0) = 0
bc[-1] = Diff(0, dx), 1.0         # u'(L) = 1.0

# Robin Example (4-tuple)
bc[-1] = (1, Diff(0, dx), 1, 3)   # u + u' = 3

# Robin Example (Custom Operator)
from findiff import Identity
robin_op = alpha * Identity() + beta * Diff(0, dx)
bc[-1] = robin_op, g

Standard (explicit) finite differences express a derivative as a weighted sum of function values:

$$f'i = \sum_k c_k , f{i+k}$$

Compact (implicit) finite differences generalize this by including derivative values on the left-hand side of the equation:

$$\sum_k \alpha_k , f'{i+k} = \sum_k c_k , f{i+k}$$

Key Characteristics:

• High Accuracy: Compact schemes (like the Lele, 1992 tridiagonal scheme) can achieve much higher accuracy (e.g., 6th-order) with a much smaller stencil width than explicit schemes.
• Computational Cost: Applying a compact operator requires solving a banded linear system (such as a tridiagonal system via the Thomas algorithm). This has a modest $O(N)$ overhead.
• Resolution: They resolve short-wavelength features significantly more accurately than explicit schemes of the same stencil width.

For implementation details, refer to the compact-fd guide.

You can implement vector calculus operators in non-Cartesian coordinate systems (like polar coordinates) by assembling general linear combinations of Diff operators with variable coefficients.

For example, to implement the 2D polar Laplacian $\nabla^2 = \frac{\partial^2}{\partial r^2} + \frac{1}{r}\frac{\partial}{\partial r} + \frac{1}{r^2}\frac{\partial^2}{\partial \varphi^2}$, you combine Diff instances for the radial ($r$) and angular ($\varphi$) components using the appropriate variable coefficients ($1/R$ and $1/R^2$).

import numpy as np
from findiff import Diff

# Setup grid
r = np.linspace(0.1, 10, 100)
phi = np.linspace(0, 2*np.pi, 100, endpoint=False)
dr, dphi = r[1] - r[0], phi[1] - phi[0]
R, Phi = np.meshgrid(r, phi, indexing='ij')

# Define function f(r, phi) = r^2
f_polar = R**2

# Assemble the polar Laplacian operator
laplace_polar = (
    Diff(0, dr)**2
    + (1/R) * Diff(0, dr)
    + (1/R**2) * Diff(1, dphi)**2
)

# Apply the operator
result = laplace_polar(f_polar)

To solve time-dependent PDEs of the form $\frac{\partial u}{\partial t} = L(u)$, use the TimeDependentPDE class. This implements the Method of Lines (MOL) by discretizing the spatial operator $L$ with findiff finite difference operators and advancing the resulting system of ODEs in time using built-in time-stepping methods.

To set up a problem, you need:

1. A spatial operator L (constructed using Diff).
2. An initial condition u0 (an array representing the state at $t=0$).
3. BoundaryConditions object defining the behavior at the edges of your grid.
4. A time array t defining the temporal discretization.

import numpy as np
from findiff import Diff, TimeDependentPDE, BoundaryConditions

# Setup spatial grid and operator
nx = 101
x = np.linspace(0, 1, nx)
dx = x[1] - x[0]
D = 0.01
L = D * Diff(0, dx)**2

# Initial condition
u0 = np.sin(np.pi * x)

# Boundary conditions
bc = BoundaryConditions((nx,))
bc[0] = 0
bc[-1] = 0

# Time steps
t = np.linspace(0, 1, 500)

# Solve
pde = TimeDependentPDE(L, u0, bc, t)
u_final = pde.solve()  # returns the solution at t=1

The findiff package approximates differential operators numerically using finite difference schemes.

1D Case

In 1D, the $n$-th derivative of a function $f$ at a grid point $x_k$ is approximated by a linear combination of function values in a surrounding stencil:

$$\left(\frac{d^n f}{dx^n}\right)k \approx \sum{j \in A} c_{j} f_{k+j}$$

Where $A$ is a set of offsets. findiff supports:

• Symmetric stencils: Uses points on both sides of $x_k$ (e.g., $p=q=1$). These are ideal for interior points but cannot be used at boundaries.
• One-sided stencils: Uses points on only one side (e.g., forward or backward stencils) to handle boundary conditions where neighbors are missing.

Multiple Dimensions

For multi-dimensional functions, partial derivatives are calculated using stencils that are typically superpositions of 1D stencils. For example, a 2D Laplacian ($\frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2}$) uses a cross-shaped stencil composed of two 1D stencils.