2017-02-17 303 views
0

我試圖在scipy中使用splev的幾個點找到樣條函數的衍生物。例如:樣條函數的衍生物:`scipy splev`

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] 

我想在指定點繪製切線:

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

我覺得這些切線都是錯誤的。我的理解是,xp_numyp_num是相對於xy的樣條的數值導數。所以要找到dy/dx,我應該簡單地將它們分開。它是否正確?

最後,我想找到一個曲線的切線像this

+1

wh在「正常化」嗎? –

+0

@PaulPanzer:它規範化每一行。我加了'normalize'的定義。 – Mahdi

+0

奇怪的是,你的代碼在我看來好像它不應該在劇情中創建任何向量。 'yp_num'和'xp_num'是1d,不是嗎?所以你的'R'應該是1x東西,所以當你切割R [1:,0]時,你應該得到一個空數組。我在這裏錯過了什麼嗎? –

回答

1

您的問題(的明顯錯誤的衍生工具),因爲你不使用它們,至少在你的代碼發佈不相關的數值導。什麼顯然是錯誤的,除非你的normalized功能做一些真正神奇的是你將yp_the通過xp_the,因爲前者是確實的增量,後者是不是應該不斷得到

dy 
-- 
dx 

,而不是你的

dy 
---- 
x dx 

你很可能從公式來進行使用的參數曲線

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

t=x,然後忽略dx/dx是恆定的。發生在我們最好的事情。

+0

太棒了,完全正確我想過參數曲線!謝謝! – Mahdi

1

你沒有包括您R = normalized(yp_the/xp_the)

我linalg.norm

然後我改變了對delta_Y的標準化衍生

做到了,放棄了對顫動

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

+0

對不起我的錯誤。我打算用'xp_num'和'yp_num'來繪製切線。 – Mahdi