Getting Started

Installation

pip install --upgrade findiff

Derivatives

First Derivatives

The first derivative along the 0-th axis (”\(x_0\)-axis”),

\[\frac{\partial}{\partial x_0}\quad,\]

can be defined by

from findiff import FinDiff

d_dx = FinDiff(0, dx)

The first argument is the axis along which to take the partial derivative. The second argument is the spacing of the (equidistant) grid along that axis.

Accordingly, the first partial derivative with respect to the k-th axis

\[\frac{\partial}{\partial x_k}\]

is

FinDiff(k, dx_k)

Let’s initialize a one-dimensional array f with some values, for example:

import numpy as np
x = np.linspace(-np.pi, np.pi, 100)
dx = x[1] - x[0]
f = np.sin(x)

and calculate the first derivative with respect of the zeroth axis.

FinDiff objects behave like operators, so in order to apply them, you can simply call them on a numpy ndarray of any shape:

d_dx = FinDiff(0, dx)
df_dx = d_dx(f)

Now df_dx is a new numpy array with the same shape as f containing the first derivative with respect to the zeroth axis:

../_images/get_started_plot_1.png

Higher Derivatives

The n-th partial derivatives, say with respect to \(x_k\),

\[\frac{\partial^n}{\partial x_k^n}\]

is

FinDiff(k, dx_k, n)

where the last argument stands for the degree of the derivative.

A mixed partial derivatives like

\[\frac{\partial^3}{\partial x^2 \partial y}\]

is defined by

FinDiff((0, dx, 2), (1, dy, 1))

where for each partial derivative, there is a tuple of the form (axis, spacing, degree) in the argument list.

General Differential Operators

FinDiff objects can be combined to describe general differential operators. For example, the wave operator

\[\frac{1}{c^2}\frac{\partial^2}{\partial t^2} - \frac{\partial^2}{\partial x^2}\]

can be written as

1 / c**2 * FinDiff(0, dt, 2) - FinDiff(1, dx, 2)

if the 0-axis represents the t-axis and the 1-axis the x-axis.

Non-constant coefficients must be wrapped as Coef objects. For instance,

\[x^2 \frac{\partial^2}{\partial x^2}\]

is written as

x = np.linspace(-1, 1, 21)
Coef(x) * FinDiff(0, dx, 2)

Finally, multiplication of two FinDiff objects means chaining differential operators, for example

\[\left(\frac{\partial}{\partial x} - \frac{\partial}{\partial y}\right) \cdot \left(\frac{\partial}{\partial x} + \frac{\partial}{\partial y}\right) = \frac{\partial^2}{\partial x^2} - \frac{\partial^2}{\partial y^2}\]

or in findiff:

d_dx = FinDiff(0, dx, 1)
d_dy = FinDiff(1, dx, 1)

(d_dx - d_dy) * (d_dx + d_dy)

Accuracy Control

By default, findiff uses finite difference schemes with second order accuracy in the grid spacing. Higher orders can be selected by setting the keyword argument acc, e.g.

FinDiff(0, dx, 2, acc=4)

for fourth order accuracy.

Finite Difference Coefficients

findiff uses finite difference schemes to calculate numerical derivatives. If needed, the finite difference coefficients can be obtained from the coefficients function, e.g. for second derivative with second order accuracy:

from findiff import coefficients
coefficients(deriv=2, acc=2)

which yields

{'backward': {'accuracy': 2,
              'coefficients': array([-1.,  4., -5.,  2.]),
              'offsets': array([-3, -2, -1,  0])},
 'center': {'accuracy': 2,
            'coefficients': array([ 1., -2.,  1.]),
            'offsets': array([-1,  0,  1])},
 'forward': {'accuracy': 2,
             'coefficients': array([ 2., -5.,  4., -1.]),
             'offsets': array([0, 1, 2, 3])}}

Matrix Representations

For a given FinDiff differential operator, you can get the matrix representation using the matrix(shape) method, e.g.

x = [np.linspace(0, 6, 7)]
d2_dx2 = FinDiff(0, x[1]-x[0], 2)
u = x**2

mat = d2_dx2.matrix(u.shape)  # this method returns a scipy sparse matrix
print(mat.toarray())

yields

[[ 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.]]

Of course this also works for general differential operators.

Stencils

Automatic Stencils

When you define a differential operator in findiff, it automatically chooses suitable stencils to apply on a given grid. For instance, consider the 2D Laplacian

\[\frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2}\]

which can be defined (in second order accuracy here) as

laplacian = FinDiff(0, dx, 2) + FinDiff(1, dy, 2)

When you then apply the Laplacian to an array, findiff applies it to all grid points. So depending on the grid point point, findiff chooses on or the other stencil.

You can inspect the stencils for a differential operator by calling its stencil method, passing the shape of the grid

laplacian.stencil(f.shape)

This returns

{('L', 'L'): {(0, 0): 4.0, (0, 1): -5.0, (0, 2): 4.0, (0, 3): -1.0, (1, 0): -5.0, (2, 0): 4.0, (3, 0): -1.0},
 ('L', 'C'): {(0, -1): 1.0, (0, 1): 1.0, (1, 0): -5.0, (2, 0): 4.0, (3, 0): -1.0},
 ('L', 'H'): {(0, -3): -1.0, (0, -2): 4.0, (0, -1): -5.0, (0, 0): 4.0, (1, 0): -5.0, (2, 0): 4.0, (3, 0): -1.0},
 ('C', 'L'): {(-1, 0): 1.0, (0, 1): -5.0, (0, 2): 4.0, (0, 3): -1.0, (1, 0): 1.0},
 ('C', 'C'): {(-1, 0): 1.0, (0, -1): 1.0, (0, 0): -4.0, (0, 1): 1.0, (1, 0): 1.0},
 ('C', 'H'): {(-1, 0): 1.0, (0, -3): -1.0, (0, -2): 4.0, (0, -1): -5.0, (0, 0): 0.0, (1, 0): 1.0},
 ('H', 'L'): {(-3, 0): -1.0, (-2, 0): 4.0, (-1, 0): -5.0, (0, 0): 4.0, (0, 1): -5.0, (0, 2): 4.0, (0, 3): -1.0},
 ('H', 'C'): {(-3, 0): -1.0, (-2, 0): 4.0, (-1, 0): -5.0, (0, -1): 1.0, (0, 0): 0.0, (0, 1): 1.0},
 ('H', 'H'): {(-3, 0): -1.0, (-2, 0): 4.0, (-1, 0): -5.0, (0, -3): -1.0, (0, -2): 4.0, (0, -1): -5.0, (0, 0): 4.0}}

In the interior of the grid (the :code:`(‘C’, ‘C’) case), the stencil looks like this:

../_images/laplace2d.png

The blue points denote the grid points used by the stencil, the tu ple below denotes the offset from the current grid point and the value inside the blue dot represents the finite different coefficient for grid point. So, this stencil evaluates the Laplacian at the center of the “cross” of blue points. Obviously, this does not work near the boundaries of the grid because that stencil would require points “outside” of the grid. So near the boundary, findiff switches to asymmetric stencils (of the same accuracy order), for example

../_images/stencil_laplace2d_corner.png

for a corner ('L', 'L'), or

../_images/stencil_laplace2d_border.png

for the left edge ('L', 'C').

The stencil method works for grids of all dimensions and not just two. But of course, it is not easy to visualize for higher dimensions.

While FinDiff object act on complete arrays, stencils can be applied to individual grid points, if you just need a numeric derivative evaluated at one point. For instance,

x = y = np.linspace(0, 1, 101)
X, Y = np.meshgrid(x, y, indexing='ij')
f = X**3 + Y**3

stencils = laplacian.stencil(f.shape)
stencils.apply(f, (100, 100)) # evaluate at f[100, 100]

returns 12, as expected. stencil returns a list of stencils and when calling apply, the appropriate stencil for the selected grid point is chosen. In the example, it chooses a stencil for the corner point.

Stencils From Scratch

There may be situations when you want to create your own stencils and do not want to use the stencils automatically created by FinDiff objects. This is mainly for exploratory work. For example, you may wonder, how the coefficients for the 2D Laplacian look like if you don’t use the cross-shaped stencil from the previous section but rather an x-shaped one:

../_images/laplace2d-x.png

This can easily be determined with findiff by using the Stencil class directly:

offsets = [(0, 0), (1, 1), (-1, -1), (1, -1), (-1, 1)] # x-shaped offsets
stencil = Stencil(offsets, {(2, 0): 1, (0, 2): 1})

returns

{(0, 0): -2.0, (1, 1): 0.5, (-1, -1): 0.5, (1, -1): 0.5, (-1, 1): 0.5}

The second argument of the Stencil constructor defines the derivative operator:

{(2, 0): 1, (0, 2): 1}

corresponds to

\[1 \cdot \frac{\partial^2}{\partial x_0} + 1 \cdot \frac{\partial^2}{\partial x_1}.\]