使用 Scipy/Numpy-NDSolve 等效求解 Python 中的矩阵微分方程? [英] Solving Matrix Differential Equation in Python using Scipy/Numpy- NDSolve equivalent?
问题描述
我有两个 numpy 数组:9x9 和 9x1.我想在离散时间点求解微分方程,但无法让 ODEInt 工作.我不确定我是否在做正确的事情.
使用 Mathematica,等式是:
Solution = {A[t]}/.NDSolve[{A'[t] == Ab.A[t] &&A[0] == A0}, {A[t]}, {t, 0, .5}, MaxSteps ->\[无限]];时间 = 0.25;增量 = 0.05;MA = Table[Solution, {t, 0, time, increment}];
其中 Ab 是 9x9 矩阵,A0 是 9x1 矩阵(初始).在这里,我解决了时间和生活的美好.
在 Python 实现中,我有以下代码给出了错误的答案:
from scipy.integrate import odeint从 numpy 导入数组,点,pidef deriv(A, t, Ab):返回点(Ab,A)def MatrixBM3(k12,k21,k13,k31,k23,k32,delta1,delta2,delta3,w1、R1、R2):K = 数组([[-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 = 数组([[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*piR = 数组([[-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]])回报(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 =阵列([ - 0.001071585381162955,-0.89153191708755708,-0.00038431516707591748,0.0,0.0,0.0,0.00054009700135979673,0.4493470361764082,0.00019370128872934646])时间 = 数组([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+0000.00000000e+000 0.00000000e+000 5.40097001e-004 4.49347036e-0011.93701289e-004][ 3.09311322e+019 9.45061860e+022 2.35327270e+019 2.11901406e+0201.63784238e+023 7.60569684e+019 2.29098804e+020 1.89766602e+0238.18752241e+019][ 9.84409730e+042 3.00774018e+046 7.48949158e+042 6.74394343e+0435.21257342e+046 2.42057805e+043 7.29126532e+043 6.03948436e+0462.60574901e+043][ 3.13296814e+066 9.57239028e+069 2.38359473e+066 2.14631766e+0671.65894606e+070 7.70369662e+066 2.32050753e+067 1.92211754e+0708.29301904e+066][ 9.97093898e+089 3.04649506e+093 7.58599405e+089 6.83083947e+0905.27973769e+093 2.45176732e+090 7.38521364e+090 6.11730342e+0932.63932422e+090][ 3.17333659e+113 9.69573101e+116 2.41430747e+113 2.17397307e+1141.68032166e+117 7.80295913e+113 2.35040739e+114 1.94688412e+1178.39987500e+113]]
但正确答案应该是:
<预> <代码> { - 0.0010733126291998989,-0.8929689437405254,-0.0003849346301906338,0.1,0.1,0.1,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}谁能指出我可能做错了什么?
在对 odeint
的调用中,尝试将 tuple(array[Ab])
更改为 (array(Ab),)
,甚至只是 (Ab,)
.即使用
MA = odeint(deriv, A0, time, (Ab,))
在没有看到您如何定义 A0
和 Ab
的情况下,我不能确定这会解决问题,但是您的代码的以下变体将起作用.我使用了 3x3 数组而不是 9x9.
将 numpy 导入为 np从 scipy.integrate 导入 odeintdef deriv(A, t, Ab):返回 np.dot(Ab, A)Ab = np.array([[-0.25, 0, 0],[ 0.25, -0.2, 0],[ 0, 0.2, -0.1]])时间 = np.linspace(0, 25, 101)A0 = np.array([10, 20, 30])MA = odeint(deriv, A0, 时间, args=(Ab,))
I have two numpy arrays: 9x9 and 9x1. I'd like to solve the differential equation at discrete time points, but am having trouble getting ODEInt to work. I do am unsure if I'm even doing the right thing.
With Mathematica, the equation is:
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}];
Where Ab is the 9x9 matrix, A0 is the 9x1 matrix (initial). Here, I solve for time and life is good.
In Python implementation I have the following code which gives me the wrong answer:
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)
Output is:
[[ -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]]
But the correct answer should be:
{-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}
Can anyone point me to what I may be doing wrong?
In the call to odeint
, try changing tuple(array[Ab])
to (array(Ab),)
, or even just (Ab,)
. That is, use
MA = odeint(deriv, A0, time, (Ab,))
Without seeing how you defined A0
and Ab
, I can't be sure that this will fix the problem, but the following variation of your code will work. I used a 3x3 array instead of 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,))
这篇关于使用 Scipy/Numpy-NDSolve 等效求解 Python 中的矩阵微分方程?的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!