python complex_ode pass matrix-valued parameters

Question

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!

Can you show us the full traceback please? The `TypeError` seems to be the same as [here](https://stackoverflow.com/questions/34577870/using-scipy-integrate-complex-ode-instead-of-scipy-integrate-ode). In short, do not use `set_f_params()` with `complex_ode()`. For solving complex ODEs, you can also use `zvode()`. That takes care of problem number 2. — Reti43, Jan 21 '16 at 10:32
I updated the post showing the full error message in both cases. I know of zvode but it is sort of slow when the matrix size gets big. Or do u mean that there is a way to combine complex_ode() and zvode just to pass the parameters? — python_freak, Jan 21 '16 at 11:23
That's the bug, alright. If you print the values that are passed in `_wrap()`, you'll see that they are all mangled up if you have used `set_f_params()`, or modified `solver.params` yourself (which is what `set_f_params()` does internally). So where you're supposed to be slicing arrays in `_wrap()`, you end up having floats, hence the `TypeError`. No, there is no combining the two integrators, you simply use `zvode` from the get go. Why shouldn't the integrator become slower when the matrix gets bigger? If the matrices are supposed to hold unknown parameters, their number increases by N**2. — Reti43, Jan 21 '16 at 11:39
I seem to have misunderstood part of your problem (through no fault of your own, no worries). Are the parameters `Omega`, `A` and `B` constant? If yes, your 2 questions boil down to *"how to pass parameters in complex_ode"*, for which the answer is in the link. To clarify, you can't use `set_f_params()` without breaking `complex_ode()`. The only solution is to create the function as `def f(t, y)` and let `Omega`, `A` and `B` be fetched from the globals. The solution in the link achieves pretty much that but by creating a wrapper function. Also, `zvode` is faster than `complex_ode()`. — Reti43, Jan 21 '16 at 14:36

score 1 · Accepted Answer · answered Jan 21 '16 at 12:35

1

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?

answered Jan 21 '16 at 12:35

python_freak

259
3
12

See my comment above about passing a function of only `t` and `y` to `complex_ode()`. If you answer my question whether `Omega`, `A` and `B` are constant parameters, then we can resolve this as a duplicate. By linking to the other question, future readers will be able to see the solution without having to post a duplicate one here. – Reti43 Jan 21 '16 at 14:42
yes, Omega, A and B a t-independent. I think this is a duplicate, so it can be removed. – python_freak Jan 21 '16 at 17:24
Don't remove the question. Duplicates are useful because they effectively introduce more keywords that people can use in their search to stumble upon the original question/answer. – Reti43 Jan 21 '16 at 17:30

python complex_ode pass matrix-valued parameters

1 Answers1