python - 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 类(class)的误解。 This PR会解决大部分问题，但它被放弃了。