Python solve_ivp vs R lsoda - PullRequest
Новая
       27

Python solve_ivp vs R lsoda

0 голосов
/ 22 февраля 2019

Я пытаюсь преобразовать код R в Python.Код R использует функцию lsoda, которая является оберткой для решателя FORTRAN DOE.Аналогом Python является solve_ivp, который является оболочкой для FORTRAN ODEPACK.Я использую method='LSODA' в Python, который должен быть эквивалентом того, что использует R.Тем не менее, мои результаты отличаются с ошибкой до 1%.В моем коде нет ничего случайного, поэтому я считаю, что смогу полностью воспроизвести результаты.

Любая идея?!

Это часть кода R (код до этого простовычисление значений для параметров: 

val = c("A1" = 1, "A2" = 1, "A3" = 1, "A4" = 1, "A5" = 1, "A6" = 1, "A7" = 1)

hamberg_ode <- function(t,val,p) { 
    dA1 = p["ktr1"]*(1 - ((p["E_MAX"] * p["C_s_gamma"])/(p["EC_50_gamma"] + p["C_s_gamma"]))) - p["ktr1"]*val["A1"]
    dA2 = p["ktr1"]*val["A1"] - p["ktr1"]*val["A2"]
    dA3 = p["ktr1"]*val["A2"] - p["ktr1"]*val["A3"]
    dA4 = p["ktr1"]*val["A3"] - p["ktr1"]*val["A4"]
    dA5 = p["ktr1"]*val["A4"] - p["ktr1"]*val["A5"]
    dA6 = p["ktr1"]*val["A5"] - p["ktr1"]*val["A6"]
    dA7 = p["ktr2"]*(1 - ((p["E_MAX"] * p["C_s_gamma"])/(p["EC_50_gamma"] + p["C_s_gamma"]))) - p["ktr2"]*val["A7"]
    cat(val["A1"], dA1, '\n')
    list(c(dA1, dA2, dA3, dA4, dA5, dA6, dA7))
}

out = lsoda(val, times, hamberg_ode, p)

Код Python: 

val = [1]*7

class hamberg_ode:
    def __init__(self, p):
        self.p = p

    def f(self, t, val, p=None):
        if p is None:
            p = self.p
        dA1=p["ktr1"]*(1 - ((p["E_MAX"] * p["C_s_gamma"]) /
                        (p["EC_50_gamma"] + p["C_s_gamma"]))) - p["ktr1"]*val[0]
        dA2=p["ktr1"]*val[0] - p["ktr1"]*val[1]
        dA3=p["ktr1"]*val[1] - p["ktr1"]*val[2]
        dA4=p["ktr1"]*val[2] - p["ktr1"]*val[3]
        dA5=p["ktr1"]*val[3] - p["ktr1"]*val[4]
        dA6=p["ktr1"]*val[4] - p["ktr1"]*val[5]
        dA7=p["ktr2"]*(1 - ((p["E_MAX"] * p["C_s_gamma"]) /
                        (p["EC_50_gamma"] + p["C_s_gamma"]))) - p["ktr2"]*val[6]
        print(val[0], dA1)

        return (dA1, dA2, dA3, dA4, dA5, dA6, dA7)


h_function = hamberg_ode(p).f
out = solve_ivp(h_function, (0, maxTime), val, t_eval=times, method='LSODA')

В качестве примера расхождения чисел ниже приведены несколько первых значений A1 и dA1 для двух кодов: R

1 -0,2289151

1 -0,2289151

0,9997726 -0,2287975

0,999927 -0,2287976

0,9995454-0,22868

0,9995455 -0,2286801

0,9901534 -0,2238221

0,9901523 -0,2238215

0,9809609 -0,2190673

0,21 809587

0,9719626 -0,214413

0,9719604 -0,2144119

0,9493722 -0,2027284

0,9493668 -0,2025255

0,927996 -0,19009 * 10,192839 * 10,1162839 * 10,10 028 * 10480,1916758

0,9078033 -0,1812272

0,9078049 -0,181228

0,8887056-0.1713491

0.8887071 -0.1713499

Python

1.0 -0.22891514470392998

0.9868338969217406 -0.22210509138758867 1066+0,22230768534978135

1069 * +0,9744278526945864 -0,2156881719597506 1071 * +0,9748085754105683 -0,2158850975024998 1073 * +0,9069550726140441 -0,18078845812498728 1075 * +0,906362742770375 -0,18048208061964116 1077 * +0,8502494750308627 -0,15145797661644517 1079 * 0,8491489787959607 -0,15088875442597866

0,8022897024620746 -0,1266511977015548

0,8013657199203642 -0,1261732756972218

0,7405758555885625 -0,094730242422152

0,7400188154524862 -0,09444211821383663

0,7145516960742005 -0,08126947025955095

0,7148846597052525 -0,08144169282733643

1 Ответ

0 голосов
/ 23 февраля 2019

Как указывает @astoeriko, относительный допуск по умолчанию (rtol) составляет 1e-3 для solve_ivp Сципи, но 1e-6 для R lsoda.

...