Python中的曲线拟合微分方程答案

【问题标题】：Curve fitting differential equations in PythonPython中的曲线拟合微分方程
【发布时间】：2014-05-28 16:42:06
【问题描述】：

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

【问题讨论】：

尝试解析求解，然后针对不同的 a、b、c、d 和初始条件绘制它。任何解析曲线都比任何计算曲线带来更多的时间信息。这种类型的 diff.eq。相当容易解决，并且在任何关于 diff.eq 的书中都有关于这类同质 eq-s 的章节。

标签： python-2.7 physics curve-fitting differential-equations

【解决方案1】：

当然，您打算在右边使用三次导数。

将您的数据分组到相对较小的箱中，可能会重叠。对于每个 bin，计算数据的三次近似值。由此计算组中心点的导数。有了所有组的导数，您现在就有了一个经典的线性回归问题。

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

【讨论】：

【解决方案2】：

您可以使用优化来找到最佳参数a、b、c、d，以最大限度地减少测量值和预测值之间的差异。这是我在gekko 中开发的带有三阶微分方程的示例代码。

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()

大多数微分方程求解器要求您将高阶导数转换为单独的一阶导数方程。这很容易做到，因为您需要为每个附加订单（二阶和三阶导数）定义一个新状态为shown here。

【讨论】：