CODEWITHSUNDEEP
pythonnumpyscipycurve-fitting
Mahdi
Mahdi

Reputation: 3238

Derivatives of a spline: `scipy splev`

I'm trying to find derivatives of a spline at several points using splev in scipy. For example:

import numpy as np
from scipy.interpolate import splprep, splev
import matplotlib.pyplot as plt

# function to normalize each row
def normalized(a, axis=-1, order=2):
    l2 = np.atleast_1d(np.linalg.norm(a, order, axis))
    l2[l2==0] = 1
    return a / np.expand_dims(l2, axis)

# points and spline
pts = np.array([[0,0],[1,1],[2,np.sqrt(2)],[4,2],[9,3]])
tck, u = splprep(pts.T, u=None, k=3, s=0.0, per=0)

# compute new points and derivatives
u_new = np.linspace(u.min(), u.max(), 5*u.shape[0])
x_new, y_new = splev(u_new, tck, der=0)
xp_num, yp_num = splev(pts, tck, der=1) # numerical derivatices
xp_the, yp_the= pts[1:,0], 0.5/np.sqrt(pts[1:,0]) # analytical derivatives
R = normalized(yp_num/xp_num)
X,Y = pts[1:,0], pts[1:,1]
U,V = X + R[1:,0], Y+R[1:,1]

I'd like to plot the tangents at the given point:

plt.plot(x_new,y_new,'-b')
plt.plot(pts[:,0],pts[:,1],'--ro')
plt.quiver(X,Y,U,V, angles='xy', scale_units='xy')

enter image description here

I think these tangent lines are wrong. My understanding was that xp_num and yp_num are numerical derivatives of the spline with respect to x and y. So to find dy/dx, I should simply divide them. Is this correct?

Eventually, I'd like to find tangents of a curve like this

Upvotes: 0

Views: 2072

Answers (2)

Paul Panzer
Paul Panzer

Reputation: 53089

Your problem (the obviously wrong derivatives) is not related to the numerical derivative since you are not using them at least in the code you posted. What is clearly wrong unless your normalized function does something truly magic is your dividing yp_the by xp_the since the former is indeed the increment, the latter is not it should be constant to get

dy
--
dx

as opposed to your

 dy
----
x dx

You were probably coming from the formula for a parametric curve

.     dy
      --
dy    dt
-- = ----
dx    dx
      --
      dt

used t=x and then overlooked that dx/dx is constant. Kind of thing that happens to the best of us.

Upvotes: 2

f5r5e5d
f5r5e5d

Reputation: 3741

you didn't include your R = normalized(yp_the/xp_the) source

I did it with linalg.norm

then I changed the delta_Y for the normalized derivative

and gave up on quiver

import numpy as np
from scipy.interpolate import splprep, splev
import matplotlib.pyplot as plt

# points and spline
pts = np.array([[0,0],[1,1],[2,np.sqrt(2)],[4,2],[9,3]])
tck, u = splprep(pts.T, u=None, k=3, s=0.0, per=0)

# compute new points and derivatives
u_new = np.linspace(u.min(), u.max(), 5*u.shape[0])
x_new, y_new = splev(u_new, tck, der=0)
xp_num, yp_num = splev(pts, tck, der=1) # numerical derivatices
xp_the, yp_the= pts[1:,0], 0.5/np.sqrt(pts[1:,0]) # analytical derivatives
#R = normalized(yp_the/xp_the)
N = np.linalg.norm(np.array([xp_the, yp_the]), axis=0)

X,Y = pts[1:,0], pts[1:,1]
#U,V = X + R[1:,0], Y+R[1:,1]
U,V = X + xp_the/N, Y + X*yp_the/N # delta Y = dy/dx * x

plt.axes().set_aspect('equal', 'datalim')

plt.plot(x_new,y_new,'-b')
plt.plot(pts[:,0],pts[:,1],'--ro')
#plt.quiver(X,Y,U,V, scale=10, pivot='mid')# angles='xy', scale_units=’xy’, scale=1

plt.plot((X, U), (Y, V), '-k', lw=3)

enter image description here

Upvotes: 1

Related Questions