1

I want to solve a matrix differential equation, like this one:

import numpy as np
from scipy.integrate import odeint


def deriv(A, t, Ab):
    return np.dot(Ab, A)


Ab = np.array([[-0.25,    0,    0],
               [ 0.25, -0.2,    0],
               [    0,  0.2, -0.1]])

time = np.linspace(0, 25, 101)
A0 = np.array([10, 20, 30])

MA = odeint(deriv, A0, time, args=(Ab,))

However, this does not work in the case of having complex matrix elements. I am looking for something similar to scipy.integrate.complex_ode but for odeint. If this is not possible, what other library should I use to perform the integration? I appreciate your help!

Vladimir Vargas
  • 1,744
  • 4
  • 24
  • 48
  • Several related questions: https://stackoverflow.com/questions/26742818/solve-ode-in-python-with-complex-matrix-as-initial-value, https://stackoverflow.com/questions/26580854/how-to-plot-the-eigenvalues-when-solving-matrix-coupled-differential-equations-i, https://stackoverflow.com/questions/19910189/scipy-odeint-with-complex-initial-values (and there are more). In those questions, I suggested using the wrapper that I wrote called [`odeintw`](https://github.com/WarrenWeckesser/odeintw). – Warren Weckesser Aug 30 '17 at 20:31
  • @WarrenWeckesser This approach requires to write the matrix as a function, right? My problem is that I have 150x150 matrices and writing each them as python functions is way too time-consuming. Also computing the jacobian demands a lot of time. Is there a way of using your wrapper similar to what you proposed in that question I mention in mine? i.e. sth like `MA = odeint(deriv, A0, time, args=(Ab,))` Thanks – Vladimir Vargas Aug 30 '17 at 21:40
  • *"This approach requires to write the matrix as a function, right?"* No, do what you are already doing, but with complex coefficients in the array `Ab`. Then also make `A0` complex, and `odeintw` should Just Work. – Warren Weckesser Aug 30 '17 at 22:23
  • *"Also computing the jacobian demands a lot of time."* You don't have to provide the Jacobian function. – Warren Weckesser Aug 30 '17 at 22:28
  • @WarrenWeckesser Thanks, at first I did not make A0 complex. That solved the issue. – Vladimir Vargas Aug 30 '17 at 22:36

1 Answers1

3

odeintw wrapper for odeint must be used in the same fashion as in the question. However, the initial value A0 must be complex-valued vector.

Vladimir Vargas
  • 1,744
  • 4
  • 24
  • 48
  • 1
    Actually, `Ab` does not have to be complex. I only suggested that in my comment to the question because you said your actual array would be complex. It is the initial condition `A0` that *must* be complex. That's how `odeintw` knows you are using complex variables. – Warren Weckesser Aug 30 '17 at 22:58