How to determine unknown parameters of a differential equation based on the best fit to a data set in Python?

Question

I am trying to fit different differential equations to a given data set with python. For this reason, I use the scipy package, respectively the solve_ivp function. This works fine for me, as long as I have a rough estimate of the parameters (b= 0.005) included in the differential equations, e.g:

import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
import numpy as np


def f(x, y, b):
    dydx= [-b[0] * y[0]]
    return dydx

xspan= np.linspace(1, 500, 25)
yinit= [5]
b= [0.005]


sol= solve_ivp(lambda x, y: f(x, y, b),
               [xspan[0], xspan[-1]], yinit, t_eval= xspan)

print(sol)
print("\n")

print(sol.t)
print(sol.y)

plt.plot(sol.t, sol.y[0], "b--")

However, what I like to achieve is, that the parameter b (or more parameters) is/are determined "automatically" based on the best fit of the solved differential equation to a given data set (x and y). Is there a way this can be done, for example by combining this example with the curve_fit function of scipy and how would this look?

Thank you in advance!

Answer 1

Yes, what you think about should work, it should be easy to plug together. You want to call

popt, pcov = scipy.optimize.curve_fit(curve, xdata, ydata, p0=[b0])
b = popt[0]

where you now have to define a function curve(x,*p) that transforms any list of point into a list of values according to the only parameter b.

def curve(x,b):
    res = solve_ivp(odefun, [1,500], [5], t_eval=x, args = [b])
    return res.y[0]

Add optional arguments for error tolerances as necessary.

---

To make this more realistic, make also the initial point a parameter. Then it also becomes more obvious where a list is expected and where single arguments. To get a proper fitting task add some random noise to the test data. Also make the fall to zero not so fast, so that the final plot still looks somewhat interesting.

from scipy.integrate import solve_ivp
from scipy.optimize import curve_fit

xmin,xmax = 1,500

def f(t, y, b):
    dydt= -b * y
    return dydt

def curve(t, b, y0):
    sol= solve_ivp(lambda t, y: f(t, y, b),
               [xmin, xmax], [y0], t_eval= t)
    return sol.y[0]


xdata = np.linspace(xmin, xmax, 25)
ydata = np.exp(-0.02*xdata)+0.02*np.random.randn(*xdata.shape)

y0 = 5
b= 0.005
p0 = [b,y0]

popt, pcov = curve_fit(curve, xdata, ydata, p0=p0)
b, y0 = popt

print(f"b={b}, y0 = {y0}")

This returns

b=0.019975693539459473, y0 = 0.9757709108115179

Now plot the test data against the fitted curve