Python中的曲线拟合微分方程

Question

我有一条 1000 点的曲线，我想将其拟合为一个微分方程，其形式为 x'' = (a*x'' + b x' + c x + d)，其中 a,b,c,d 是常数。在使用Python 2.7时，我将如何进行操作？

Answer 1

你可以使用优化来寻找最佳参数 a, b, c, d 使测量值和预测值之间的差异最小化。下面是一个三阶微分方程的示例代码，它是在 ekko 我开发的。

from gekko import GEKKO

t_data = [0,0.1,0.2,0.4,0.8,1,1.5,2,2.5,3,3.5,4]
x_data = [2.0,1.6,1.2,0.7,0.3,0.15,0.1,\
          0.05,0.03,0.02,0.015,0.01]

m = GEKKO()
m.time = t_data

# states
x = m.CV(value=x_data); x.FSTATUS = 1  # fit to measurement
y,z = m.Array(m.Var,2,value=0)

# adjustable parameters
a,b,c,d = m.Array(m.FV,4)
a.STATUS=1; b.STATUS=1; c.STATUS=1; d.STATUS=1 

# differential equation
#      Original:  x''' = a*x'' + b x' + c x + d
#      Transform: y = x'
#                 z = y'
#                 z' = a*z + b*y + c*x + d
m.Equations([y==x.dt(),z==y.dt()])
m.Equation(z.dt()==a*z+b*y+c*x+d) # differential equation

m.options.IMODE = 5   # dynamic estimation
m.options.NODES = 3   # collocation nodes
m.solve(disp=False)   # display solver output
print(a.value[0],b.value[0],c.value[0],d.value[0])

import matplotlib.pyplot as plt  # plot solution
plt.plot(m.time,x.value,'bo',label='Predicted')
plt.plot(m.time,x_data,'rx',label='Measured')
plt.legend(); plt.xlabel('Time'), plt.ylabel('Value'); plt.show()

大多数微分方程求解器要求你将高阶导数转化为独立的一阶导数方程。这很容易做到，因为你需要为每一个额外的阶数(二阶和三阶导数)定义一个新的状态，作为 此处.

Answer 2

当然，你打算把第三个导数放在右边。

将你的数据分成相对较小的仓，可能会重叠。对于每一个bin，计算一个数据的立方体近似值。由此计算出该组中心点的导数。有了所有组的导数，你现在就有了一个经典的线性回归问题。

如果样本间距相等，你可以尝试通过FFT将问题移到频率空间。对数据进行合理的截断可能是这里的一个问题。在频率空间中，任务简化为多项式线性回归。