Оптимизация констант в дифференциальных уравнениях в Python

Итак, как мне подойти к написанию кода для оптимизации констант a и b в дифференциальном уравнении, например dy/dt = a*y^2 + b, с помощью curve_fit? Я бы использовал odeint для решения ODE, а затем curve_fit для оптимизации a и b. Если бы вы могли, пожалуйста, внести свой вклад в эту ситуацию, я был бы очень признателен!


python scipy curve-fitting ode odeint
person user2199360    schedule 25.04.2013    source источник

Ответы (3)


arrow_upward
4
arrow_downward

Возможно, вам будет лучше взглянуть на ODE с Sympy. Scipy/Numpy - это в основном числовые пакеты, которые на самом деле не предназначены для выполнения алгебраических/символических операций.

person BrenBarn    schedule 25.04.2013
comment
Спасибо, это интересно и может помочь, например, то, что я хочу сделать с дифференциальным уравнением, однажды решенным с помощью odeint, - это оптимизировать одну или две константы с помощью curve_fit. Как вы думаете, мне придется использовать sympy, чтобы сделать это, или есть другой способ? - person user2199360; 25.04.2013
comment
Более подходящей ссылкой SymPy является docs.sympy.org/dev/modules/solvers/ ode.html - person Warren Weckesser; 25.04.2013
comment
@ user2199360: Хотя Scipy не настроен для символьных операций, Sympy не настроен для числовых операций, таких как оптимизация. Если вы используете подгонку кривой, тип кривой, которую вы подгоняете, может дать вам интерпретацию числовых результатов, которые вы можете использовать для создания символьной функции (например, если параметры являются полиномиальными коэффициентами, вы можете использовать их для записи многочлен в символьной форме). Возможно, если вы отредактируете свой вопрос, чтобы описать общую задачу, которую вы пытаетесь выполнить, люди смогут внести более полезный вклад. - person BrenBarn; 25.04.2013
comment
На самом деле эта задача не требует аналитического решения ОДУ (которого Sympy не нашел бы, если бы RHS была более сложной). - person pv.; 26.04.2013

arrow_upward
4
arrow_downward

Вы определенно можете сделать это:

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

person pv.    schedule 26.04.2013

arrow_upward
0
arrow_downward

Чтобы решить именно этот тип проблемы, я решил написать пакет-оболочку, который объединяет sympy и scipy. Он называется symfit. Подгонка к вашему ODE будет выглядеть так:

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

Как вы можете видеть из того, как он определяется как dict, подгонка к системам ОДУ (первого порядка) не проблема. Ознакомьтесь с документами, чтобы узнать больше!

person tBuLi    schedule 23.10.2016