如何从文本文件中读取微分方程组，以使用scipy.odeint求解该系统? [英] How to read a system of differential equations from a text file to solve the system with scipy.odeint?
问题描述
我有一个大型(> 2000个方程)的ODE系统，我想用python scipy的odeint来解决.
I have a large (>2000 equations) system of ODE's that I want to solve with python scipy's odeint.
我要解决三个问题(也许我将不得不问三个不同的问题?). 为简单起见，我将在此处用玩具模型对其进行解释，但是请记住，我的系统很大. 假设我具有以下ODE系统:
I have three problems that I want to solve (maybe I will have to ask 3 different questions?). For simplicity, I will explain them here with a toy model, but please keep in mind that my system is large. Suppose I have the following system of ODE's:
dS/dt = -beta*S
dI/dt = beta*S - gamma*I
dR/dt = gamma*I
with beta = c p I
with beta = cpI
其中c，p和gamma是我要传递给odeint的参数.
where c, p and gamma are parameters that I want to pass to odeint.
odeint期望这样的文件:
odeint is expecting a file like this:
def myODEs(y, t, params):
c,p, gamma = params
beta = c*p
S = y[0]
I = y[1]
R = y[2]
dydt = [-beta*S*I,
beta*S*I - gamma*I,
- gamma*I]
return dydt
然后可以像这样传递给odeint:
that then can be passed to odeint like this:
myoutput = odeint(myODEs, [1000, 1, 0], np.linspace(0, 100, 50), args = ([c,p,gamma], ))
我在Mathematica中生成了一个文本文件，例如myOdes.txt，其中文件的每一行都对应于我的ODE系统的RHS，所以看起来像这样
I generated a text file in Mathematica, say myOdes.txt, where each line of the file corresponds to the RHS of my system of ODE's, so it looks like this
#myODEs.txt
-beta*S*I
beta*S*I - gamma*I
- gamma*I
我的文本文件看起来与odeint所期望的相似，但是我还没有到那里. 我有三个主要问题:
My text file looks similar to what odeint is expecting, but I am not quite there yet. I have three main problems:
- 如何传递我的文本文件，以便odeint理解这是系统的RHS?
- 如何以一种聪明的方式(即系统的方式)定义变量?由于它们有2000多个，因此我无法手动定义它们.理想情况下，我会在一个单独的文件中定义它们并阅读.
- 我如何也可以将参数(很多)作为文本文件传递?
我阅读了这个问题接近我的问题1和2并尝试将其复制(我直接为参数添加了值，这样我就不必担心上面的第3点):
I read this question that is close to my problems 1 and 2 and tried to copy it (I directly put values for the parameters so that I didn't have to worry about my point 3 above):
systemOfEquations = []
with open("myODEs.txt", "r") as fp :
for line in fp :
systemOfEquations.append(line)
def dX_dt(X, t):
vals = dict(S=X[0], I=X[1], R=X[2], t=t)
return [eq for eq in systemOfEquations]
out = odeint(dX_dt, [1000,1,0], np.linspace(0, 1, 5))
但是我得到了错误:
odepack.error: Result from function call is not a proper array of floats.
ValueError: could not convert string to float: -((12*0.01/1000)*I*S),
我将代码修改为:
systemOfEquations = []
with open("SIREquationsMathematica2.txt", "r") as fp :
for line in fp :
pattern = regex.compile(r'.+?\s+=\s+(.+?)$')
expressionString = regex.search(pattern, line)
systemOfEquations.append( sympy.sympify( expressionString) )
def dX_dt(X, t):
vals = dict(S=X[0], I=X[1], R=X[2], t=t)
return [eq for eq in systemOfEquations]
out = odeint(dX_dt, [1000,1,0], np.linspace(0, 100, 50), )
，这有效(我不太了解for循环的前两行在做什么).但是，我想更自动地定义变量，并且我仍然不知道如何使用此解决方案并在文本文件中传递参数.同样，如何在dX_dt函数中定义参数(取决于变量)?
and this works (I don't quite get what the first two lines of the for loop are doing). However, I would like to do the process of defining the variables more automatic, and I still don't know how to use this solution and pass parameters in a text file. Along the same lines, how can I define parameters (that will depend on the variables) inside the dX_dt function?
提前谢谢！
推荐答案
这不是一个完整的答案，而是一些观察/问题，但是对于评论来说太长了.
This isn't a full answer, but rather some observations/questions, but they are too long for comments.
dX_dt被odeint用1d数组y和元组t多次调用.您可以通过args参数提供t. y由odeint生成，并随每个步骤而变化. dX_dt应该简化以便快速运行.
dX_dt is called many times by odeint with a 1d array y and tuple t. You provide t via the args parameter. y is generated by odeint and varies with each step. dX_dt should be streamlined so it runs fast.
通常，像[eq for eq in systemOfEquations]这样的表达式可以简化为systemOfEquations. [eq for eq...]没有做任何有意义的事情.但是systemOfEquations可能有一些要求.
Usually an expresion like [eq for eq in systemOfEquations] can be simplified to systemOfEquations. [eq for eq...] doesn't do anything meaningful. But there may be something about systemOfEquations that requires it.
为您和我们的利益，建议您打印出systemOfEquations(对于这种3行小的情况).您正在使用sympy将文件中的字符串转换为方程式.我们需要看看它会产生什么.
I'd suggest you print out systemOfEquations (for this small 3 line case), both for your benefit and ours. You are using sympy to translated the strings from the file into equations. We need to see what it produces.
请注意，myODEs是函数，而不是文件.它可以从模块(当然是文件)中导入.
Note that myODEs is a function, not a file. It may be imported from a module, which of course is a file.
指向vals = dict(S=X[0], I=X[1], R=X[2], t=t)的目的是产生一个sympy表达式可以使用的字典.一个更直接(我认为更快)的dX_dt函数如下所示:
The point to vals = dict(S=X[0], I=X[1], R=X[2], t=t) is to produce a dictionary that the sympy expressions can work with. A more direct (and I think faster) dX_dt function would look like:
def myODEs(y, t, params):
c,p, gamma = params
beta = c*p
dydt = [-beta*y[0]*y[1],
beta*y[0]*y[1] - gamma*y[1],
- gamma*y[1]]
return dydt
我怀疑运行sympy生成的表达式的dX_dt会比这样的硬编码"表达式慢很多.
I suspect that the dX_dt that runs sympy generated expressions will be a lot slower than a 'hardcoded' one like this.
我要添加sympy标记，因为按照书面要求，这是将文本文件转换为odeint可以使用的功能的关键.
I'm going add sympy tag, because, as written, that is the key to translating your text file into a function that odeint can use.
我倾向于将方程可变性放在t参数中，而不是一列sympy表达式.
I'd be inclined to put the equation variability in the t parameters, rather a list of sympy expressions.
那是替换:
dydt = [-beta*y[0]*y[1],
beta*y[0]*y[1] - gamma*y[1],
- gamma*y[1]]
类似
arg12=np.array([-beta, beta, 0])
arg1 = np.array([0, -gamma, -gamma])
arg0 = np.array([0,0,0])
dydt = arg12*y[0]*y[1] + arg1*y[1] + arg0*y[0]
一旦这是对的，则argxx定义可以移到dX_dt之外，并通过args传递.现在dX_dt只是一个简单，快速的计算.
Once this is right, then the argxx definitions can be move outside dX_dt, and passed via args. Now dX_dt is just a simple, and fast, calculation.
整个sympy方法可能很好用，但恐怕在实践中会很慢.但是具有更多sympy经验的人可能会有其他见解.
This whole sympy approach may work fine, but I'm afraid that in practice it will be slow. But someone with more sympy experience may have other insights.
这篇关于如何从文本文件中读取微分方程组，以使用scipy.odeint求解该系统?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
