1

I want to solve a complex matrix differential equation y' = Ay.

import numpy as np
from scipy.integrate import solve_ivp


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


A = np.array([[-0.25 + 0.14j,    0,    0.33 + 0.44j],
              [ 0.25 + 0.58j, -0.2 + 0.14j,    0],
              [    0,  0.2 + 0.4j, -0.1 + 0.97j]])

time = np.linspace(0, 25, 101)
y0 = np.array([[2, 3, 4], [5, 6 , 7], [9, 34, 78]])

result = solve_ivp(deriv, y0, time, args=(A,))

There already seems to be an answer in case of 'odeint'. https://stackoverflow.com/a/45970853/7952027

https://stackoverflow.com/a/26320130/7952027

https://stackoverflow.com/a/26747232/7952027

https://stackoverflow.com/a/26582411/7952027

I am curious as to whether it can be done with any of the new API of Scipy?

Tejas Shetty
  • 685
  • 6
  • 30

1 Answers1

1

I have updated your snippet, have a look below. You should carefully check the doc as, I believe, everything is well detailed there.

import numpy as np
from scipy.integrate import solve_ivp


def deriv_vec(t, y):
    return A @ y


def deriv_mat(t, y):
    return (A @ y.reshape(3, 3)).flatten()


A = np.array([[-0.25 + 0.14j, 0, 0.33 + 0.44j],
              [0.25 + 0.58j, -0.2 + 0.14j, 0],
              [0, 0.2 + 0.4j, -0.1 + 0.97j]])

result = solve_ivp(deriv_vec, [0, 25], np.array([10 + 0j, 20 + 0j, 30 + 0j]),
                   t_eval=np.linspace(0, 25, 101))
print(result.y[:, 0])
# [10.+0.j 20.+0.j 30.+0.j]
print(result.y[:, -1])
# [18.46+45.25j 10.01+36.23j -4.98+80.07j]

y0 = np.array([[2 + 0j, 3 + 0j, 4 + 0j],
               [5 + 0j, 6 + 0j, 7 + 0j],
               [9 + 0j, 34 + 0j, 78 + 0j]])
result = solve_ivp(deriv_mat, [0, 25], y0.flatten(),
                   t_eval=np.linspace(0, 25, 101))
print(result.y[:, 0].reshape(3, 3))
# [[ 2.+0.j  3.+0.j  4.+0.j]
#  [ 5.+0.j  6.+0.j  7.+0.j]
#  [ 9.+0.j 34.+0.j 78.+0.j]]
print(result.y[:, -1].reshape(3, 3))
# [[ 5.67+12.07j 17.28+31.03j 37.83+63.25j]
#  [ 3.39+11.82j 21.32+44.88j 53.17+103.80j]
#  [ -2.26+22.19j -15.12+70.191j -38.34+153.29j]]
Patol75
  • 4,342
  • 1
  • 17
  • 28
  • Actually, I want to know if the new API can handle a complex matrix differential equation, here you have done it for a complex vector y. – Tejas Shetty Jan 22 '21 at 09:52
  • 1
    Well, `y` has to be a vector, but you can always work with matrices within the function, then flatten the result in the return call, and finally reshape the solution after `solve_ivp` finishes. Or am I missing something? – Patol75 Jan 23 '21 at 14:16
  • You mean first, flatten the matrix, use solve_ivp and then reshape the result? – Tejas Shetty Jan 24 '21 at 05:55
  • 1
    As long as you provide `solve_ivp` with a vector both in the original call and in the function's return call, you will be good. You can do any reshape within the function to perform matrix operations. – Patol75 Jan 24 '21 at 22:19
  • What can I do if I have to use a matrix itself? The differential equation is not amenable to this flattening of the matrix to the vector. – Tejas Shetty Jan 25 '21 at 06:40
  • 1
    I am not sure you understood my point. I have updated the snippet. Please let me know if it is clearer now. – Patol75 Jan 25 '21 at 08:06
  • much clearer @Patol75 But can you read the comments of https://stackoverflow.com/q/65880726/7952027 and let me know what you think. – Tejas Shetty Jan 25 '21 at 09:51
  • 1
    From my experience, I have always benefited from comments and answers by both @Lutz Lehmann and @Warren Weckesser, so I think you should listen to them. My additional advice is for you to read the source code of the `solve_ivp` framework, I found it quite intuitive, and it contains references to published works. – Patol75 Jan 26 '21 at 04:27
  • So the odeintw package achieves the same objective as your @patol75 answer without using solve_ivp. Am I right? – Tejas Shetty Jan 26 '21 at 04:39