且构网

分享程序员开发的那些事...
且构网 - 分享程序员编程开发的那些事

如何使用solve_ivp求解复杂的矩阵微分方程?

更新时间:1970-01-01 07:55:24

我已经更新了你的代码片段,看看下面.您应该仔细检查 doc 作为,我相信,那里的一切都很详细.

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]]