Учет данных для численного решения оды в python

Question

У меня есть система из двух ОДЭ первого порядка, которые являются нелинейными и, следовательно, трудно решить аналитически в замкнутой форме. Я хочу подогнать численное решение этой системы ODE к набору данных. Мой набор данных относится только к одной из двух переменных, которые являются частью системы ODE. Как мне это сделать? This не помогло, потому что там есть только одна переменная.

Мой код, который в настоящее время ведет к ошибке:

import numpy as np 
from scipy.integrate import odeint 
from scipy.optimize import curve_fit 

def f(y, t, a, b, g): 
    S, I = y # S, I are supposed to be my variables 
    Sdot = -a * S * I 
    Idot = (a - b) * S * I + (b - g - b * I) * I 
    dydt = [Sdot, Idot] 
    return dydt 

def y(t, a, b, g, y0): 
    y = odeint(f, y0, t, args=(a, b, g)) 
    return y.ravel() 

I_data =[] # I have data only for I, not for S 
file = open('./ratings_showdown.csv') 
for e_raw in file.read().split('\r\n'): 
    try: 
     e=float(e_raw); I_data.append(e) 
    except ValueError: 
     continue 

data_t = range(len(I_data)) 
popt, cov = curve_fit(y, data_t, I_data, [.05, 0.02, 0.01, [0.99,0.01]]) 
#want to fit I part of solution to data for variable I 
#ERROR here, ValueError: setting an array element with a sequence 
a_opt, b_opt, g_opt, y0_opt = popt 

print("a = %g" % a_opt) 
print("b = %g" % b_opt) 
print("g = %g" % g_opt) 
print("y0 = %g" % y0_opt) 

import matplotlib.pyplot as plt 
t = np.linspace(0, len(data_y), 2000) 
plt.plot(data_t, data_y, '.', 
     t, y(t, a_opt, b_opt, g_opt, y0_opt), '-') 
plt.gcf().set_size_inches(6, 4) 
#plt.savefig('out.png', dpi=96) #to save the fit result 
plt.show() 

Answer 1

Так я понял, проблема. Функция curve_fit(), по-видимому, возвращает список, поскольку это второе возвращаемое значение. Итак, вместо передачи начальных условий в виде списка [0.99,0.01] я передал их отдельно как 0.99 и 0.01.

Answer 2

Этот тип установки ODE становится намного проще в symfit, который я написал специально как удобную для пользователя упаковку до scipy. Я думаю, что это будет очень полезно для вашей ситуации, потому что уменьшенное количество котельного кода упрощает многое.

Из документов и применяется примерно к вашей проблеме:

from symfit import variables, parameters, Fit, D, ODEModel 

S, I, t = variables('S, I, t') 
a, b, g = parameters('a, b, g') 

model_dict = { 
    D(S, t): -a * S * I, 
    D(I, t): (a - b) * S * I + (b - g - b * I) * I 
} 

ode_model = ODEModel(model_dict, initial={t: 0.0, S: 0.99, I: 0.01}) 

fit = Fit(ode_model, t=tdata, I=I_data, S=None) 
fit_result = fit.execute() 

Отъезд docs для более :)