Python solve_ivp 与 R lsoda

Question

我正在尝试将 R 代码转换为 Python。R 代码使用lsoda函数，它是 FORTRANDOE求解器的包装器。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

        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.9997727 -0.2287976

0.9995454 -0.22868

0.9995455 -0.2286801

0.9901534 -0.2238221

0.9901523 -0.2238215

0.9809609 -0.2190673

0.9809587 -0.2190662

0.9719626 -0.214413

0.9719604 -0.2144119

0.9493722 -0.2027284

0.9493668 -0.2027255

0.927996 -0.1916717

0.9280039 -0.1916758

0.9078033 -0.1812272

0.9078049 -0.181228

0.8887056 -0.1713491

0.8887071 -0.1713499

Answer 1

正如@astoeriko 所指出的，默认相对容差 (&nbsp;rtol) 在 scipy 中为 1e-3solve_ivp而在 R 中为 1e-6&nbsp;lsoda。