Derivatives for a single point in a vector space

[Nickleaton]

I have a vector problem with 5 components.

I would like to calculate a full set of derivatives for a given pair of points. Ideally up to order 4.

The function is numeric and not symbolic.

The step size that is appropriate, I can calculate using the standard h/x0 where h = x0 * sqrt (ulp)
where ulp is the unit of least precision.

For the 4th derivate order, it ties up if values are calculated at -2,-1,0,1 and 2 steps. Based off this https://en.wikipedia.org/wiki/Finite_difference_coefficient

I think I get 4 orders of accurace for first and second. 2 degrees of accuracy for the 3rd and 4th order. That could be changed but that is more than likely to be good enough for my purposes.

I need mixed partial derivatives too. This is to feed into a Taylor expansion. This is for set points in the vector space. A regular grid is not what is needed.

That gives quite a lot of FinDiff objects [780], but easy to create in a loop.

I'm not too worried about efficiency, since the function to be evaluated caches well with standard functools.cache

So here's a test function with just 3 components.

def ff(x, y, z) -> float:
result = x ** 2 + .1 * y * z
print(f"{x:0.4f} {y:0.4f} {z:0.4f} -> {result:0.4f}")
return result

This can be vectorised

fv = np.vectorize(ff, signature="(),(),()->()")

A bit of test code

x, y, z = [np.linspace(0, 1, 4)] * 3

res = fv(x,y,z)

print (res)

gives the right sort of answer

0.0000 0.0000 0.0000 -> 0.0000
0.3333 0.3333 0.3333 -> 0.1222
0.6667 0.6667 0.6667 -> 0.4889
1.0000 1.0000 1.0000 -> 1.1000
[0.0000 0.1222 0.4889 1.1000]

I'm not sure how to proceed.

How do I get a full set of derivatives up to order 4 for the point <1,1,1>?

Thanks

[maroba]

FinDiff applies the derivative operator to the whole grid, not to a single point. But of course, since you say that performance is no problem for you, you can retrieve that single point afterwards from the grid.

For getting the Taylor expansion, you'd have to define all the terms manually, probably something along these lines:

from findiff import FinDiff, Identity
from itertools import permutations, product

xyz = 0, 1, 2

taylor_terms = []

h = 0.1


for deriv in range(1, 5):

    diffs = list(product(*[xyz]*deriv))

    for axes in diffs:
        fd = Identity()
        for ax in axes:        
             fd = FinDiff(ax, h, 1) * fd
        taylor_terms.append(fd)