CODEWITHSUNDEEP
python
python_freak
python_freak

Reputation: 259

python complex_ode pass matrix-valued parameters

I'm having some trouble with python's complex_ode solver.

I'm trying to solve the following equation:

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

where A and B are NxN arrays and the unknown y is an Nx1 array, i is the imaginary unit and Omega is a parameter.

Here's my code:

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()

Running this script produces the error

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)

If instead I use solver = ode(f), ie. the real-valued solver, it runs fine. Except that it doesn't solve the ODE I want which is complex-valued :(

I then tried to reduce the number of parameters by making the matrices A and B global variables. This way the only parameter the function f accepts is Omega. The error changes to

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

where I figured out that float refers to the parameter Omega [by trying an integer]. Again, "ode" alone works in this case as well.

Last, I tried the same complex valued equation, but now A and B are just numbers. I tried to pass them both as parameters, i.e. params = [Omega,A,B], as well as making them global variables in which case params=[Omega]. The error is the 

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

error - the full error is the same as above. And once again this problem does not occur for the real-valued "ode".

I know zvode is an alternative, but it appears to become quite slow for large N. In the real problem I have, A is a diagonal matrix but B is a non-sparse full matrix.

Any insights are much appreciated! I'm interested both in (i) alternative ways to solve this complex-valued ODE with array-valued parameters, and (ii) how to get complex_ode to run :)

Thanks!

Upvotes: 3

Views: 515

Answers (1)

python_freak
python_freak

Reputation: 259

It seems like the link that Reti43 posted contains the answer, so let me put it here for the benefit of future users:

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)

"""

Could maybe someone comment on why this trick does the job?

Upvotes: 1

Related Questions