使用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中,表达式自动被假定为等于零,所以为了简化事情,让我们不要打扰= 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

现在,如果我们更换x'用y,

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'.所以我们y'用z以下代替:

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()得到我们的取代ODE来讲z'.

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

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)])作为练习给读者.