蟒蛇complex_ode通矩阵值参数

Question

我遇到一些麻烦python的complex_ode求解。

我试图解决以下等式：

DY/DT = -iAy - ICOS（Omegat）乙ý

其中A和B是N×N个阵列和未知y是Nx1阵列，i是虚数单位，Ω是参数。

这里是我的代码：

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接受的唯一参数是欧米茄，错误修改

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__' 

，我想通了，浮指参数欧米茄[试图通过一个整数。同样，在这种情况下，“颂”单独的作品也是如此。

末，我试图值方程相同的配合物，但现在A和B仅仅是数字。我试图他们两个作为参数传递，即PARAMS =ω，A，B]为WEL l使它们成为全局变量，在这种情况下，params =Ω。错误是

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

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

我知道zvode是一种选择，但它似乎成为大N.在真正的问题我已经很慢，A是对角矩阵，但B是一个非稀疏全矩阵。

任何见解都非常感谢！我在（我）的替代方法来解决与数组值参数此复值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) 

""" 

莫非为什么这招没有工作也许有人对此有何评论？