How to set up a time-dependant matrix for an ODE to be solved using python?

Question

I want to solve a problem numerically in python like this:

$$ y(t)' = \mathbf{M}(t)y ,\\ y(0) = (1,0,0,0 ...) $$

where $y$ is an $n$-dimensional vector and $\mathbf{M}(t)$ is a time-dependant $n \times n$ matrix.

The matrix $\mathbf{M}(t)$ is of the shape

$$ \mathbf{M}(t)= \left[ \begin{array}{ccccc} 0 & T_1(t) & 0 & 0 & 0 & 0\\ T_1(t) & 0 & T_2(t) & 0 & 0 & 0\\ 0 & T_2(t) & 0 & T_1(t) & 0 & 0\\ 0 & 0 & T_1(t) & 0 & T_2(t) & 0\\ 0 & 0 & 0 & T_2(t) & 0 & T_1(t)\\ 0 & 0 & 0 & 0 & T_1(t) & 0\\ \end{array} \right] $$

or in numpy/scipy notation:

[[0. T1 0. 0. 0. 0.]
 [T1 0. T2 0. 0. 0.]
 [0. T2 0. T1 0. 0.]
 [0. 0. T1 0. T2 0.]
 [0. 0. 0. T2 0. T1]
 [0. 0. 0. 0. T1 0.]]

where $T_1$ and $T_2$ can be arbitrarily time dependant.

My problem is, that for the routine solve_ivp in python I need to put in a function of $y$ into the solve_ivp routine. I don't know how to create this matrix $\mathbf{M}(t)$ above so that solve_ivp can use it.

Has someone some hints or suggestions?

Answer 1

What about avoiding to construct the matrix by using its structure?

def odefunc(t,u):
    dotu = zeros_like(u)
    T1 = T1func(t)
    T2 = T2func(t)
    dotu[0::2] += T1*u[1::2]
    dotu[1::2] += T1*u[0::2]
    dotu[1:-1:2] += T2*u[2::2]
    dotu[2::2] += T2*u[1:-1:2]
    return dotu

This only works if the matrix size is even and not divisible by 4, like the current size 6.