我有一个同意写的ODE系统:
from sympy.parsing.sympy_parser import parse_expr
xs = symbols('x1 x2')
ks = symbols('k1 k2')
strs = ['-k1 * x1**2 + k2 * x2', 'k1 * x1**2 - k2 * x2']
syms = [parse_expr(item) for item in strs]
我想将其转换为向量值函数,接受x值的1D numpy数组,k值的1D numpy数组,返回在这些点处评估的1D numpy数组 . 签名看起来像这样:
import numpy as np
x = np.array([3.5, 1.5])
k = np.array([4, 2])
xdot = my_odes(x, k)
我想要这样的东西的原因是将此功能赋予 scipy.integrate.odeint
,因此它需要快速 .
Attempt 1: subs
当然,我可以在 subs
周围写一个包装器:
def my_odes(x, k):
all_dict = dict(zip(xs, x))
all_dict.update(dict(zip(ks, k)))
return np.array([sym.subs(all_dict) for sym in syms])
但这是超级慢的,特别是对于我的真实系统,它更大,并且运行了很多次 . 我需要将此操作编译为C代码 .
Attempt 2: theano
我可以接近sympy's integration with theano:
from sympy.printing.theanocode import theano_function
f = theano_function(xs + ks, syms)
def my_odes(x, k):
return np.array(f(*np.concatenate([x,k]))))
这会编译每个表达式,但输入和输出的所有打包和解包都会减慢它的速度 . theano_function
返回的函数接受numpy数组作为参数,但每个符号需要一个数组,而不是每个符号一个元素 . 这与 functify
和 ufunctify
的行为相同 . 我不需要广播行为;我需要它将数组的每个元素解释为一个不同的符号 .
Attempt 3: DeferredVector
如果我使用 DeferredVector
我可以创建一个接受numpy数组的函数,但我无法将其编译为C代码或返回一个numpy数组而不自行打包 .
import numpy as np
import sympy as sp
from sympy import DeferredVector
x = sp.DeferredVector('x')
k = sp.DeferredVector('k')
deferred_syms = [s.subs({'x1':x[0], 'x2':x[1], 'k1':k[0], 'k2':k[1]}) for s in syms]
f = [lambdify([x,k], s) for s in deferred_syms]
def my_odes(x, k):
return np.array([f_i(x, k) for f_i in f])
使用 DeferredVector
我不需要解压缩输入,但我仍然需要打包输出 . 此外,我可以使用 lambdify
,但 ufuncify
和 theano_function
版本会消失,因此不会生成快速C代码 .
from sympy.utilities.autowrap import ufuncify
f = [ufuncify([x,k], s) for s in deferred_syms] # error
from sympy.printing.theanocode import theano_function
f = theano_function([x,k], deferred_syms) # error
2 回答
您可以使用sympy功能lambdify . 例如,
运行脚本时会打印
True
.速度要快得多(注意时序结果单位的变化):
我写了a module named JiTCODE,它适合你的问题 . 它需要符号表达式,将它们转换为C代码,围绕它包装Python扩展,编译它,并加载它以用于
scipy.integrate.ode
或scipy.integrate.solve_ivp
.您的示例如下所示:
然后,您可以使用
ODE
,就像scipy.integrate.ode
的实例一样 .虽然您不需要为您的应用程序使用它,但您也可以提取并使用已编译的函数:
请注意,与您的规范不同,
k
不作为NumPy数组传递 . 对于大多数ODE应用程序,这应该不相关,因为对控制参数进行硬编码会更有效 .最后,请注意,对于这个小例子,由于
scipy.integrate.ode
或scipy.integrate.solve_ivp
(也见SciPy Issue #8257或this answer of mine)的开销,您可能无法获得最佳性能 . 对于大的微分方程(如你所知),这种开销变得无关紧要 .