CODEWITHSUNDEEP
scipyinterpolationdifferentiation
annon
annon

Reputation: 63

How to get the first and second derivatives after interpolating with scipy.interpolate.CloughTocher2DInterpolator?

I came across this really cool demonstration of CloughTocher2DInterpolator as part of the response to this question: How to get a non-smoothing 2D spline interpolation with scipy, and am interested in getting the derivatives along the x and y axis. According to scipy.org, the extrapolant is continuously differentiable.

Couldn't find any help in the documentation. Getting dx and dy seems more straightforward with scipy's bivariatesplines but cloughTocher seems to fit my purpose best.

Any help appreciated!

Here is the example code from the link:

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

import numpy as np
from scipy.interpolate import CloughTocher2DInterpolator
from scipy.spatial import Delaunay

# Example unstructured mesh:
nodes = np.array([[-1.        , -1.        ],
       [ 1.        , -1.        ],
       [ 1.        ,  1.        ],
       [-1.        ,  1.        ],
       [ 0.        ,  0.        ],
       [-1.        ,  0.        ],
       [ 0.        , -1.        ],
       [-0.5       ,  0.        ],
       [ 0.        ,  1.        ],
       [-0.75      ,  0.4       ],
       [-0.5       ,  1.        ],
       [-1.        , -0.6       ],
       [-0.25      , -0.5       ],
       [-0.5       , -1.        ],
       [-0.20833333,  0.5       ],
       [ 1.        ,  0.        ],
       [ 0.5       ,  1.        ],
       [ 0.36174242,  0.44412879],
       [ 0.5       , -0.03786566],
       [ 0.2927264 , -0.5411368 ],
       [ 0.5       , -1.        ],
       [ 1.        ,  0.5       ],
       [ 1.        , -0.5       ]])

# Theoretical function:
def F(x, y):
    return x + y -  x*y - (x*y)**2 - 2*x*y**2 + x**2*y + 3*np.exp( -((x+1)**2 + (y+1)**2)*5 )

z = F(nodes[:, 0], nodes[:, 1])

# Finer regular grid:
N2 = 19
x2, y2 = np.linspace(-1, 1, N2), np.linspace(-1, 1, N2)
X2, Y2 = np.meshgrid(x2, y2)

# Interpolation:
tri = Delaunay(nodes)
CT_interpolator = CloughTocher2DInterpolator(tri, z)
z_interpolated = CT_interpolator(X2, Y2)

# Plot
fig = plt.figure(1, figsize=(8,14))

ax = fig.add_subplot(311, projection='3d')
ax.scatter3D(nodes[:, 0], nodes[:, 1], z, s=15, color='red', label='points')

ax.plot_wireframe(X2, Y2, z_interpolated, color='black', label='interpolated')
plt.legend();

Upvotes: 5

Views: 689

Answers (1)

tfnn
tfnn

Reputation: 1

First derivatives (gradient) values at nodes are accessible as CT_interpolator.grad.

Derivatives away from nodes appear to not be accessible using scipy, rather unfortunate because the interpolating spline coefficients needed are stored somewhere.

I would suggest that your best option is to use finite differences if the derivatives of the interpolation is what you want:

h = 2e-8
grad_x = (0.5 * CT_interpolator(X2 + h, Y2) - 0.5 * CT_interpolator(X2 - h, Y2))/h
grad_y = (0.5 * CT_interpolator(X2, Y2 + h) - 0.5 * CT_interpolator(X2, Y2 - h))/h

Since the interpolation is cubic, the second order accurate central difference should be fairly accurate.

For second derivatives, you would use finite differences for second derivatives. (For a second-order-accurate, second-derivative finite difference, you can simply apply the same second-order-accurate first-derivative central differences to grad_x and grad_y.)

Upvotes: 0

Related Questions