python - 如何使用 curve_fit 和 solve_ivp 优化微分方程组？

Question

我想优化微分方程组中的参数。举一个最小的例子，我有微分方程:
dx1/dt=-k*x1 和
dx2/dt=k*x1。

因此，我使用 solve_ivp 对 ODE 进行积分，并使用 curve_fit 来拟合参数。

t_end=5
t_eval=np.linspace(0,t_end,t_end+1) #define timespan of 

def fun(t,x,k): #System of differential equations with parameter k to optimize
       return [-k*x[0],k*x[0]]

def fun2(t,k): #Solving differential equations and return x2(t)
       res=solve_ivp(fun, [0,t_end], [1,0],t_eval=t_eval,args=(k,))
       return res.y[1]

z=(1-np.exp(-0.2*t_eval))*np.ones(len(t_eval)) #creating an array with the analytic solution for the above statet differential equations with k=0.2
popt, pcov = curve_fit(fun2, t_eval, z)

这行得通，但如果我增加 t_end 或使微分方程更复杂，它总是会出现相同的错误:

ValueError: operands could not be broadcasted together with shapes (xxx,) (xxx,)

例如t_end=1200 我得到错误

ValueError: operands could not be broadcasted together with shapes (857,) (1201,)

或者，使用 ODE 和 t_end = 10 的稍微复杂的系统:

def fun(t,x,k1,k2): #System of differential equations
   return [-k1*x[0],k1*x[0]-k2*x[1],k2*x[1]]

def fun2(t,k1,k2): #Solving differential equations and return x3(t)
   res=solve_ivp(fun, [0,t_end], [1,0,0],t_eval=t_eval,args=(k1,k2))
   return res.y[2]

我收到以下错误:

ValueError: operands could not be broadcasted together with shapes (8,) (11,)

我已经在 solve_ivp 中找到了两个选项，它们让我可以使用更高的 t_end 和更复杂的 ODE 系统:选择 method='DOP853' 或 method='LSODA' 与 min_step=1e-3 或 1e-2。

但最终，这对于我想要使用的系统来说还不够(4-5 个具有 7-8 个参数的微分方程可以拟合并且 t_end 为 6000+)。

那么，我的方法是否存在任何一般性问题或我错过了什么技巧？

Answer 1

我同意 Lutz Lehmann 的观点，浮点溢出导致 ValueError 上升。为避免这种情况，请设置参数“k”的最大界限。对于 64 位数字，最大可能数约为 10^308，因此对于指数解，我们需要取上界 < 308/t_end 和下界 > -307/t_end:

t_end=6200
n=100
t_eval=np.linspace(0,t_end,n)
def fun(t,x,k1,k2):
   return [-k1*x[0],k1*x[0]-k2*x[1],k2*x[1]]
def fun2(t,k1,k2):
   res=solve_ivp(fun, [0,t_end], [1,0,0],t_eval=t_eval,args=(k1,k2))
   return res.y[2]

z=1-0.5*np.exp(-0.2*t_eval)-0.5*np.exp(-0.4*t_eval)

popt, pcov = curve_fit(fun2, t_eval, z,bounds=[(-307/t_end,-307/t_end),(308/t_end,308/t_end)])
res=solve_ivp(fun, [0,t_end], [1,0,0],t_eval=t_eval,args=(popt[0],popt[1]))

plt.plot(t_eval,res.y[2],'r', label="optimized data")
plt.plot(t_eval,z,'k--',label="data")
plt.legend()
plt.show()