Hi,
I just wrote an example demonstrating the way how to write and solve matrix equations
import numpy as np
from scipy.integrate import odeint
def equations(y, t, shape=None): # in case of rectangular matrix/equations shape is (n, m), in case of rectangular shape is (n, n); we will pass shape parameter below
'''
This defines matrix equations:
dA/dt = A^2 - 4 * A + t,
where A is a matrix with n rows and n columns
(the shape of A is assumed to be n x n).
'''
A = y.reshape(*shape) #lets transform our flat array (y) to a matrix (This allows to avoid elemntwise writing of the equations)
result = A @ A - 4 * A + t # matrix equations, @ -- denotes matrix multiplication (see PEP 465)
return result.ravel() # reverse transform (to flat array)
# integration process
n = 10 # matrix size
A0 = np.random.rand(n, n) # initial conditions, i.e. A0.shape = (n, n)
y0 = A0.ravel() # initial conditions should be flat (this is needed for odeint, or ode from scipy package)
t = np.linspace(0, 10, 200) # array of points where we are planing to get the solution of the matrix equation.
solution = odeint(equations, y0, t, args=((n, n),)) # args=((n, n),) sets the shape to (n, n) in the equations above
def A(solution, j, shape):
'''Returns A(t[j])'''
return solution[j, :].reshape(shape)
print('Result at {} is {}'.format(t[10], A(solution, 10, shape=(n, n))))
However, your equation needs additional clarifications:
Suppose we have matrices S and I. The equation refers to S[s-1, i+1] elements. These elements (probably, boundary conditions) should be defined, e.g. when s=1 this requires S[0, i+1] to be defined, and when i=N (N denotes equation's size) this requires S[s-1, N + 1].
So, you need to define boundary conditions, e.g. S[0, i], S[s, N + 1], I[0, i], I[s, M+1] and shapes of the your matrices (M, N).
if S.shape is NxM and I.shape is NxM, then
1) y.shape will be 2*N*M, y = [elements of S (size=N*M), elements of I (size=N*M)]
2) transformations inside the function that defines equations should be: S = y[:N*M].reshape(N, M), I = y[N*M:].reshape(N, M)
3) you probably will need to rewrite your equations using Cython (scipy.weave, C). If loops in your equations (i.e. double sums) write in pure Python, computational process will be quite slow.