使用 Scipy/Numpy-NDSolve 等效求解 Python 中的矩阵微分方程？答案

【问题标题】：Solving Matrix Differential Equation in Python using Scipy/Numpy- NDSolve equivalent?使用 Scipy/Numpy-NDSolve 等效求解 Python 中的矩阵微分方程？
【发布时间】：2020-11-19 22:05:30
【问题描述】：

我有两个 numpy 数组：9x9 和 9x1。我想在离散时间点求解微分方程，但无法让 ODEInt 工作。我不确定自己是否做对了。

使用 Mathematica，方程为：

Solution = {A[t]} /. NDSolve[{A'[t] == Ab.A[t] && A[0] == A0}, {A[t]}, {t, 0, .5}, MaxSteps -> \[Infinity]];

time = 0.25;
increment = 0.05;
MA = Table[Solution, {t, 0, time, increment}];

其中 Ab 是 9x9 矩阵，A0 是 9x1 矩阵（初始）。在这里，我解决时间和生活是美好的。

在 Python 实现中，我有以下代码给了我错误的答案：

from scipy.integrate import odeint
from numpy import array, dot, pi

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

def MatrixBM3(k12,k21,k13,k31,k23,k32,delta1,delta2,delta3,
              w1, R1, R2):

  K = array([[-k21 -k23, k12, k32, 0., 0., 0., 0., 0., 0.],
                [k21, -k12 - k13, k31, 0., 0., 0., 0., 0., 0.],
                [k23, k13, -k31 - k32, 0., 0., 0., 0., 0., 0.],
                [0., 0., 0., -k21 - k23, k12, k32, 0., 0., 0.],
                [0., 0., 0., k21, -k12 - k13, k31, 0., 0., 0.],
                [0., 0., 0., k23, k13, -k31 - k32, 0., 0., 0.],
                [0., 0., 0., 0., 0., 0., -k21 - k23, k12, k32],
                [0., 0., 0., 0., 0., 0., k21, -k12 - k13, k31],
                [0., 0., 0., 0., 0., 0., k23, k13, -k31 - k32]])

  Der = array([[0., 0., 0., -delta2, 0., 0., 0., 0., 0.],
                  [0., 0., 0., 0., -delta1, 0., 0., 0., 0.],
                  [0., 0., 0., 0., 0., -delta3, 0., 0., 0.],
                  [delta2, 0., 0., 0., 0., 0., 0., 0., 0.],
                  [0., delta1, 0., 0., 0., 0., 0., 0., 0.],
                  [0., 0., delta3, 0., 0., 0., 0., 0., 0.],
                  [0., 0., 0., 0., 0., 0., 0., 0., 0.],
                  [0., 0., 0., 0., 0., 0., 0., 0., 0.],
                  [0., 0., 0., 0., 0., 0., 0., 0., 0.]])

  W = array([[0., 0., 0., 0., 0., 0., 0., 0., 0.],
                [0., 0., 0., 0., 0., 0., 0., 0., 0.],
                [0., 0., 0., 0., 0., 0., 0., 0., 0.], 
                [0., 0., 0., 0., 0., 0., w1, 0., 0.],
                [0., 0., 0., 0., 0., 0., 0., w1, 0.],
                [0., 0., 0., 0., 0., 0., 0., 0., w1],
                [0., 0., 0., w1, 0., 0., 0., 0., 0.], 
                [0., 0., 0., 0., w1, 0., 0., 0., 0.],
                [0., 0., 0., 0., 0., w1, 0., 0., 0.]])*2*pi

  R = array([[-R2, 0., 0., 0., 0., 0., 0., 0., 0.],
                [0., -R2, 0., 0., 0., 0., 0., 0., 0.],
                [0., 0., -R2, 0., 0., 0., 0., 0., 0.],
                [0., 0., 0., -R2, 0., 0., 0., 0., 0.],
                [0., 0., 0., 0., -R2, 0., 0., 0., 0.], 
                [0., 0., 0., 0., 0., -R2, 0., 0., 0.], 
                [0., 0., 0., 0., 0., 0., -R1, 0., 0.],
                [0., 0., 0., 0., 0., 0., 0., -R1, 0.],
                [0., 0., 0., 0., 0., 0., 0., 0., -R1]]) 
  return(K + Der + W + R)

Ab = MatrixBM3(21.218791062154633, 17653.497151475527, 40.50203461096454, 93956.36617049483, 0.0, 0.0, -646.4238856161137, 6727.748368359598, 20919.132768439955, 200.0, 2.36787, 5.39681)
A0 = array([-0.001071585381162955, -0.89153191708755708, -0.00038431516707591748, 0.0, 0.0, 0.0, 0.00054009700135979673, 0.4493470361764082, 0.00019370128872934646])
time = array([0.0,0.05,0.1,0.15,0.2,0.25])
MA = odeint(deriv, A0, time, args=(Ab,), maxsteps=2000)

输出是：

[[ -1.07158538e-003  -8.91531917e-001  -3.84315167e-004   0.00000000e+000
    0.00000000e+000   0.00000000e+000   5.40097001e-004   4.49347036e-001
    1.93701289e-004]
 [  3.09311322e+019   9.45061860e+022   2.35327270e+019   2.11901406e+020
    1.63784238e+023   7.60569684e+019   2.29098804e+020   1.89766602e+023
    8.18752241e+019]
 [  9.84409730e+042   3.00774018e+046   7.48949158e+042   6.74394343e+043
    5.21257342e+046   2.42057805e+043   7.29126532e+043   6.03948436e+046
    2.60574901e+043]
 [  3.13296814e+066   9.57239028e+069   2.38359473e+066   2.14631766e+067
    1.65894606e+070   7.70369662e+066   2.32050753e+067   1.92211754e+070
    8.29301904e+066]
 [  9.97093898e+089   3.04649506e+093   7.58599405e+089   6.83083947e+090
    5.27973769e+093   2.45176732e+090   7.38521364e+090   6.11730342e+093
    2.63932422e+090]
 [  3.17333659e+113   9.69573101e+116   2.41430747e+113   2.17397307e+114
    1.68032166e+117   7.80295913e+113   2.35040739e+114   1.94688412e+117
    8.39987500e+113]]

但正确答案应该是：

{-0.0010733126291998989, -0.8929689437405254, -0.0003849346301906338, 0., 0., 0., 0.0005366563145999495, 0.4464844718702628, 0.00019246731509531696}
{-0.000591095648651598, -0.570032546156741, -0.00023381082725213798, -0.00024790706920038567, 0.00010389803046880286, -0.00005361569187144767, 0.0003273277204077012, 0.2870035216110215, 0.00012300339326137006}
{-0.0003770535829276868, -0.364106358478121, -0.0001492324135668267, -0.0001596072774600538, -0.0011479989178276948, -0.000034744485507007025, 0.00020965172928479557, 0.18378613639965447, 0.00007876820247280559}
{-0.00024100792803807562, -0.23298939195213314, -0.00009543704274825206, -0.00010271831380730501, -0.0013205519868311284, -0.000022472380871477824, 0.00013326471695185768, 0.11685506361394844, 0.00005008078740423007}
{-0.00015437993249587976, -0.1491438843823813, -0.00006111736454518403, -0.00006545797627466387, -0.0005705018939767294, -0.000014272382451480663, 0.00008455890984798549, 0.0741820536557778, 0.00003179071165818503}
{-0.00009882799610556456, -0.09529950309336405, -0.00003909275555213336, -0.00004138741286392128, 0.00006303116741431477, -8.944610716890746*^-6, 0.00005406263888971806, 0.04743157303933772, 0.00002032674776723143}

谁能指出我做错了什么？

【问题讨论】：

如果您提供一个可以运行以重现问题的最小完整示例 (stackoverflow.com/help/mcve)，将会有所帮助。例如，在你的 python 代码 sn-p 中，A0 和 Ab 没有定义。
谢谢。我只是尝试将代码编辑得尽可能少，但现在 ODEInt 函数根本不起作用。
错字：mxstep，不是maxsteps。

标签： python numpy scipy wolfram-mathematica odeint

【解决方案1】：

在调用odeint 时，尝试将tuple(array[Ab]) 更改为(array(Ab),)，甚至只是(Ab,)。也就是说，使用

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

没有看到您如何定义A0 和Ab，我不能确定这是否能解决问题，但您的代码的以下变体将起作用。我使用 3x3 数组而不是 9x9。

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,))

【讨论】：

您的示例与同一事物的 Mathematica 版本完全一致，但我没有。 " lsoda-- 在当前 t (=r1), mxstep (=i1) 在到达 tout 之前在这个调用上采取的步骤在上面的消息中，I1 = 500 在上面的消息中，R1 = 0.2155051525133E-01 在这个调用上完成了额外的工作（也许Dfun 类型错误）。”
如果你想用我们可以复制和运行的东西更新你的问题中的代码，我们可以自己重现问题，并且可能很快找到修复。同时：您的系统可能非常僵硬；尝试将参数 mxstep=2000（或更大的参数，如有必要）放入对 odeint 的调用中。
谢谢，我刚刚将它编辑为一个“工作”示例，它给出了帖子中包含的错误消息。
我现在只能建议仔细检查数组 Ab 是如何创建的。 K、Der、W 和 R 中的所有系数是否都在正确的位置？标志是否正确？在 Python 和 Mathematica 中打印 Ab，并检查它们是否相同。我怀疑 something 不同，因为我得到Ab 的特征值之一是 1082.3，这意味着系统不稳定。