How to optimize a system of differential equations with curve_fit and solve_ivp?

Question

I want to optimize parameters in a system of differential equations. For a minimal example, I have the differential equations:
dx1/dt=-k*x1 and
dx2/dt=k*x1.

Therefore, I use solve_ivp to integrate the ODE's and curve_fit to fit the parameters.

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)

This works, but if I increase t_end or make the differential equation more complicated, it always gets the same errors:

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

e.g. with t_end=1200 i get the error

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

or, with the slightly more complex system of ODE's and 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]

I get the below error:

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

I already figured out two options in solve_ivp which let me use higher t_end's and more complex systems of ODE's: Choosing method='DOP853' or method='LSODA' with min_step=1e-3 or 1e-2.

But in the end, it's not enough for the System I want to use (4-5 differential equations with 7-8 parameters to fit and t_end of 6000+).

So, Is there any general problem with my approach or any trick I missed?

Answer 1

I agree with Lutz Lehmann, the ValueError rise due to floating point overflow. To avoid it, set the maximum bounds on the parameters "k". For 64 bit numbers, maximum possible number is approximately 10^308, so for exponential solutions we need to take the upper bound < 308/t_end and lower bound > -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()