Approximations of higher-order derivatives

[asross]

Is your feature request related to a problem? Please describe.

I'm working on a project where it's important to estimate higher-order derivatives (e.g. 2nd, 3rd, 4th, and potentially mixed) of quantities in xarray datasets. xarray only has a helper for first derivatives (from #2398).

Describe the solution you'd like

I'd like to be able to call differentiate with a list of variables (e.g., data_array.differentiate(['x','x','y','y'])) and get the most accurate finite difference approximation of the derivative of the data array with respect to all of those variables.

Describe alternatives you've considered

It's tempting to chain together differentiate to approximate higher-order derivatives (e.g. data_array.differentiate('x').differentiate('x').differentiate('y').differentiate('y')). However, although I'm not an expert in finite difference approximations, I'm not sure if this is equivalent to computing the most accurate nth-order central (or average of forward and backward) differences, especially wrt. edges or with irregular grids. I could be wrong, though.

Additional context

The extent of my knowledge of how to estimate higher-order differences.

[TomNicholas]

estimate higher-order derivatives

You could wrap scipy's function for nth-order derivatives, which uses central differencing. I think that would be an appropriate generalisation to include in xarray.

get the most accurate finite difference approximation

AFAIK there isn't really a single "best" way to do derivatives in general though - you might want forward, backwards, central, etc. depending on the use case. If someone wanted to wrap a bunch of options for numerical finite differences then this might be a sensible place to register an xarray accessor, so that that functionality could live outside xarray whilst still being easily accessible.

derivative of the data array with respect to all of those variables.

It's tempting to chain together differentiate

I'm not an expert either, but it seems to me (based on these notes) that mixed partial derivatives could in theory be computed by chaining together single variable derivatives, and it would be analytically equivalent, but you would be wasting a lot of computation because most of the terms cancel out (see page 122).

There are other complications with taking derivatives too, such as "what grid does the data lie on" (e.g. Arakawa grids), "what should the boundary conditions be", and "is there a metric that encodes the distance between points". These kind of complications are addressed in the xGCM library, which you might be interested in.