2

Okay so how would i approach to writing a code to optimize the constants a and b in a differential equation, like dy/dt = a*y^2 + b, using curve_fit? I would be using odeint to solve the ODE and then curve_fit to optimize a and b. If you could please provide input on this situation i would greatly appreciate it!

user2199360
  • 31
  • 2
  • 4

3 Answers3

4

You might be better served by looking at ODEs with Sympy. Scipy/Numpy are fundamentally numerical packages and aren't really set up to do algebraic/symbolic operations.

BrenBarn
  • 242,874
  • 37
  • 412
  • 384
  • Thanks thats interesting and may help, like what i want to do with the differential equation, once solved by odeint, is to optimize one or two constants using curve_fit. do you think i will have to use sympy to do this or is there another way? – user2199360 Apr 25 '13 at 19:41
  • A more appropriate SymPy link is http://docs.sympy.org/dev/modules/solvers/ode.html – Warren Weckesser Apr 25 '13 at 19:46
  • @user2199360: While Scipy isn't set up for symbolic operations, Sympy isn't set up for numerical operations like optimizing. If you use curvefit, the type of curve you're fitting may provide you with an interpretation of the numerical results that you can use to create a symbolic function (e.g., if the parameters are polynomial coefficients, then you can use those to write down the polynomial in symbolic form). Perhaps if you edit your question to describe the overall task you are trying to accomplish, people can give more useful input. – BrenBarn Apr 25 '13 at 19:47
  • 2
    This problem actually doesn't require an analytical solution of the ODE (which Sympy would not find if the RHS was more complicated). – pv. Apr 26 '13 at 17:27
4

You definitely can do this:

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()

pv.
  • 33,875
  • 8
  • 55
  • 49
0

To address specifically this type of problem, I decided to write a wrapper package which unifies sympy and scipy. It's called symfit. Fitting to your ODE would then look like this:

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()

As you can see from the way it is defined as a dict, fitting to systems of (first order) ODEs is no problem. Check out the docs for more!

tBuLi
  • 2,295
  • 2
  • 16
  • 16