pythonnumpyscipy# Python finite difference functions?

I've been looking around in Numpy/Scipy for modules containing finite difference functions. However, the closest thing I've found is `numpy.gradient()`

, which is good for 1st-order finite differences of 2nd order accuracy, but not so much if you're wanting higher-order derivatives or more accurate methods. I haven't even found very many specific modules for this sort of thing; most people seem to be doing a "roll-your-own" thing as they need them. So my question is if anyone knows of any modules (either part of Numpy/Scipy or a third-party module) that are specifically dedicated to higher-order (both in accuracy and derivative) finite difference methods. I've got my own code that I'm working on, but it's currently kind of slow, and I'm not going to attempt to optimize it if there's something already available.

Note that I am talking about finite differences, not derivatives. I've seen both `scipy.misc.derivative()`

and Numdifftools, which take the derivative of an analytical function, which I don't have.

Solution

One way to do this quickly is by convolution with the derivative of a gaussian kernel. The simple case is a convolution of your array with `[-1, 1]`

which gives exactly the simple finite difference formula. Beyond that, `(f*g)'= f'*g = f*g'`

where the `*`

is convolution, so you end up with your derivative convolved with a plain gaussian, so of course this will smooth your data a bit, which can be minimized by choosing the smallest reasonable kernel.

```
import numpy as np
from scipy import ndimage
import matplotlib.pyplot as plt
#Data:
x = np.linspace(0,2*np.pi,100)
f = np.sin(x) + .02*(np.random.rand(100)-.5)
#Normalization:
dx = x[1] - x[0] # use np.diff(x) if x is not uniform
dxdx = dx**2
#First derivatives:
df = np.diff(f) / dx
cf = np.convolve(f, [1,-1]) / dx
gf = ndimage.gaussian_filter1d(f, sigma=1, order=1, mode='wrap') / dx
#Second derivatives:
ddf = np.diff(f, 2) / dxdx
ccf = np.convolve(f, [1, -2, 1]) / dxdx
ggf = ndimage.gaussian_filter1d(f, sigma=1, order=2, mode='wrap') / dxdx
#Plotting:
plt.figure()
plt.plot(x, f, 'k', lw=2, label='original')
plt.plot(x[:-1], df, 'r.', label='np.diff, 1')
plt.plot(x, cf[:-1], 'r--', label='np.convolve, [1,-1]')
plt.plot(x, gf, 'r', label='gaussian, 1')
plt.plot(x[:-2], ddf, 'g.', label='np.diff, 2')
plt.plot(x, ccf[:-2], 'g--', label='np.convolve, [1,-2,1]')
plt.plot(x, ggf, 'g', label='gaussian, 2')
```

Since you mentioned `np.gradient`

I assumed you had at least 2d arrays, so the following applies to that: This is built into the `scipy.ndimage`

package if you want to do it for ndarrays. Be cautious though, because of course this doesn't give you the full gradient but I believe the product of all directions. Someone with better expertise will hopefully speak up.

Here's an example:

```
from scipy import ndimage
x = np.linspace(0,2*np.pi,100)
sine = np.sin(x)
im = sine * sine[...,None]
d1 = ndimage.gaussian_filter(im, sigma=5, order=1, mode='wrap')
d2 = ndimage.gaussian_filter(im, sigma=5, order=2, mode='wrap')
plt.figure()
plt.subplot(131)
plt.imshow(im)
plt.title('original')
plt.subplot(132)
plt.imshow(d1)
plt.title('first derivative')
plt.subplot(133)
plt.imshow(d2)
plt.title('second derivative')
```

Use of the `gaussian_filter1d`

allows you to take a directional derivative along a certain axis:

```
imx = im * x
d2_0 = ndimage.gaussian_filter1d(imx, axis=0, sigma=5, order=2, mode='wrap')
d2_1 = ndimage.gaussian_filter1d(imx, axis=1, sigma=5, order=2, mode='wrap')
plt.figure()
plt.subplot(131)
plt.imshow(imx)
plt.title('original')
plt.subplot(132)
plt.imshow(d2_0)
plt.title('derivative along axis 0')
plt.subplot(133)
plt.imshow(d2_1)
plt.title('along axis 1')
```

The first set of results above are a little confusing to me (peaks in the original show up as peaks in the second derivative when the curvature should point *down*). Without looking further into how the 2d version works, I can only really recomend the 1d version. If you want the magnitude, simply do something like:

```
d2_mag = np.sqrt(d2_0**2 + d2_1**2)
```

- Activating Anaconda Environment in VsCode
- How to implement the Softmax function in Python?
- Dataframe - Select the minimum set of rows to cover all possible values of each columns
- Non-equi join in polars
- Tkinter and matplotlib - NavigationToolbar2Tk - cursor position precision
- Media Image is not loaded in the output of Django Project
- selenium upload part not working on windows. Path related?
- How to encode UTF8 filename for HTTP headers? (Python, Django)
- Split columns in Csv file
- In the Django admin site, how do I change the display format of time fields?
- Condition in pytest fixture depending on some test-suite results
- How to make a class JSON serializable
- Apache Airflow unit and integration test
- TypeError: 'MultiPolygon' object is not iterable with Plotly Choropleth
- How do I get the current time in milliseconds in Python?
- Python: Decoding base64 encoded strings
- How to fix the error when try to plot data with pandas?
- Python Image - Finding largest branch from image skeleton
- Enumerate is Returning A None Type
- Python finite difference functions?
- How to extract hex color codes from a colormap
- Regex: Get "London, UK" from "Evolution Recruitment (Agency) (London, UK)"
- Visual Studio Code - How to add multiple paths to python path?
- What are the Spark transformations that causes a Shuffle?
- How to find the second lowest lists into a nested list by their second value?
- python-docx: Parse a table to Panda Dataframe
- Deposit and Withdraw on Account
- Why does this code snippet give the error list object is not callable?
- How to find out if argparse argument has been actually specified on command line?
- Print current call stack from a method in code