Resolviendo la ecuación diferencial de matrices en Python usando Scipy / Numpy-NDSolve equivalente?

Tengo dos matrices numpy: 9×9 y 9×1. Me gustaría resolver la ecuación diferencial en puntos de tiempo discretos, pero tengo problemas para hacer que ODEInt funcione. No estoy seguro de si estoy haciendo lo correcto.

Con Mathematica, la ecuación es:

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

Donde Ab es la matriz de 9×9, A0 es la matriz de 9×1 (inicial). Aquí, resuelvo por tiempo y la vida es buena.

En la implementación de Python tengo el siguiente código que me da la respuesta incorrecta:

 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) 

La salida es:

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

Pero la respuesta correcta debe ser:

 {-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} 

¿Puede alguien señalarme lo que puedo estar haciendo mal?

En la llamada a odeint , intente cambiar la tuple(array[Ab]) por (array(Ab),) , o incluso solo (Ab,) . Es decir, usar

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

Sin ver cómo definió A0 y Ab , no puedo estar seguro de que esto solucionará el problema, pero la siguiente variación de su código funcionará. Usé una matriz de 3×3 en lugar de 9×9.

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