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)

如果我改为使用 solver = 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)

"""

也许有人可以评论为什么这个技巧能起作用？