examples-basic.rst

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)