Basic Examples of findiff

findiff works in any dimension. But for the sake of demonstration, let’s concentrate on the cases 1D and 3D. We are using uniform, i.e. equidistant, grids here. The non-uniform case will be shown in another notebook.

Preliminaries

Our imports:

import numpy as np
from findiff import Diff, coefficients

Simple 1D Cases

Suppose we want to differentiate two 1D-arrays f and g, which are filled with values from a function

\[f(x) = \sin(x) \quad \mbox{and}\quad g(x) = \cos(x)\]

and we want to take the 2nd derivative. This is easy done analytically:

\[\frac{d^2f}{dx^2} = -\sin(x) \quad \mbox{and}\quad \frac{d^2g}{dx^2} = -\cos(x)\]

Let’s do this numerically with findiff. First we set up the grid and the arrays:

x = np.linspace(0, 10, 100)
dx = x[1] - x[0]
f = np.sin(x)
g = np.cos(x)

Then we construct the derivative object, which represents the differential operator \(\frac{d^2}{dx^2}\):

d_dx = Diff(0, dx)

The first parameter is the axis along which to take the derivative. Since we want to apply it to the one and only axis of the 1D array, this is a 0. The second parameter describes the grid to be used. In our case, we have an equidistant grid point, so a single number (the grid spacing along the desired axis) suffices. Diff always returns a first derivative. If you need higher order derivatives, use exponentiation:

d2_dx2 = Diff(0, dx) ** 2

Then we apply the operator to f and g, respectively:

result_f = d2_dx2(f)
result_g = d2_dx2(g)

That’s it! The arrays result_fand result_g have the same shape as the arrays f and g and contain the values of the second derivatives.

Finite Difference Coefficients

By default the Diff class uses second order accuracy. For the second derivative, it uses the following finite difference coefficients:

coefficients(deriv=2, acc=2)

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

But Diff can handle any accuracy order. For instance, have you ever wondered, what the 10th order accurate coefficients look like? Here they are:

coefficients(deriv=2, acc=10)

{'backward': {'coefficients': array([  -0.53253968,    6.42373016,  -35.55158728,  119.41369042,
         -271.26190464,  439.39444427, -521.11333314,  457.02976176,
         -295.51984119,  138.59325394,  -44.43730158,    7.56162698]),
  'offsets': array([-11, -10,  -9,  -8,  -7,  -6,  -5,  -4,  -3,  -2,  -1,   0])},
 'center': {'coefficients': array([ 3.17460317e-04, -4.96031746e-03,  3.96825397e-02, -2.38095238e-01,
          1.66666667e+00, -2.92722222e+00,  1.66666667e+00, -2.38095238e-01,
          3.96825397e-02, -4.96031746e-03,  3.17460317e-04]),
  'offsets': array([-5, -4, -3, -2, -1,  0,  1,  2,  3,  4,  5])},
 'forward': {'coefficients': array([   7.56162876,  -44.43731776,  138.59331976, -295.52000468,
          457.03003946, -521.1136706 ,  439.39474213, -271.26209495,
          119.41377646,  -35.55161345,    6.42373497,   -0.53254009]),
  'offsets': array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11])}}

Accuracy order

If you want to use for example 10th order accuracy, just tell the Diff constructor to use it:

d2_dx2 = Diff(0, dx, acc=10) ** 2
result = d2_dx2(f)

Simple 3D Cases

Now let’s differentiate a 3D-array f representing the function

\[f(x, y, z) = \sin(x) \cos(y) \sin(z)\]

x, y, z = [np.linspace(0, 10, 100)]*3
dx, dy, dz = x[1] - x[0], y[1] - y[0], z[1] - z[0]
X, Y, Z = np.meshgrid(x, y, z, indexing='ij')
f = np.sin(X) * np.cos(Y) * np.sin(Z)

The partial derivatives \(\frac{\partial f}{\partial x}\) or \(\frac{\partial f}{\partial z}\) are given by

d_dx = Diff(0, dx)
d_dz = Diff(2, dz)

The x-axis is the 0th axis, y, the first, z the 2nd, etc. The mixed partial derivative \(\frac{\partial^2 f}{\partial x \partial y}\) is specified by multiplying the two first order partial derivatives:

d2_dxdy = Diff(0, dx) * Diff(1, dy)
result = d2_dxdy(f)

Of course, the accuracy order can be specified the same way as for 1D.

General Linear Differential Operators

Diff objects can bei added and easily multiplied by numbers. For example, to express

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

we can say

linear_op = Diff(0, dx)**2 + 2 * Diff(0, dx) * Diff(1, dy) + Diff(1, dy)**2

If you want to multiply by variables instead of plain numbers, it works the same way. For example,

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

is

linear_op = X * Diff(0, dx) + Y**2 * Diff(1, dy)

Applying those general operators works the same way as for the simple derivatives:

result = linear_op(f)