使用 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'。所以让我们用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() 来获取我们的替换 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)]) 的任务留给读者作为练习。

Answer 2

我写了一个实验库来处理常微分方程组：

https://github.com/bjodah/symodesys

它基于 sympy，不幸的是我没有写太多文档，但我提供了很多示例。我将按如下方式处理您的方程式：

from sympy import *
from symodesys.odesys import AnyOrderODESystem

t = Symbol('t')
x = Function('x')(t)

D1x = x.diff(t)
D2x = x.diff(t, 2)
D3x = x.diff(t, 3)

expr = Eq(D3x, 2*D2x - D1x)

odesys = AnyOrderODESystem.from_list_of_eqs([expr])
print(odesys.all_depv)
redsys = odesys.reduce_to_sys_of_first_order()
print(redsys.all_depv)
print(redsys.f)

哪些输出：

[x(t)]
[x(t), x_h1(t), x_h2(t)]
OrderedDict([(x(t), x_h1(t)), (x_h1(t), x_h2(t)), (x_h2(t), -x_h1(t) + 2*x_h2(t))])

添加一些额外的行给你一个 gui 来试验一个初始值问题 （见解曲线作为初始值的函数）

from symodesys.gui import get_chaco_viewer
from collections import defaultdict
y0 = defaultdict(int)
y0[redsys['x']] = 3.14
params = {}
t0, tend, N = 0, 10, 100
viewer = get_chaco_viewer(redsys, y0, params, t0, tend, N)
viewer.configure_traits()
viewer.clean_up()

给你：

安装一些依赖有点棘手，需要帮助的可以添加评论！

Answer 3

问题是，首先如何将微分方程编码为字符串。因为编码已经比写下一阶系统复杂多了。

一般的手动程序是设置 x1=x, x2=x', x3=x'' 然后注意

x1'=x'=x2
x2'=x''=x3
x3'=x'''= 2*x'' - x' = 2*x3 - x2

然后将得到的系统转换为矩阵形式。

另请参阅多项式的伴随矩阵，这也是（直到转置）您为高阶线性微分方程的系统矩阵获得的一般形式。

Answer 4

我不使用 Sympy 或 Sagemath。但在他们的 API 文档中查找 Laplace 或 Z 变换。

• 如果找到了，那么你要做的工作就更少了
• 如果没有，您将不得不自己查找 lib 或编写代码

利用拉普拉斯求解微分系统

• 我有一段时间没用这个了，所以用一些数学书检查一下！！！
• 无论如何，如果我没记错的话
• 拉普拉斯变换将积分函数转换为线性函数（时域到s域）
• 您的微分函数只需要具有连续的导数/积分即可完成这项工作

• 要解决您的问题，请执行以下操作：

  1. 通过积分整个事物将所有微分转换为积分
  2. 应用拉普拉斯变换
     • 这会将微分系统转换为多项式系统
     • wiki example
  3. 求解多项式方程组
  4. 应用拉普拉斯逆变换
     • 这会将部分结果转换为您的解决方案的结果
  5. 通过探针定义的边缘情况求解积分常数

• 另请参阅here 以获取其他示例

• 关于这个主题的内容很多，只需 google 一下

用Z求解微分系统

• 从未这样做过，但如果与拉普拉斯变换的解决方案不同，则应该是相似的