Sympy 和矩阵微分 ODE

Question

如何使用 Sympy 求解矩阵微分方程？

我有一个形式为 y'(t) = A*y(t) + B 的方程，其中 A 是 3x3 矩阵，y(t) 是 1x3 向量，B 是 1x3 向量。

更具体地说，我正在研究一个使用微分方程移动 3D 点的计算机图形问题。我在 3D 空间 y(t) 中有一个点、一个 3x3 旋转矩阵和一个平移向量。 y'=Ay+B 方程是我正在研究的实际问题的简化，但我无法使用 Sympy 来解决甚至 y'=Ay+B。我正在 Sympy 中寻找封闭形式的解决方案（不是数字答案）。

我知道如何解方程 y'=Ay+B，但我想使用 Sympy 找到相同的解，然后调整代码以解决我要解决的更复杂的问题。

我当前的代码是这样的：

from sympy import *

y0 = Function('y0')
y1 = Function('y1')
y2 = Function('y2')
t = symbols(('t'))
b0,b1,b2 = symbols(('b0:3'))

y = Matrix([y0(t), y1(t), y2(t)])
B = Matrix([b0,b1,b2])

ode = Eq(y.diff(t), y)

soln = dsolve(ode, y0(t),y1(t),y2(t))

但这会导致 Python 错误

TypeError: cannot add <class 'sympy.matrices.immutable.ImmutableDenseMatrix'> and <class 'sympy.core.symbol.Dummy'>

在上面的示例中，它使用了更简化的 y'=y+B 示例，但即使这样也不行。

在 Sympy 中设置此类问题的最佳方法是什么？

Answer 1

dsolve 需要平面列表或平面元组，而不是矩阵。要将列矩阵 A 转换为平面列表，可以使用咒语A.T.tolist()[0] - 即转置，转换为嵌套列表[[x, y, z]]，然后取第 0 个条目[x, y, z]。所以，你的代码应该是

ode = (y.diff(t) - y).T.tolist()[0]
soln = dsolve(ode, y.T.tolist()[0])

通常传递差异lhs-rhs而不是Eq(lhs, rhs)更方便，所以我这样做了。

遗憾的是，输出是

[Eq(y0(t), C1*exp(t)), False, False] 

（在 1.1.1 和 1.2 中），这显然是一个错误。 SymPy 中具有两个以上方程的线性 ODE 系统的解目前是基于一些大学生对他们所学的 ODE 课程的误解。 This PR 可以解决大部分问题，但已被放弃。