从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

我不知道.solve_ivp（）是如何工作的，但是如果你只想得到方程的字符串表示形式，我会做以下工作

import pandas as pd

df = pd.read_csv("test.csv")  # create DataFrame

def get_dt(some_variable):
    equation = ["d{}/dt =".format(some_variable)]  # store str representation of our eqn in a list
    variables = []
    for index, row in df.iterrows():  # iterate over all rows to check for a match
        if (row["From"] == some_variable) or (row["To"] == some_variable):  # if variable is found in row
            if len(variables):
                variables.append("- {}*{}".format(row["Rate"], row["From"]))  # add it to our equation
            else:
                variables.append("{}*{}".format(row["Rate"], row["From"]))  # the first number is positive for some reason
    equation.extend(variables)  # combine left and right hand sides into one list
    return " ".join(equation)  # combine our list or strings and return a nicely formatted single string

eqn = get_dt("Liver_1")
print(eqn)

输出

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

另外，如果你想一次得到所有的方程式

all_vars = set(df["From"])  # make a "list" of the From variables that doesn't contain duplicates
all_vars.update(set(df["To"]))  # add in all To variables still without duplicates
for var in all_vars:
    eqn = get_dt(var)
    print(eqn)

这给了我们输出

dT-bone-S/dt = 0.12474*Blood_5
dT-bone-V/dt = 0.013859999999999999*Blood_5
dUB-cont/dt = 0.0154*Blood_5 - 3.5*Blood_4 - 0.017329*Kidneys_1
dRC-cont/dt = 0.01155*Blood_5
dSI-cont/dt = 0.0009242*Liver_1
dOther_3/dt = 300.0*Blood - 0.099*Other_3 - 28.95*Blood_4
dLiver_2/dt = 0.045286*Liver_1
dOther_5/dt = 0.0231*Blood_5 - 0.0001266*Other_5
dLiver_1/dt = 0.462*Blood_5 - 0.0009242*Liver_1 - 0.045286*Liver_1
dOvaries/dt = 8.47e-05*Blood_5
dKidneys_1/dt = 0.0077*Blood_5 - 0.017329*Kidneys_1
dKidneys_2/dt = 0.000385*Blood_5 - 0.0001266*Kidneys_2
dBlood/dt = 300.0*Blood - 700.0*Blood
dOther_4/dt = 0.018511*Blood_5 - 0.0013859999999999999*Other_4
dBlood_5/dt = 700.0*Blood - 0.462*Blood_5 - 0.08778*Blood_5 - 0.00462*Blood_5 - 0.12474*Blood_5 - 0.013859999999999999*Blood_5 - 0.0154*Blood_5 - 0.0077*Blood_5 - 0.000385*Blood_5 - 0.01155*Blood_5 - 0.0002695*Blood_5 - 8.47e-05*Blood_5 - 0.018511*Blood_5 - 0.0231*Blood_5 - 0.099*Other_3 - 67.55*Blood_4
dBlood_4/dt = 3.5*Blood_4 - 67.55*Blood_4 - 28.95*Blood_4 - 0.0001266*Kidneys_2 - 0.0013859999999999999*Other_4 - 0.0001266*Other_5
dTestes/dt = 0.0002695*Blood_5
dC-bone-S/dt = 0.08778*Blood_5
dC-bone-V/dt = 0.00462*Blood_5

我认为这足以让你走，这样你就可以找到它，但它可能不完全是你想要的，因为奇怪的线条，如：

dBlood/dt = 300.0*Blood - 700.0*Blood

这对我来说似乎没什么意义

Answer 2

我编写了一个名为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))

一些注意事项：

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

• 我在这里打印输出，但您可以用不同的方式存储和处理输出

• 您还可以使用^{}使用您的标识符作为字典访问解决方案（通过dynvars）。这样，您就不需要直接处理动态变量的编号

• 基本的简化如

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

  将由SymEngine（象征性主干）在引擎盖下动态完成

• JiTCODE在引擎盖下使用scipy.integrate.ode或solve_ivp