Question

This is a specialized application of curve_fit , in which the curve to be fitted is defined implcitly by an ordinary differentail equation

\[\frac{dx}{dt} = -kx\]

and we want to use observed data to estiamte the parameters \(k\) and the initial value \(x_0\). Of course this can be explicitly solved but the same approach can be used to find multiple paraemters for \(n\)-dimensional systems of ODEs.

A more elaborate example for fitting a system of ODEs to model the zombie apocalypse

from scipy.integrate import odeint

def f(x, t, k): “”“Simple exponential decay.”“” return -k*x

def x(t, k, x0): “”” Solution to the ODE x’(t) = f(t,x,k) with initial condition x(0) = x0 “”” x = odeint(f, x0, t, args=(k,)) return x.ravel()

# True parameter values
x0_ = 10
k_ = 0.1*np.pi

# Some random data genererated from closed form soltuion plus Gaussian noise
ts = np.sort(np.random.uniform(0, 10, 200))
xs = x0_*np.exp(-k_*ts) + np.random.normal(0,0.1,200)

popt, cov = curve_fit(x, ts, xs)
k_opt, x0_opt = popt

print("k = %g" % k_opt)
print("x0 = %g" % x0_opt)

k = 0.314062
x0 = 9.754

import matplotlib.pyplot as plt
t = np.linspace(0, 10, 100)
plt.plot(ts, xs, '.', t, x(t, k_opt, x0_opt), '-');