MATLAB 的 ode15 如何处理包含值矩阵的方程的积分？

Question

我正在将 MATLAB 代码翻译成 Python，但在担心翻译之前，我想了解 MATLAB（特别是其 ODE15s 求解器）如何解释方程。

我有一个函数脚本，在主脚本中调用它，该函数脚本包含方程：

function testFun=testFunction(t,f,dmat,releasevec)

    testFun=(dmat*f)+(releasevec.');

在 testFunction 中，t 指时间，f 指我正在求解的值，dmat 指我感兴趣的常量矩阵，releasevec 指附加常量的向量。

主脚本中的 ODE15s 求解器通过以下几行发挥其魔力：

for i=1:1461

    [~,f]=ode15s(@(t, f) testFunction(t, f, ...
    [dAremoval(i), dFWtoA(i), dSWtoA(i), dStoA(i), dFSedtoA(i), dSSedtoA(i); ...
    dAtoFW(i), dFWremoval(i), dSWtoFW(i), dStoFW(i), dFSedtoFW(i), dSSedtoFW(i); ...
    dAtoSW(i), dFWtoSW(i), dSWremoval(i), dStoSW(i), dFSedtoSW(i), dSSedtoSW(i); ...
    dAtoS(i), dFWtoS(i), dSWtoS(i), dSremoval(i), dFSedtoS(i), dSSedtoS(i); ...
    dAtoFSed(i), dFWtoFSed(i), dSWtoFSed(i), dStoFSed(i), dFSedremoval(i), dSSedtoFSed(i); ...
    dAtoSSed(i), dFWtoSSed(i), dSWtoSSed(i), dStoSSed(i), dFSedtoSSed(i), dSSedremoval(i)], ...
    [Arelease(i), FWrelease(i), SWrelease(i), Srelease(i), FSedrelease(i), SSedrelease(i)]), [i, i+1], fresults(:, i),options);
    fresults(:, i + 1) = f(end, :).';

fresults

是一个最初由零组成的表，其中包含 f 结果。这些选项调用 odeset 来获取“非负”值。上面的 d 值矩阵是一个 6x6 矩阵。我已经计算了所有 d 值和释放值。我的问题是：ode15s 如何与 testfunction 方程中给出的 6x6 矩阵进行积分？我试图手动解决这个问题，但没有成功。任何帮助将不胜感激！！

#

def func(y, t, params):
    f = 5.75e-16
    f = y
    dmat, rvec = params
    derivs = [(dmat*f)+rvec]
    return derivs

# Parameters
dmat = np.array([[-1964977.10876756, 58831.976165, 39221.31744333,  1866923.81515922, 0.0, 0.0],
             [58831.976165, -1.89800738e+09, 0.0, 1234.12447489, 21088.06180415, 14058.70786944],
             [39221.31744333, 0.84352331, -7.59182852e+09, 0.0, 0.0, 0.0],
             [1866923.81515922, 0.0, 0.0, -9.30598884e+08, 0.0, 0.0],
             [0.0, 21088.10183616, 0.0, 0.0, -1.15427010e+09, 0.0],
             [0.0, 0.0, 14058.73455744, 0.0, 0.0, -5.98519566e+09]], np.float)
new_d = np.ndarray.flatten(dmat)

rvec = np.array([[0.0], [0.0], [0.0], [0.0], [0.0], [0.0]])

f = 5.75e-16

# Initial conditions for ODE
y0 = f

# Parameters for ODE 
params = [dmat, rvec]

# Times
tStop = 2.0
tStart = 0.0
tStep = 1.0

t = np.arange(tStart, tStop, tStep)

# Call the ODE Solver
soln = odeint(func, y0, t, args=(params,))
#y = odeint(lambda y, t: func(y,t,params), y0, t)

Answer 1

它说here ode15s 使用后向差分公式进行微分。 你的微分方程是（据我所知）

f' = testFunc(t,f)

，它在函数内部有一些向量矩阵计算。 然后你可以用后向差分公式来代替微分，即：

f_next = f_prev + h*testFunc(t,f_next);

其中 f_prev 是向量的初始值。这里的计算没有重要的区别，只是因为 testFunc(t,f) 函数包含一个 6x6 矩阵。每次它解决逆问题时，都会通过数字创建雅可比矩阵来找到 f_next。

然而，尝试像 matlab 那样编写算法可能比我们想象的更难，因为 matlab 对问题有一些（与优化相关或无关的）特殊处理。您应该小心您获得的每个值。

Answer 2

本质上，你需要改变的东西很少。使用

numpy.ndarray

表示向量和矩阵。时间步进可以使用

scipy.integrate.ode

来完成。对于 ODE 函数的每次更改，您都需要重新初始化积分器，或者通过

set_f_parameter

提供矩阵和参数作为附加函数参数。

更接近 matlab 界面但仅限于 lsoda 的是

scipy.integrate.odeint

。但是，由于您使用求解器来解决刚性问题，因此这可能正是您所需要的。