我想将n阶常微分方程简化为一阶方程组。这是为数值分析做准备。我将Sympy和Sagemath都用于计算机代数,但是我没有在其中找到任何函数来进行这种简化。我不确定是否有人可以指出此功能是否在Sympy或Sagemath中存在。
这样的一个例子是简化等式:
x''' - 2x'' + x' = 0
到一阶方程组:
[0 1 0
0 0 1
0 -1 2]
最佳答案
据我所知,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()
来获取根据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)])
的任务留给读者练习。关于math - 使用sagemath或sympy将n阶微分方程简化为一阶方程组,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/25590481/