问题描述
我正在尝试计算添加白噪声的洛伦兹系统模拟的导数。例如，如果我想计算 x 数据 w.r.t. 的导数。时间 t。 我故意添加了一些噪音来检查导数计算器对实验数据的稳健性。所以，我不能使用诸如有限差分之类的数值导数技术。 我选择了具有更好条件多项式（例如切比雪夫多项式）的多项式回归来获得多项式方程，然后对其进行简单微分。 当然，对于这么复杂的曲线，我需要做多次局部多项式拟合，才能捕捉到曲线局部的剧烈变化。 但是，我从多项式回归中获得了非常高的导数值。请看一下 Python 代码。运行代码需要 5 分钟。
import numpy as np # matrix calc
import matplotlib.pyplot as plt
from methods import lorenz
from scipy.integrate import odeint # scientific computation (ode solver)
import scipy.sparse as sparse
# %% Initializing the chaotic Lorenz-system data
print('Initializing the chaotic Lorenz-system data')
sigma = 10
beta = 8/3
rho = 28
dt = 0.01
tspan = np.arange(0,50.01,dt,dtype=float) # time span
x0 = np.array([-8,8,27]) # initial conditions
x = odeint(lambda x,t: lorenz(x,t,sigma,beta,rho),x0,tspan,atol=1e-10,rtol=1e-10)
rows = np.shape(x)[0]
col = np.shape(x)[1]
# Visualizing the data
plt.figure(1)
plt.plot(tspan,x[:,0])
# %% True derivative using clean data
print('True derivative using clean data')
deriv = np.zeros((rows,col))
for i in range(1,rows): # for each row
deriv[i,:] = lorenz(x[i,:],rho)
# %% Tikhonov Differentiation (for clean data)
# returns argmin of (g) |Ag-f|_2^2 + Lambda*|Dg|_2^2
# d nth derivative
# f The data
# dx/dt choose step length dx/ time step dt. For uniform grid only such as the rectangular grid
# Lambda The hyperparameter or the trust over the regularization term
# Intializing paramters
print('Tikhonov Differentiation (for clean data)')
Lambda = 0.025 # threshold
# Addition and subtraction won't affect the derivative/slope but simplify the computation
f = np.array(x-x[0])
g = np.zeros((rows,col): # for each columns
A = np.zeros((rows,rows))
for j in range(1,rows):
A[j,j] = dt/2 # for diagonal elements
A[j,0] = dt/2 # For first column only
for k in range(1,j):
A[j,k] = dt
# Since k<=j,=> A is a lower triangular matrix
e = np.ones(rows-1)
D = sparse.diags([e,-e],[1,0],shape=(rows-1,rows)).todense()/dt
# Invert to compute the derivative
g[:,i] = np.linalg.lstsq(A.T.dot(A)+Lambda*D.T.dot(D),A.T.dot(x[:,i]),rcond=None)[0]
# %% polynomial differentiation at data
# u The data near the point
# x variable with respect to which the differentiation is to be done
# deg degree of polynomial to use
# diff maximum order of derivative we want
# width number of data points neighbouring to the data point of interest
# Since the derivatives in dynamical systems are really complex in nature,therefore choosing the whole
# data for polynomial derivative might not capture the real derivative,therefore a series of
# polynomial derivatives are taken with their 5 neighbouring data points
# to capture the local gradient/ slope.
print('polynomial differentiation at data')
rows = np.shape(x)[0]
#index = int((rows-1)/2) # where you want to compute your derivative
deg = 3
diff = 1
width = 100 # determined on the basis of true data
# fit to a Chebyshev polynomial as they are better conditioned than normal polynomials
# so small purturbation in the coefficients won't result in drastic change of the solution
# Making width array
width_array=np.arange(0,rows,width)
derivatives=np.zeros((rows,3))
# Here,width_array contains index of the input,such that index and index+1 accommodates exactly width number of data points
# inside the loop width_index is used to loop over each element of the width_array
# so
#np.shape(width_array)[0]-1 otherwise upper_bound would exceed the array size
for width_index in range(0,np.shape(width_array)[0]-1):
#lower_bound=x[width_array[width_index]][0]
#upper_bound=x[width_array[width_index+1]][0]
lower_bound=width_array[width_index]
upper_bound=width_array[width_index+1]
# here lower_bound and upper_bound gives the index of input data in between which width number of datapoint lies
u = np.transpose(x[lower_bound:upper_bound,0])
time= np.transpose(tspan[lower_bound:upper_bound])
poly = np.polynomial.chebyshev.Chebyshev.fit(time,u,deg)
# Take derivative
derivatives_window = np.zeros((width))
for l in range(0,width-1): # for each data point in the current window
for d in range(1,diff+1): # diff+1 is used so that d goes to diff
derivatives_window[l] = (poly.deriv(m=d)(u[l])) # consecutive differentiation
derivatives[lower_bound:upper_bound,0]=derivatives_window
# %% Plotting
print('Plotting')
plt.figure(2)
#for i in range(0,col):plt.plot(tspan,deriv[:,i])
plt.plot(tspan,"black",g[:,'blue',derivatives[:,"green")
plt.legend(['True','Tikhonov','polynomial'])
# Conclusion: Here we can't fit 1 curve as 1 curve can't capture the complex shape of the dynamics.
解决方法
暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！
如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。
小编邮箱:dio#foxmail.com (将#修改为@）