使用 scipy.integrate 一段时间后,我需要更多功能,如分岔分析或参数估计。这就是为什么我对使用 PyDSTool 感兴趣,但是从文档中我无法弄清楚如何使用 ModelSpec 以及这是否真的会引导我找到解决方案。
这是我正在尝试做的一个玩具示例:我有一个包含两个节点的网络,两个节点都具有相同的 (SIR) 动态,由两个 ODE 描述,但初始条件不同。方程通过 Epsilon 在节点之间耦合(见下面的公式)。
公式作为图片以便更好地阅读,“n”和“m”是指数,而不是指数~>
http://image.noelshack.com/fichiers/2014/28/1404918182-odes.png
(遗憾的是,无法使用堆栈上传)
在两个节点的情况下,我的代码(使用 PyDSTool)如下所示:
#multiple SIR metapopulations
#parameter and initial condition definition; a dict is a must
import PyDSTool as pdt
params={'alpha': 0.7, 'beta':0.1, 'epsilon1':0.5,'epsilon2':0.5}
ini={'s1':0.99,'s2':1,'i1':0.01,'i2':0.00}
DSargs=pdt.args(name='SIRtest_multi',
ics=ini,
pars=params,
tdata=[0,20],
#the for-macro generates formulas for s1,s2 and i1,i2;
#sum works similar but sums over the expressions in it
varspecs={'s[o]':'for(o,1,2,-alpha*s[o]*sum(k,1,2,epsilon[k]*i[k]))',
'i[l]':'for(l,1,2,alpha*s[l]*sum(m,1,2,epsilon[m]*i[m]))'})
#generator
DS = pdt.Generator.Vode_ODEsystem(DSargs)
#computation, a trajectory object is generated
trj=DS.compute('test')
#extraction of the points for plotting
pts=trj.sample()
#plotting; pylab is imported along with PyDSTool as plt
pdt.plt.plot(pts['t'],pts['s1'],label='s1')
pdt.plt.plot(pts['t'],pts['i1'],label='i1')
pdt.plt.plot(pts['t'],pts['s2'],label='s2')
pdt.plt.plot(pts['t'],pts['i2'],label='i2')
pdt.plt.legend()
pdt.plt.xlabel('t')
pdt.plt.show()
但是在我最初的问题中,每个节点有 1000 多个节点和 5 个 ODE,每个节点都耦合到不同数量的其他节点,并且所有节点的 epsilon 值并不相等。所以修改这个语法并没有让我接近解决方案。
我实际上在想的是一种为每个节点构建单独的子模型/求解器(?)的方法,具有自己的参数(epsilons,因为它们对于每个节点都不同)。然后将它们相互链接。这就是我不知道在 PyDSTool 中是否可行以及它是否是处理此类问题的方法的地方。
我查看了 PyDSTool 的示例和文档,但不知道该怎么做,因此非常感谢您的帮助!如果我尝试做事的方式不正统或愚蠢,欢迎您提出如何更有效地做事的建议。 (实际上哪种方法更有效/更快/更好地解决此类问题:将其分割为许多小的(仍未解耦的)模型/求解器,还是一次包含所有 ODE 的模型/求解器?)
(我既不是数学家也不是程序员,但愿意学习,所以请耐心等待!)
最佳答案
解决方案绝对不是建立单独的仿真模型。那是行不通的,因为子模型之间会持续耦合这么多变量。您绝对必须将所有 ODE 放在一个地方。
听起来您需要的解决方案是使用 ModelSpec 对象构造。这些使您可以从符号片段中分层构建子模型定义。它们可以有自己的“epsilon”参数等。完成后声明所有部分,让 PyDSTool 为您制作包含 ODE 定义的最终字符串。我建议您查看以下教程示例:
http://www.ni.gsu.edu/~rclewley/PyDSTool/Tutorial/Tutorial_compneuro.html
以及提供的示例:ModelSpec_test.py、MultiCompartments.py。但是,请记住,您仍然必须拥有参数和耦合数据的来源(即从文件加载的大矩阵或字典)才能使构建模型的过程自动化,否则您仍然需要编写它全部手动。
您必须为想要拥有的组件构建一些类。您还可以创建一个工厂函数(比较neuralcomp.py 工具箱中的'makeSoma'),它将获取您的所有子组件并根据从每个声明的组件中总结出来的内容来创建一个ODE。最后,您可以通过参数在层次结构中的位置来引用参数。一个可能是“s1.epsilon”,而另一个可能是“i4.epsilon”。
不幸的是,要有效地构建这样的模型,您必须学习进行一些更复杂的编程!因此,首先要了解本教程中的所有步骤。一旦您开始并有具体问题,您可以通过 SourceForge 支持讨论或电子邮件直接给我发送电子邮件。
关于python - 使用 PyDSTool 解决网络上的 ODE,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/24657820/