python complex_ode 传递矩阵值参数

Question

我在使用 python 的 complex_ode 求解器时遇到了一些问题。

我正在尝试解决以下等式：

dy/dt = -iAy - icos(Omegat)By

其中 A 和 B 是 NxN 数组，未知 y 是 Nx1 数组，i 是虚数单位，Omega 是参数。

这是我的代码：

import numpy as np
from scipy.integrate import ode,complex_ode


N = 3 #linear matrix dim
Omega =  1.0 #parameter

# define symmetric matrices A and B
A = np.random.ranf((N,N))
A = (A + A.T)/2.0

B = np.random.ranf((N,N))
B = (B + B.T)/2.0

# define RHS of ODE
def f(t,y,Omega,A,B):
     return -1j*A.dot(y)-1j*np.cos(Omega*t)*B.dot(y)

# define list of parameter
params=[Omega,A,B]

# choose solver: need complex_ode for this ODE
#solver = ode(f)
solver = complex_ode(f)
solver.set_f_params(*params)
solver.set_integrator("dop853")
# set initial value
v0 = np.zeros((N,),dtype=np.float64)
v0[0] = 1.0

# check that the function f works properly
print f(0,v0,Omega,A,B)

# solve-check the ODE
solver.set_initial_value(v0)
solver.integrate(10.0)

print solver.successful()

运行此脚本会产生错误

capi_return is NULL
Call-back cb_fcn_in___user__routines failed.
Traceback (most recent call last):
File "ode_test.py", line 37, in <module>
   solver.integrate(10.0)
File "/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python/scipy/integrate/_ode.py", line 515, in integrate
y = ode.integrate(self, t, step, relax)
File "/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python/scipy/integrate/_ode.py", line 388, in integrate
self.f_params, self.jac_params)
File "/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python/scipy/integrate/_ode.py", line 946, in run
tuple(self.call_args) + (f_params,)))
File "/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python/scipy/integrate/_ode.py", line 472, in _wrap
f = self.cf(*((t, y[::2] + 1j * y[1::2]) + f_args))
TypeError: f() takes exactly 5 arguments (2 given)

如果我改为使用求解器 = ode(f)，即。实值求解器，它运行良好。除了它不能解决我想要的复值 ODE :(

然后我尝试通过使矩阵 A 和 B 成为全局变量来减少参数的数量。这样，函数 f 接受的唯一参数是 Omega。错误变为

capi_return is NULL
Call-back cb_fcn_in___user__routines failed.
Traceback (most recent call last):
File "ode_test.py", line 37, in <module>
solver.integrate(10.0)
File "/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python/scipy/integrate/_ode.py", line 515, in integrate
y = ode.integrate(self, t, step, relax)
File "/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python/scipy/integrate/_ode.py", line 388, in integrate
self.f_params, self.jac_params)
File "/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python/scipy/integrate/_ode.py", line 946, in run
tuple(self.call_args) + (f_params,)))
File "/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python/scipy/integrate/_ode.py", line 472, in _wrap
f = self.cf(*((t, y[::2] + 1j * y[1::2]) + f_args))
TypeError: 'float' object has no attribute '__getitem__'

我发现 float 指的是参数 Omega [通过尝试整数]。同样，“ode”也适用于这种情况。

最后，我尝试了相同的复值方程，但现在 A 和 B 只是数字。我尝试将它们都作为参数传递，即 params = [Omega,A,B]，以及在这种情况下将它们设为全局变量 params=[Omega]。错误是

TypeError: 'float' object has no attribute '__getitem__'

error - 完整的错误与上面相同。再一次，对于实值“ode”，这个问题不会发生。

我知道 zvode 是一种替代方案，但对于大 N 来说它似乎变得相当慢。在我遇到的真正问题中，A 是对角矩阵，但 B 是非稀疏完整矩阵。

非常感谢任何见解！我对 (i) 使用数组值参数解决这个复值 ODE 的替代方法以及 (ii) 如何让 complex_ode 运行感兴趣:)

谢谢！

Answer 1

Reti43 发布的链接似乎包含了答案，所以让我把它放在这里，以方便未来的用户：

from scipy.integrate import complex_ode
import numpy as np

N = 3
Omega = 1.0;

class myfuncs(object):
    def __init__(self, f, fargs=[]):
        self._f = f
        self.fargs=fargs

    def f(self, t, y):
        return self._f(t, y, *self.fargs)

def f(t, y, Omega,A,B):
    return -1j*(A+np.cos(Omega*t)*B).dot(y)    


A = np.random.ranf((N,N))
A = (A + A.T)/2.0

B = np.random.ranf((N,N))
B = (B + B.T)/2.0

v0 = np.zeros((N,),dtype=np.float64)
v0[0] = 1.0

t0 = 0

case = myfuncs(f, fargs=[Omega, A, B] )
solver = complex_ode(case.f)
solver.set_initial_value(v0, t0)

solver.integrate([10.0])

print solver.successful()

"""
t1 = 10
dt = 1
while solver.successful() and solver.t < t1:
    solver.integrate(solver.t+dt)
    print(solver.t, solver.y)

"""

也许有人可以评论一下为什么这个技巧能奏效吗？