在 Python 中优化微分方程中的常数答案

【问题标题】：Optimize constants in differential equations in Python在 Python 中优化微分方程中的常数
【发布时间】：2013-04-19 20:25:51
【问题描述】：

好的，那么我将如何编写代码来优化微分方程中的常数 a 和 b，例如 dy/dt = a*y^2 + b，使用curve_fit？我将使用 odeint 来求解 ODE，然后使用 curve_fit 来优化 a 和 b。如果您能就这种情况提供意见，我将不胜感激！

【问题讨论】：

标签： python scipy curve-fitting ode odeint

【解决方案1】：

查看ODEs with Sympy 可能会为您提供更好的服务。 Scipy/Numpy 基本上是数值包，并没有真正设置为进行代数/符号运算。

【讨论】：

谢谢，这很有趣，可能会有所帮助，就像我想要对微分方程做的事情一样，一旦被 odeint 解决，就是使用curve_fit优化一个或两个常数。你认为我必须使用 sympy 来执行此操作还是有其他方法？
更合适的 SymPy 链接是docs.sympy.org/dev/modules/solvers/ode.html
@user2199360：虽然 Scipy 没有针对符号运算进行设置，但 Sympy 也没有针对优化等数值运算进行设置。如果您使用curvefit，您正在拟合的曲线类型可能会为您提供可用于创建符号函数的数值结果的解释（例如，如果参数是多项式系数，那么您可以使用这些来写下符号形式的多项式）。也许如果你编辑你的问题来描述你试图完成的整体任务，人们可以提供更多有用的输入。
这个问题实际上不需要 ODE 的解析解（如果 RHS 更复杂，Sympy 不会找到）。

【解决方案2】：

你绝对可以这样做：

import numpy as np
from scipy.integrate import odeint
from scipy.optimize import curve_fit

def f(y, t, a, b):
    return a*y**2 + b

def y(t, a, b, y0):
    """
    Solution to the ODE y'(t) = f(t,y,a,b) with initial condition y(0) = y0
    """
    y = odeint(f, y0, t, args=(a, b))
    return y.ravel()

# Some random data to fit
data_t = np.sort(np.random.rand(200) * 10)
data_y = data_t**2 + np.random.rand(200)*10

popt, cov = curve_fit(y, data_t, data_y, [-1.2, 0.1, 0])
a_opt, b_opt, y0_opt = popt

print("a = %g" % a_opt)
print("b = %g" % b_opt)
print("y0 = %g" % y0_opt)

import matplotlib.pyplot as plt
t = np.linspace(0, 10, 2000)
plt.plot(data_t, data_y, '.',
         t, y(t, a_opt, b_opt, y0_opt), '-')
plt.gcf().set_size_inches(6, 4)
plt.savefig('out.png', dpi=96)
plt.show()

【讨论】：

【解决方案3】：

为了专门解决这类问题，我决定编写一个统一sympy 和scipy 的包装包。它被称为symfit。适合您的 ODE 将如下所示：

tdata = np.array([10, 26, 44, 70, 120])
ydata = 10e-4 * np.array([44, 34, 27, 20, 14])
y, t = variables('y, t')
a, b = parameters('a, b')

model_dict = {
    D(y, t): a*y^2 + b
}

ode_model = ODEModel(model_dict, initial={t: 0.0, y: 0.0})

fit = Fit(ode_model, t=tdata, y=ydata)
fit_result = fit.execute()

从定义为 dict 的方式可以看出，适合（一阶）ODE 系统是没有问题的。查看docs 了解更多信息！

【讨论】：