python - 从 CSV 文件编写微分方程(与 solve_ivp 一起使用)

Question

我正在尝试使用 solve_ivp 求解大型微分方程组。

    from scipy import integrate
    sol = integrate.solve_ivp(func_system, (0,100), initial_value_array, t_eval)

func_system 是一个微分方程组，必须从 CSV 文件中推断出来，其示例如下所示:

From,To,Rate
Blood,Other_3,3.0000E+02
Blood,Blood_5,7.0000E+02
Blood_5,Liver_1,4.6200E-01
Blood_5,C-bone-S,8.7780E-02
Blood_5,C-bone-V,4.6200E-03
Blood_5,T-bone-S,1.2474E-01
Blood_5,T-bone-V,1.3860E-02
Blood_5,UB-cont,1.5400E-02
Blood_5,Kidneys_1,7.7000E-03
Blood_5,Kidneys_2,3.8500E-04
Blood_5,RC-cont,1.1550E-02
Blood_5,Testes,2.6950E-04
Blood_5,Ovaries,8.4700E-05
Blood_5,Other_4,1.8511E-02
Blood_5,Other_5,2.3100E-02
Other_3,Blood_5,9.9000E-02
Blood_4,UB-cont,3.5000E+00
Blood_4,Blood_5,6.7550E+01
Blood_4,Other_3,2.8950E+01
Kidneys_1,UB-cont,1.7329E-02
Kidneys_2,Blood_4,1.2660E-04
Other_4,Blood_4,1.3860E-03
Other_5,Blood_4,1.2660E-04
Liver_1,SI-cont,9.2420E-04
Liver_1,Liver_2,4.5286E-02
[...]

例如，Liver_1 隔室的微分方程为:

dLiver_1/dt = 0.462*Blood_5 - 0.000924*Liver_1 - 0.045286*Liver_1

必须在 func_system 中编写成类似的东西

dLiver_1_dt = 0.462*Blood_5 - 0.000924*Liver_1 - 0.045286*Liver_1

其中 0.462 是从 Blood_5 到 Liver_1 的速率，0.000924 和 0.045286 是远离 Liver_1 的速率。

有没有一种方法可以在不实际编写的情况下创建所有这些方程式(我将总共有 150 多个)？ 我可以使用矩阵方法，但我还会有另一个需要添加的非线性微分方程组。

Answer 1

我写了一个名为 JiTCODE 的模块它允许您以符号形式指定微分方程。这对您来说很方便，因为它主要负责将方程式转换为函数。此外，在实际集成期间，您的函数评估将非常高效，因为它是硬编码在 C 中的。如果您有 150 个方程，这可能很重要。

对于你的具体问题，我会这样做:

from jitcode import jitcode,y
import numpy as np
from pandas import read_csv

# Dictionary describing the RHS of the ODE:
f = {}
def add_to_equation(dynvar,to_be_added):
    try:
        f[dynvar] += to_be_added
    except KeyError:
        f[dynvar]  = to_be_added

# Dictionary mapping identifiers to JiTCODE’s dynamic variables
dynvars = {}
def get_dynvar(name):
    try:
        return dynvars[name]
    except KeyError:
        new_dynvar = y(len(dynvars))
        dynvars[name] = new_dynvar
        return new_dynvar

# Parsing the file
for source_name,target_name,rate in read_csv("flows.csv").values:
    rate = float(rate)
    source = get_dynvar(source_name)
    target = get_dynvar(target_name)

    add_to_equation(target, rate*source)
    add_to_equation(source,-rate*source)

initial_condition = np.random.random(len(dynvars))

# Your test case:
Liver_1 = get_dynvar("Liver_1")
Blood_5 = get_dynvar("Blood_5")
assert f[Liver_1] == 0.462*Blood_5-0.0009242*Liver_1-0.045286*Liver_1

# Setup integrator
ODE = jitcode(f)
ODE.set_integrator("dopri5")
ODE.set_initial_value(initial_condition,0.0)

# Actual integration
times = np.linspace(0,100,50)
for time in times:
    print(ODE.integrate(time))

一些注意事项:

• 包含不同的非线性组件应该很简单。

• 我在这里打印输出，但您存储以不同方式进一步处理它。

• 您也可以使用 ODE.y_dict使用您的标识符(通过 dynvars)以字典形式访问解决方案。这样您就无需直接处理动态变量的编号。

• 基本的简化如

  -0.000924*Liver_1 - 0.045286*Liver_1 = -0.04621*Liver_1
  

  将由 SymEngine(符号主干)在引擎盖下动态完成。

• JiTCODE 在底层使用 scipy.integrate.ode 或 solve_ivp。