使用Scipy/Numpy解决Python中的矩阵微分方程 - NDSolve相当于？

我有两个numpy数组：9x9和9x1。我想在离散时间点求解微分方程，但是在使ODEnt工作时遇到问题。我不确定我是否做对了。使用Scipy/Numpy解决Python中的矩阵微分方程 - NDSolve相当于？

随着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.00} 
{-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}

任何人都可以指出我可能做错了什么？

来源

2014-10-11 StuckOnDiffyQ

如果您提供了可以运行以重现问题的最小完整示例（http://stackoverflow.com/help/mcve），它将有所帮助。例如，在你的python代码片段中，'A0'和'Ab'没有定义。 – 2014-10-11 23:02:49

谢谢。我只是试着编写代码尽可能小，但现在ODEInt函数根本不起作用。 – StuckOnDiffyQ 2014-10-12 02:30:29

错字：'mxstep'，而不是'maxsteps'。 – 2014-10-13 15:47:11

在致电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,))

来源

2014-10-11 22:59:20

你的例子与Mathematica版本的完全一致，但我没有。 “lsoda--在当前吨（= R1），mxstep（= I1）所采取的步骤在此呼叫到达TOUT 在上述消息之前，I1 = 500 在上述消息，R1 = 0.2155051525133E-01 过量工作完成上这个调用（也许是错误的Dfun类型）。“ – StuckOnDiffyQ 2014-10-12 17:51:18

如果您想用我们可以复制和运行的内容更新问题中的代码，我们可以自己重现问题，并且可能会很快找到修复程序。同时：你的系统可能非常僵硬;尝试在调用'odeint'时将参数'mxstep = 2000'（或者更大的参数，如果需要的话）。 – 2014-10-12 19:47:41

谢谢，我刚刚编辑它为“工作”的例子，它给出了帖子中包含的错误信息。 – StuckOnDiffyQ 2014-10-13 11:52:11

使用Scipy/Numpy解决Python中的矩阵微分方程 - NDSolve相当于？

回答

相关问题