math - 使用 sagemath 或 sympy 将 n 阶微分方程简化为一阶方程组

Question

我想将 n 阶常微分方程化简为一阶方程组。这是为数值分析做准备。我将 Sympy 和 Sagemath 用于计算机代数，但我没有在它们中找到任何函数来进行此类归约。我不确定是否还有其他人可以指出此功能是否存在于 Sympy 或 Sagemath 中。

这方面的一个例子是减少等式:

x''' - 2x'' + x' = 0 

一阶方程组:

[0  1 0 
 0  0 1 
 0 -1 2]

Answer 1

据我所知，SymPy 没有直接执行此操作的函数，但可以直接手动执行。

我假设您总是希望您的两个附加方程的形式为 y = x' 和 z = y'。

首先，让我们创建 ODE(在 SymPy 中，表达式自动假定为零，因此为简化起见，我们不要理会 with = 0 部分)。我假设您的自变量是 t。

In [4]: t = symbols('t')

In [7]: x, y, z = symbols('x y z', cls=Function)

In [6]: ode = x(t).diff(t, t) - 2*x(t).diff(t) + x(t)

In [13]: ode = x(t).diff(t, 3) - 2*x(t).diff(t, t) + x(t).diff(t)

In [14]: ode
Out[14]:
               2           3
d             d           d
──(x(t)) - 2⋅───(x(t)) + ───(x(t))
dt             2           3
             dt          dt

现在，如果我们用 y 替换 x'，

In [15]: ode.subs(x(t).diff(t), y(t))
Out[15]:
                      2
         d           d
y(t) - 2⋅──(y(t)) + ───(y(t))
         dt           2
                    dt

我们看到它还将 x'' 替换为 y'。因此，让我们用 z 替换 y':

In [16]: ode = ode.subs(x(t).diff(t), y(t)).subs(y(t).diff(t), z(t))

In [17]: ode
Out[17]:
                d
y(t) - 2⋅z(t) + ──(z(t))
                dt

现在我们的系统[x' y' z']是

In [20]: Matrix([y(t), z(t), solve(ode, z(t).diff(t))[0]])
Out[20]:
⎡     y(t)     ⎤
⎢              ⎥
⎢     z(t)     ⎥
⎢              ⎥
⎣-y(t) + 2⋅z(t)⎦

请注意，我们已经知道 x' = y 和 y' = z，所以我们直接使用它们，但我们使用 solve() 以根据 z' 获取我们的替代 ODE。

如果你想要系数，一个简单的技巧是采用雅可比行列式:

In [23]: M = Matrix([y(t), z(t), solve(ode, z(t).diff(t))[0]]).jacobian([x(t), y(t), z(t)])

In [24]: M
Out[24]:
⎡0  1   0⎤
⎢        ⎥
⎢0  0   1⎥
⎢        ⎥
⎣0  -1  2⎦

我将把它打包成一个函数 ode_to_system(ode, [x(t), y(t), z(t)]) 作为练习读者。