有一組4096長度的資料,需要找到一階導數從正到負的點,和三階導數從負到正的點,擷取了一小段。
394.0
388.0
389.0
388.0
388.0
392.0
393.0
395.0
395.0
394.0
394.0
390.0
392.0
按照之前所了解的,對離散值求導其實就是求差分,例如第i點的導數(差分)為: y(
s)i=
c−my
i−m+
c−m+
1yi−
m−1+
...+
cm−1
yi+m
−1+c
myi+
m 即在乙個寬度為2m+1的視窗內通過計算前後m個值加權後的和得到。但是在實際使用過程中效果不是很好。於是想到了同樣在乙個寬度為2k+1的視窗內,將這2k+1個點擬合成乙個函式,然後求導就可以得到任意階數的導數值。
首先是函式擬合,使用from scipy.optimize import leastsq即最小二乘擬合
from scipy.optimize import leastsq
class
search
(object):
def__init__
(self, filename):
self.filename = filename
deffunc
(self, x, p):
f = np.poly1d(p)
return f(x)
defresiduals
(self, p, x, y, reg):
regularization = 0.1
# 正則化係數lambda
ret = y - self.func(x, p)
if reg == 1:
return ret
defleastsquare
(self, data, k=100, order=4, reg=1, show=1):
# k為求導視窗寬度,order為多項式階數,reg為是否正則化
l = self.len
step = 2 * k + 1
p = [1] * order
for i in range(0, l, step):
if i + step < l:
y = data[i:i + step]
x = np.arange(i, i + step)
else:
y = data[i:]
x = np.arange(i, l)
try:
r = leastsq(self.residuals, p, args=(x, y, reg))
except:
print("error - curve_fit failed")
fun = np.poly1d(r[0]) # 返回擬合方程係數
df_1 = np.poly1d.deriv(fun) # 求得導函式
df_2 = np.poly1d.deriv(df_1)
df_3 = np.poly1d.deriv(df_2)
df_value = df_1(x)
df3_value = df_3(x)
fun = np.poly1d(r[0]),fun返回的是乙個 polynomial class,具體使用可以見官方文件numpy.poly1d
polynomial物件可以使用deriv方法求導數,求得的依然是 polynomial物件。 df_value = df_1(x)所得到的就是x這個幾個點求得的導數值。
看似大功告成,但是求導的結果並不是很好,如下圖,實際最高點在100左右,但是擬合出來的曲線最高點在120左右,而原因在於使用多項式擬合很難準確擬合曲線。
於是想用高斯函式來實現對曲線的擬合,在matlab中試了下,三階高斯擬合可以很好的擬合曲線,
但是numpy以及sicpy中沒有找到類似poly1d這種物件,雖然可以自己定義高斯函式,如下
def
gaussian
(self, x, *param):
fun = param[0]*np.exp(-np.power(x - param[2], 2.) / (2 * np.power(param[4], 2.)))+param[1]*np.exp(-np.power(x - param[3], 2.) / (2 * np.power(param[5], 2.)))
return fun
但是,在通過最小二乘擬合得到函式引數後只能得到擬合後的點,無法直接求導數..所以並不適合。
所以還是只能回到多項式擬合,如果4階多項式不能表徵的話,更高階的呢
總體來說,效果還是可以接受的。
如果下階段找到好的高斯函式擬合方法,會繼續更新。
mathematica 驗證離散導數格式
數值計算中,需要用差分代替差商,離散常微分方程或者偏微分方程的導數。比如已知 f f f 求 f i 使得 f i 0 近似成立,誤差為四階精度 已知 f f 2 f f 使上述條件成立 如何驗證該公式,在mathematica中輸入如下命令 a1 series f i h a2 series f ...
序列的 離散化
問題描述 對於乙個資料元素很大的序列 a 1 a 2 a n 我們在處理這個序列時,往往只會關心序列中元素的大小關係。這時候為了處理方便,可把序列對映為1.n範圍內的數,其中最小的元素對映為 第 小的數對映為2,例如 原 序 列 19827345 2000000053 73925 98203456 ...
導數法求引數範圍時注意問題
2019屆鳳翔中學高三理科單調性與最值課時作業第17題改編 已知 f x cfrac x neq a 當 a 0 時,函式 f x 在區間 1,infty 上單調遞減,求實數 a 的範圍。法1 影象法,f x 1 cfrac 其對稱中心為 a,1 用變換法做出其函式影象如下,由影象可得,要使得函式 ...