时间戳

首页 > 解决方案 > Sympy 和矩阵微分 ODE

python - Sympy 和矩阵微分 ODE

问题描述

如何使用 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 中设置此类问题的最佳方法是什么？

标签： pythonsympysymbolic-math

解决方案

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 课程的误解。这个 PR可以解决大部分问题，但它被放弃了。

推荐阅读