SymPy 1.13 中文文档（二）

原文：docs.sympy.org/latest/index.html

解算器

原文：docs.sympy.org/latest/tutorials/intro-tutorial/solvers.html

注意

对于解决常见类型方程的初学者友好指南，请参阅解方程。

>>> from sympy import *
>>> x, y, z = symbols('x y z')
>>> init_printing(use_unicode=True) 

方程的注意事项

从本教程的陷阱部分回想起，SymPy 中的符号方程不是用=或==表示，而是用Eq表示。

>>> Eq(x, y)
x = y 

不过，还有一种更简单的方法。在 SymPy 中，解函数会自动假设任何不在Eq中的表达式等于零。因此，由于(a = b)当且仅当(a - b = 0)，这意味着，不需要使用x == y，只需使用x - y。例如

>>> solveset(Eq(x**2, 1), x)
{-1, 1}
>>> solveset(Eq(x**2 - 1, 0), x)
{-1, 1}
>>> solveset(x**2 - 1, x)
{-1, 1} 

如果要解的方程已经等于 0，则无需输入solveset(Eq(expr, 0), x)，可以直接使用solveset(expr, x)。

代数解方程

解代数方程的主要函数是solveset。solveset的语法是solveset(equation, variable=None, domain=S.Complexes)，其中equations可以是Eq实例或被假定为等于零的表达式。

请注意还有另一个名为solve的函数，也可用于解方程。其语法是solve(equations, variables)，但推荐使用solveset。

当解单个方程时，solveset的输出是解的FiniteSet或Interval或ImageSet。

>>> solveset(x**2 - x, x)
{0, 1}
>>> solveset(x - x, x, domain=S.Reals)
ℝ
>>> solveset(sin(x) - 1, x, domain=S.Reals)
⎧        π │      ⎫
⎨2⋅n⋅π + ─ │ n ∊ ℤ⎬
⎩        2 │      ⎭ 

如果没有解，则返回EmptySet，如果无法找到解，则返回ConditionSet。

>>> solveset(exp(x), x)     # No solution exists
∅
>>> solveset(cos(x) - x, x)  # Not able to find solution
{x │ x ∊ ℂ ∧ (-x + cos(x) = 0)} 

在solveset模块中，使用linsolve来解线性方程组。将来我们将能够直接从solveset使用linsolve。以下是linsolve语法的示例。

• 方程的列表形式：

  >>> linsolve([x + y + z - 1, x + y + 2*z - 3 ], (x, y, z))
  {(-y - 1, y, 2)} 
  

• 增广矩阵形式：

  >>> linsolve(Matrix(([1, 1, 1, 1], [1, 1, 2, 3])), (x, y, z))
  {(-y - 1, y, 2)} 
  

• A*x = b 形式

  >>> M = Matrix(((1, 1, 1, 1), (1, 1, 2, 3)))
  >>> system = A, b = M[:, :-1], M[:, -1]
  >>> linsolve(system, x, y, z)
  {(-y - 1, y, 2)} 
  

注意

解的顺序对应于给定符号的顺序。

在solveset模块中，使用nonlinsolve来解非线性方程组。以下是nonlinsolve的示例。

1. 当只有实数解时：

   >>> a, b, c, d = symbols('a, b, c, d', real=True)
   >>> nonlinsolve([a**2 + a, a - b], [a, b])
   {(-1, -1), (0, 0)}
   >>> nonlinsolve([x*y - 1, x - 2], x, y)
   {(2, 1/2)} 
   

2. 当只有复数解时：

   >>> nonlinsolve([x**2 + 1, y**2 + 1], [x, y])
   {(-ⅈ, -ⅈ), (-ⅈ, ⅈ), (ⅈ, -ⅈ), (ⅈ, ⅈ)} 
   

3. 当既有实数解又有复数解时：

   >>> from sympy import sqrt
   >>> system = [x**2 - 2*y**2 -2, x*y - 2]
   >>> vars = [x, y]
   >>> nonlinsolve(system, vars)
   {(-2, -1), (2, 1), (-√2⋅ⅈ, √2⋅ⅈ), (√2⋅ⅈ, -√2⋅ⅈ)} 
   

   >>> system = [exp(x) - sin(y), 1/y - 3]
   >>> nonlinsolve(system, vars)
   {({2⋅n⋅ⅈ⋅π + log(sin(1/3)) │ n ∊ ℤ}, 1/3)} 
   

4. 当系统是正维度系统（有无限多个解）时：

   >>> nonlinsolve([x*y, x*y - x], [x, y])
   {(0, y)} 
   

   >>> system = [a**2 + a*c, a - b]
   >>> nonlinsolve(system, [a, b])
   {(0, 0), (-c, -c)} 
   

注意

1. 解的顺序对应于给定符号的顺序。

2. 目前nonlinsolve不会以LambertW形式返回解（如果解以LambertW形式存在）。

solve可以用于这些情况：

>>> solve([x**2 - y**2/exp(x)], [x, y], dict=True)
⎡⎧         ____⎫  ⎧        ____⎫⎤
⎢⎨        ╱  x ⎬  ⎨       ╱  x ⎬⎥
⎣⎩y: -x⋅╲╱  ℯ  ⎭, ⎩y: x⋅╲╱  ℯ  ⎭⎦
>>> solve(x**2 - y**2/exp(x), x, dict=True)
⎡⎧      ⎛-y ⎞⎫  ⎧      ⎛y⎞⎫⎤
⎢⎨x: 2⋅W⎜───⎟⎬, ⎨x: 2⋅W⎜─⎟⎬⎥
⎣⎩      ⎝ 2 ⎠⎭  ⎩      ⎝2⎠⎭⎦ 

3. 目前nonlinsolve无法正确解决具有三角函数的方程组。

solve可以用于这些情况（但不提供所有解）：

>>> solve([sin(x + y), cos(x - y)], [x, y])
⎡⎛-3⋅π   3⋅π⎞  ⎛-π   π⎞  ⎛π  3⋅π⎞  ⎛3⋅π  π⎞⎤
⎢⎜─────, ───⎟, ⎜───, ─⎟, ⎜─, ───⎟, ⎜───, ─⎟⎥
⎣⎝  4     4 ⎠  ⎝ 4   4⎠  ⎝4   4 ⎠  ⎝ 4   4⎠⎦ 

solveset仅报告每个解一次。要获取多重性的多项式解，请使用roots。

>>> solveset(x**3 - 6*x**2 + 9*x, x)
{0, 3}
>>> roots(x**3 - 6*x**2 + 9*x, x)
{0: 1, 3: 2} 

输出 {0: 1, 3: 2} 的 roots 意味着 0 是多重根为 1，而 3 是多重根为 2。

注意

目前 solveset 无法解决以下类型的方程：

• 可以用 LambertW（超越方程求解器）解的方程。

对于这类情况，可以使用 solve：

>>> solve(x*exp(x) - 1, x )
[W(1)] 

解微分方程

要解微分方程，请使用 dsolve。首先通过将 cls=Function 传递给 symbols 函数创建一个未定义函数。

>>> f, g = symbols('f g', cls=Function) 

f 和 g 现在是未定义的函数。我们可以调用 f(x)，它将代表一个未知函数。

>>> f(x)
f(x) 

f(x) 的导数是未计算的。

>>> f(x).diff(x)
d
──(f(x))
dx 

（详见导数部分，了解更多关于导数的内容）。

要表示微分方程 (f''(x) - 2f'(x) + f(x) = \sin(x))，我们可以使用以下方式：

>>> diffeq = Eq(f(x).diff(x, x) - 2*f(x).diff(x) + f(x), sin(x))
>>> diffeq
 2
 d           d
f(x) - 2⋅──(f(x)) + ───(f(x)) = sin(x)
 dx           2
 dx 

要解 ODE，请将其和要解的函数传递给 dsolve。

>>> dsolve(diffeq, f(x))
 x   cos(x)
f(x) = (C₁ + C₂⋅x)⋅ℯ  + ──────
 2 

dsolve 返回 Eq 的一个实例。这是因为一般来说，微分方程的解不能显式地解出函数。

>>> dsolve(f(x).diff(x)*(1 - sin(f(x))) - 1, f(x))
x - f(x) - cos(f(x)) = C₁ 

解中的任意常数来自 dsolve 的解，是形如 C1、C2、C3 等符号。

矩阵

原文：docs.sympy.org/latest/tutorials/intro-tutorial/matrices.html

>>> from sympy import *
>>> init_printing(use_unicode=True) 

要在 SymPy 中创建矩阵，请使用Matrix对象。通过提供构成矩阵的行向量列表来构造矩阵。例如，要构造矩阵

[\begin{split}\left[\begin{array}{cc}1 & -1\3 & 4\0 & 2\end{array}\right]\end{split}]

使用

>>> Matrix([[1, -1], [3, 4], [0, 2]])
⎡1  -1⎤
⎢     ⎥
⎢3  4 ⎥
⎢     ⎥
⎣0  2 ⎦ 

要轻松创建列向量，列表中的元素被视为列向量。

>>> Matrix([1, 2, 3])
⎡1⎤
⎢ ⎥
⎢2⎥
⎢ ⎥
⎣3⎦ 

矩阵可以像 SymPy 或 Python 中的任何其他对象一样进行操作。

>>> M = Matrix([[1, 2, 3], [3, 2, 1]])
>>> N = Matrix([0, 1, 1])
>>> M*N
⎡5⎤
⎢ ⎥
⎣3⎦ 

SymPy 矩阵的一个重要特点是，与 SymPy 中的其他对象不同，它们是可变的。这意味着它们可以就地修改，如下面将看到的。这样做的缺点是，Matrix不能用于需要不可变性的地方，例如 SymPy 表达式内部或作为字典的键。如果需要一个不可变版本的Matrix，请使用ImmutableMatrix。

基本操作

这里是对Matrix的一些基本操作。

形状

要获取矩阵的形状，请使用shape()函数。

>>> from sympy import shape
>>> M = Matrix([[1, 2, 3], [-2, 0, 4]])
>>> M
⎡1   2  3⎤
⎢        ⎥
⎣-2  0  4⎦
>>> shape(M)
(2, 3) 

访问行和列

要获取矩阵的单独行或列，请使用row或col。例如，M.row(0)将获取第一行。M.col(-1)将获取最后一列。

>>> M.row(0)
[1  2  3]
>>> M.col(-1)
⎡3⎤
⎢ ⎥
⎣4⎦ 

删除和插入行和列

要删除行或列，请使用row_del或col_del。这些操作会就地修改矩阵。

>>> M.col_del(0)
>>> M
⎡2  3⎤
⎢    ⎥
⎣0  4⎦
>>> M.row_del(1)
>>> M
[2  3] 

要插入行或列，请使用row_insert或col_insert。这些操作不会在原地执行。

>>> M
[2  3]
>>> M = M.row_insert(1, Matrix([[0, 4]]))
>>> M
⎡2  3⎤
⎢    ⎥
⎣0  4⎦
>>> M = M.col_insert(0, Matrix([1, -2]))
>>> M
⎡1   2  3⎤
⎢        ⎥
⎣-2  0  4⎦ 

除非明确说明，下文提到的方法均不在原地操作。通常情况下，不在原地操作的方法将返回一个新的Matrix，而在原地操作的方法将返回None。

基本方法

如上所述，简单的操作如加法、乘法和乘幂只需使用+、*和**。要找到矩阵的逆，只需将其提升到-1次幂。

>>> M = Matrix([[1, 3], [-2, 3]])
>>> N = Matrix([[0, 3], [0, 7]])
>>> M + N
⎡1   6 ⎤
⎢      ⎥
⎣-2  10⎦
>>> M*N
⎡0  24⎤
⎢     ⎥
⎣0  15⎦
>>> 3*M
⎡3   9⎤
⎢     ⎥
⎣-6  9⎦
>>> M**2
⎡-5  12⎤
⎢      ⎥
⎣-8  3 ⎦
>>> M**-1
⎡1/3  -1/3⎤
⎢         ⎥
⎣2/9  1/9 ⎦
>>> N**-1
Traceback (most recent call last):
...
NonInvertibleMatrixError: Matrix det == 0; not invertible. 

要对矩阵进行转置，请使用T。

>>> M = Matrix([[1, 2, 3], [4, 5, 6]])
>>> M
⎡1  2  3⎤
⎢       ⎥
⎣4  5  6⎦
>>> M.T
⎡1  4⎤
⎢    ⎥
⎢2  5⎥
⎢    ⎥
⎣3  6⎦ 

矩阵构造函数

存在多个构造函数用于创建常见矩阵。要创建单位矩阵，请使用eye。eye(n)将创建一个大小为(n\times n)的单位矩阵。

>>> eye(3)
⎡1  0  0⎤
⎢       ⎥
⎢0  1  0⎥
⎢       ⎥
⎣0  0  1⎦
>>> eye(4)
⎡1  0  0  0⎤
⎢          ⎥
⎢0  1  0  0⎥
⎢          ⎥
⎢0  0  1  0⎥
⎢          ⎥
⎣0  0  0  1⎦ 

要创建全零矩阵，请使用zeros。zeros(n, m)创建一个大小为(n\times m)的全为(0)的矩阵。

>>> zeros(2, 3)
⎡0  0  0⎤
⎢       ⎥
⎣0  0  0⎦ 

类似地，ones创建一个全为 1 的矩阵。

>>> ones(3, 2)
⎡1  1⎤
⎢    ⎥
⎢1  1⎥
⎢    ⎥
⎣1  1⎦ 

要创建对角矩阵，请使用diag。diag的参数可以是数字或矩阵。数字被解释为大小为(1\times 1)的矩阵。矩阵按对角线堆叠。剩余元素填充为(0)。

>>> diag(1, 2, 3)
⎡1  0  0⎤
⎢       ⎥
⎢0  2  0⎥
⎢       ⎥
⎣0  0  3⎦
>>> diag(-1, ones(2, 2), Matrix([5, 7, 5]))
⎡-1  0  0  0⎤
⎢           ⎥
⎢0   1  1  0⎥
⎢           ⎥
⎢0   1  1  0⎥
⎢           ⎥
⎢0   0  0  5⎥
⎢           ⎥
⎢0   0  0  7⎥
⎢           ⎥
⎣0   0  0  5⎦ 

高级方法

行列式

要计算矩阵的行列式，请使用det。

>>> M = Matrix([[1, 0, 1], [2, -1, 3], [4, 3, 2]])
>>> M
⎡1  0   1⎤
⎢        ⎥
⎢2  -1  3⎥
⎢        ⎥
⎣4  3   2⎦
>>> M.det()
-1 

行阶梯形式

要将矩阵转换为行阶梯形式，请使用rref。rref返回一个包含两个元素的元组。第一个元素是行阶梯形式，第二个是主元列的索引的元组。

>>> M = Matrix([[1, 0, 1, 3], [2, 3, 4, 7], [-1, -3, -3, -4]])
>>> M
⎡1   0   1   3 ⎤
⎢              ⎥
⎢2   3   4   7 ⎥
⎢              ⎥
⎣-1  -3  -3  -4⎦
>>> M.rref()
⎛⎡1  0   1    3 ⎤        ⎞
⎜⎢              ⎥        ⎟
⎜⎢0  1  2/3  1/3⎥, (0, 1)⎟
⎜⎢              ⎥        ⎟
⎝⎣0  0   0    0 ⎦        ⎠ 

注意

函数rref返回的元组的第一个元素是Matrix类型。第二个元素是tuple类型。

零空间

要找到矩阵的零空间，请使用nullspace。nullspace返回一个列向量列表，这些向量跨越矩阵的零空间。

>>> M = Matrix([[1, 2, 3, 0, 0], [4, 10, 0, 0, 1]])
>>> M
⎡1  2   3  0  0⎤
⎢              ⎥
⎣4  10  0  0  1⎦
>>> M.nullspace()
⎡⎡-15⎤  ⎡0⎤  ⎡ 1  ⎤⎤
⎢⎢   ⎥  ⎢ ⎥  ⎢    ⎥⎥
⎢⎢ 6 ⎥  ⎢0⎥  ⎢-1/2⎥⎥
⎢⎢   ⎥  ⎢ ⎥  ⎢    ⎥⎥
⎢⎢ 1 ⎥, ⎢0⎥, ⎢ 0  ⎥⎥
⎢⎢   ⎥  ⎢ ⎥  ⎢    ⎥⎥
⎢⎢ 0 ⎥  ⎢1⎥  ⎢ 0  ⎥⎥
⎢⎢   ⎥  ⎢ ⎥  ⎢    ⎥⎥
⎣⎣ 0 ⎦  ⎣0⎦  ⎣ 1  ⎦⎦ 

列空间

要找到矩阵的列空间，请使用columnspace。columnspace返回一个列向量列表，这些向量跨越矩阵的列空间。

>>> M = Matrix([[1, 1, 2], [2 ,1 , 3], [3 , 1, 4]])
>>> M
⎡1  1  2⎤
⎢       ⎥
⎢2  1  3⎥
⎢       ⎥
⎣3  1  4⎦
>>> M.columnspace()
⎡⎡1⎤  ⎡1⎤⎤
⎢⎢ ⎥  ⎢ ⎥⎥
⎢⎢2⎥, ⎢1⎥⎥
⎢⎢ ⎥  ⎢ ⎥⎥
⎣⎣3⎦  ⎣1⎦⎦ 

特征值、特征向量和对角化

要找到矩阵的特征值，请使用eigenvals。eigenvals返回一个字典，包含特征值: 代数重数对（类似于 roots 的输出）。

>>> M = Matrix([[3, -2,  4, -2], [5,  3, -3, -2], [5, -2,  2, -2], [5, -2, -3,  3]])
>>> M
⎡3  -2  4   -2⎤
⎢             ⎥
⎢5  3   -3  -2⎥
⎢             ⎥
⎢5  -2  2   -2⎥
⎢             ⎥
⎣5  -2  -3  3 ⎦
>>> M.eigenvals()
{-2: 1, 3: 1, 5: 2} 

这意味着M具有特征值-2、3 和 5，并且特征值-2 和 3 的代数重数为 1，特征值 5 的代数重数为 2。

要找到矩阵的特征向量，请使用eigenvects。eigenvects返回一个元组列表，形式为(特征值，代数重数，[特征向量])。

>>> M.eigenvects()
⎡⎛       ⎡⎡0⎤⎤⎞  ⎛      ⎡⎡1⎤⎤⎞  ⎛      ⎡⎡1⎤  ⎡0 ⎤⎤⎞⎤
⎢⎜       ⎢⎢ ⎥⎥⎟  ⎜      ⎢⎢ ⎥⎥⎟  ⎜      ⎢⎢ ⎥  ⎢  ⎥⎥⎟⎥
⎢⎜       ⎢⎢1⎥⎥⎟  ⎜      ⎢⎢1⎥⎥⎟  ⎜      ⎢⎢1⎥  ⎢-1⎥⎥⎟⎥
⎢⎜-2, 1, ⎢⎢ ⎥⎥⎟, ⎜3, 1, ⎢⎢ ⎥⎥⎟, ⎜5, 2, ⎢⎢ ⎥, ⎢  ⎥⎥⎟⎥
⎢⎜       ⎢⎢1⎥⎥⎟  ⎜      ⎢⎢1⎥⎥⎟  ⎜      ⎢⎢1⎥  ⎢0 ⎥⎥⎟⎥
⎢⎜       ⎢⎢ ⎥⎥⎟  ⎜      ⎢⎢ ⎥⎥⎟  ⎜      ⎢⎢ ⎥  ⎢  ⎥⎥⎟⎥
⎣⎝       ⎣⎣1⎦⎦⎠  ⎝      ⎣⎣1⎦⎦⎠  ⎝      ⎣⎣0⎦  ⎣1 ⎦⎦⎠⎦ 

这表明，例如，特征值 5 也具有几何重数 2，因为它有两个特征向量。由于所有特征值的代数和几何重数相同，因此M是可对角化的。

要对角化矩阵，请使用diagonalize。diagonalize返回一个元组((P, D))，其中(D)是对角化的，且(M = PDP^{-1})。

>>> P, D = M.diagonalize()
>>> P
⎡0  1  1  0 ⎤
⎢           ⎥
⎢1  1  1  -1⎥
⎢           ⎥
⎢1  1  1  0 ⎥
⎢           ⎥
⎣1  1  0  1 ⎦
>>> D
⎡-2  0  0  0⎤
⎢           ⎥
⎢0   3  0  0⎥
⎢           ⎥
⎢0   0  5  0⎥
⎢           ⎥
⎣0   0  0  5⎦
>>> P*D*P**-1
⎡3  -2  4   -2⎤
⎢             ⎥
⎢5  3   -3  -2⎥
⎢             ⎥
⎢5  -2  2   -2⎥
⎢             ⎥
⎣5  -2  -3  3 ⎦
>>> P*D*P**-1 == M
True 

注意，由于eigenvects也包含了特征值，如果你还想要特征向量，应该使用它而不是eigenvals。然而，由于计算特征向量可能非常耗时，如果只想找特征值，应优先选择eigenvals。

如果你只想得到特征多项式，请使用charpoly。这比使用eigenvals更有效率，因为有时符号根可能计算代价高昂。

>>> lamda = symbols('lamda')
>>> p = M.charpoly(lamda)
>>> factor(p.as_expr())
 2
(λ - 5) ⋅(λ - 3)⋅(λ + 2) 

可能出现的问题

零测试

如果你的矩阵操作失败或返回错误答案，常见原因可能是由于零测试不正确。如果表达式未经适当的零测试，可能会导致高斯消元中找不到主元，或者决定矩阵是否可逆，或者依赖先前过程的任何高级函数可能存在问题。

目前，SymPy 的默认零测试方法_iszero仅在某些有限的数值和符号域内保证准确性，对于其无法决策的复杂表达式，则被视为None，其行为类似逻辑False。

使用零测试过程的方法列表如下：

echelon_form，is_echelon，rank，rref，nullspace，eigenvects，inverse_ADJ，inverse_GE，inverse_LU，LUdecomposition，LUdecomposition_Simple，LUsolve

它们具有属性iszerofunc，供用户指定零测试方法，可以接受具有单一输入和布尔输出的任何函数，其默认值为_iszero。

这里是一个解决由于未经充分测试的零值引起的问题的示例。尽管这个特定矩阵的输出已经得到改进，但以下技术仍然具有一定的兴趣。[1] [2] [3]

>>> from sympy import *
>>> q = Symbol("q", positive = True)
>>> m = Matrix([
... [-2*cosh(q/3),      exp(-q),            1],
... [      exp(q), -2*cosh(q/3),            1],
... [           1,            1, -2*cosh(q/3)]])
>>> m.nullspace() 
[] 

您可以通过启用警告来追踪哪些表达式被低估了，通过注入自定义的零测试。

>>> import warnings
>>>
>>> def my_iszero(x):
...     result = x.is_zero
...
...     # Warnings if evaluated into None
...     if result is None:
...         warnings.warn("Zero testing of {} evaluated into None".format(x))
...     return result
...
>>> m.nullspace(iszerofunc=my_iszero) 
__main__:9: UserWarning: Zero testing of 4*cosh(q/3)**2 - 1 evaluated into None
__main__:9: UserWarning: Zero testing of (-exp(q) - 2*cosh(q/3))*(-2*cosh(q/3) - exp(-q)) - (4*cosh(q/3)**2 - 1)**2 evaluated into None
__main__:9: UserWarning: Zero testing of 2*exp(q)*cosh(q/3) - 16*cosh(q/3)**4 + 12*cosh(q/3)**2 + 2*exp(-q)*cosh(q/3) evaluated into None
__main__:9: UserWarning: Zero testing of -(4*cosh(q/3)**2 - 1)*exp(-q) - 2*cosh(q/3) - exp(-q) evaluated into None
[] 

在这种情况下，(-exp(q) - 2*cosh(q/3))*(-2*cosh(q/3) - exp(-q)) - (4*cosh(q/3)**2 - 1)**2应该得到零，但零测试未能捕捉到。这可能意味着应引入更强的零测试。对于这个特定的示例，重写为指数函数并应用简化将使零测试对双曲线函数更强，同时对其他多项式或超越函数无害。

>>> def my_iszero(x):
...     result = x.rewrite(exp).simplify().is_zero
...
...     # Warnings if evaluated into None
...     if result is None:
...         warnings.warn("Zero testing of {} evaluated into None".format(x))
...     return result
...
>>> m.nullspace(iszerofunc=my_iszero) 
__main__:9: UserWarning: Zero testing of -2*cosh(q/3) - exp(-q) evaluated into None
⎡⎡  ⎛   q         ⎛q⎞⎞  -q         2⎛q⎞    ⎤⎤
⎢⎢- ⎜- ℯ  - 2⋅cosh⎜─⎟⎟⋅ℯ   + 4⋅cosh ⎜─⎟ - 1⎥⎥
⎢⎢  ⎝             ⎝3⎠⎠              ⎝3⎠    ⎥⎥
⎢⎢─────────────────────────────────────────⎥⎥
⎢⎢          ⎛      2⎛q⎞    ⎞     ⎛q⎞       ⎥⎥
⎢⎢        2⋅⎜4⋅cosh ⎜─⎟ - 1⎟⋅cosh⎜─⎟       ⎥⎥
⎢⎢          ⎝       ⎝3⎠    ⎠     ⎝3⎠       ⎥⎥
⎢⎢                                         ⎥⎥
⎢⎢           ⎛   q         ⎛q⎞⎞            ⎥⎥
⎢⎢          -⎜- ℯ  - 2⋅cosh⎜─⎟⎟            ⎥⎥
⎢⎢           ⎝             ⎝3⎠⎠            ⎥⎥
⎢⎢          ────────────────────           ⎥⎥
⎢⎢                   2⎛q⎞                  ⎥⎥
⎢⎢             4⋅cosh ⎜─⎟ - 1              ⎥⎥
⎢⎢                    ⎝3⎠                  ⎥⎥
⎢⎢                                         ⎥⎥
⎣⎣                    1                    ⎦⎦ 

在注入备用零测试后，您可以清楚地看到nullspace返回了正确的结果。

请注意，这种方法仅适用于仅包含数值、双曲线和指数函数的某些矩阵情况。对于其他矩阵，应选择其特定领域的不同方法。

可能的建议是利用重写和简化的方式，以牺牲速度为代价[4]，或者使用随机数测试的方式，以牺牲精确度为代价[5]。

如果您想知道为什么没有通用的零测试算法可以与任何符号实体一起使用，那是因为存在零测试不可判定的常数问题[6]，而不仅仅是 SymPy，其他计算代数系统[7] [8]也会面临同样的根本性问题。

然而，发现任何零测试失败的案例可以提供一些优化 SymPy 的良好例子，因此如果您遇到了其中一个问题，可以将问题报告给 SymPy 问题跟踪器[9]，以获取社区的详细帮助。

脚注

高级表达式操作

原文：docs.sympy.org/latest/tutorials/intro-tutorial/manipulation.html

在本节中，我们讨论了一些进行高级表达式操作的方法。

理解表达式树

在我们进行此操作之前，我们需要了解 SymPy 中表达式的表示方法。数学表达式被表示为一棵树。让我们看一下表达式 (x² + xy)，即 x**2 + x*y。我们可以使用 srepr 看到这个表达式的内部结构。

>>> from sympy import *
>>> x, y, z = symbols('x y z') 

>>> expr = x**2 + x*y
>>> srepr(expr)
"Add(Pow(Symbol('x'), Integer(2)), Mul(Symbol('x'), Symbol('y')))" 

最简单的方法是查看表达式树的图示：

digraph{ # Graph style "ordering"="out" "rankdir"="TD" ######### # Nodes # ######### "Add(Pow(Symbol('x'), Integer(2)), Mul(Symbol('x'), Symbol('y')))()" ["color"="black", "label"="Add", "shape"="ellipse"]; "Pow(Symbol('x'), Integer(2))(0,)" ["color"="black", "label"="Pow", "shape"="ellipse"]; "Symbol('x')(0, 0)" ["color"="black", "label"="Symbol('x')", "shape"="ellipse"]; "Integer(2)(0, 1)" ["color"="black", "label"="Integer(2)", "shape"="ellipse"]; "Mul(Symbol('x'), Symbol('y'))(1,)" ["color"="black", "label"="Mul", "shape"="ellipse"]; "Symbol('x')(1, 0)" ["color"="black", "label"="Symbol('x')", "shape"="ellipse"]; "Symbol('y')(1, 1)" ["color"="black", "label"="Symbol('y')", "shape"="ellipse"]; ######### # Edges # ######### "Add(Pow(Symbol('x'), Integer(2)), Mul(Symbol('x'), Symbol('y')))()" -> "Pow(Symbol('x'), Integer(2))(0,)"; "Add(Pow(Symbol('x'), Integer(2)), Mul(Symbol('x'), Symbol('y')))()" -> "Mul(Symbol('x'), Symbol('y'))(1,)"; "Pow(Symbol('x'), Integer(2))(0,)" -> "Symbol('x')(0, 0)"; "Pow(Symbol('x'), Integer(2))(0,)" -> "Integer(2)(0, 1)"; "Mul(Symbol('x'), Symbol('y'))(1,)" -> "Symbol('x')(1, 0)"; "Mul(Symbol('x'), Symbol('y'))(1,)" -> "Symbol('y')_(1, 1)"; }

注意

上述图表是使用 Graphviz 和 dotprint 函数创建的。

首先，让我们看一下这棵树的叶子节点。符号是类 Symbol 的实例。虽然我们一直在做

>>> x = symbols('x') 

我们也可以这样做

>>> x = Symbol('x') 

无论哪种方式，我们都会得到一个名为“x”的符号 [1]。在表达式中的数字 2，我们得到了 Integer(2)。Integer 是 SymPy 中整数的类。它类似于 Python 内置类型 int，不过 Integer 与其他 SymPy 类型协作更好。

当我们写 x**2 时，这将创建一个 Pow 对象。Pow 是“power”的缩写。

>>> srepr(x**2)
"Pow(Symbol('x'), Integer(2))" 

我们可以通过调用 Pow(x, 2) 创建相同的对象。

>>> Pow(x, 2)
x**2 

请注意，在 srepr 输出中，我们看到 Integer(2)，这是 SymPy 版本的整数，尽管从技术上讲，我们输入了 Python 的 int 类型的 2。通常情况下，当您通过某个函数或操作将 SymPy 对象与非 SymPy 对象组合时，非 SymPy 对象将被转换为 SymPy 对象。执行此操作的函数是 sympify [2]。

>>> type(2)
<... 'int'>
>>> type(sympify(2))
<class 'sympy.core.numbers.Integer'> 

我们已经看到 x**2 表示为 Pow(x, 2)。那么 x*y 呢？正如我们所预期的那样，这是 x 和 y 的乘积。SymPy 中用于乘法的类是 Mul。

>>> srepr(x*y)
"Mul(Symbol('x'), Symbol('y'))" 

因此，我们可以通过编写 Mul(x, y) 来创建相同的对象。

>>> Mul(x, y)
x*y 

现在我们来到我们的最终表达式，x**2 + x*y。这是我们最后两个对象 Pow(x, 2) 和 Mul(x, y) 的加法。SymPy 中用于加法的类是 Add，因此，正如你所预期的那样，要创建这个对象，我们使用 Add(Pow(x, 2), Mul(x, y))。

>>> Add(Pow(x, 2), Mul(x, y))
x**2 + x*y 

SymPy 表达式树可以有许多分支，可以非常深或非常宽。这里是一个更复杂的例子。

>>> expr = sin(x*y)/2 - x**2 + 1/y
>>> srepr(expr)
"Add(Mul(Integer(-1), Pow(Symbol('x'), Integer(2))), Mul(Rational(1, 2),
sin(Mul(Symbol('x'), Symbol('y')))), Pow(Symbol('y'), Integer(-1)))" 

这里是一个图示。

digraph{ # 图形样式 "rankdir"="TD" ######### # 节点 # ######### "Half()(0, 0)" ["color"="black", "label"="有理数(1, 2)", "shape"="ellipse"]; "Symbol(y)(2, 0)" ["color"="black", "label"="符号('y')", "shape"="ellipse"]; "Symbol(x)(1, 1, 0)" ["color"="black", "label"="符号('x')", "shape"="ellipse"]; "Integer(2)(1, 1, 1)" ["color"="black", "label"="整数(2)", "shape"="ellipse"]; "NegativeOne()(2, 1)" ["color"="black", "label"="整数(-1)", "shape"="ellipse"]; "NegativeOne()(1, 0)" ["color"="black", "label"="整数(-1)", "shape"="ellipse"]; "Symbol(y)(0, 1, 0, 1)" ["color"="black", "label"="符号('y')", "shape"="ellipse"]; "Symbol(x)(0, 1, 0, 0)" ["color"="black", "label"="符号('x')", "shape"="ellipse"]; "Pow(Symbol(x), Integer(2))(1, 1)" ["color"="black", "label"="Pow", "shape"="ellipse"]; "Pow(Symbol(y), NegativeOne())(2,)" ["color"="black", "label"="Pow", "shape"="ellipse"]; "Mul(Symbol(x), Symbol(y))(0, 1, 0)" ["color"="black", "label"="Mul", "shape"="ellipse"]; "sin(Mul(Symbol(x), Symbol(y)))(0, 1)" ["color"="black", "label"="sin", "shape"="ellipse"]; "Mul(Half(), sin(Mul(Symbol(x), Symbol(y))))(0,)" ["color"="black", "label"="Mul", "shape"="ellipse"]; "Mul(NegativeOne(), Pow(Symbol(x), Integer(2)))(1,)" ["color"="black", "label"="Mul", "shape"="ellipse"]; "Add(Mul(Half(), sin(Mul(Symbol(x), Symbol(y)))), Mul(NegativeOne(), Pow(Symbol(x), Integer(2))), Pow(Symbol(y), NegativeOne()))()" ["color"="black", "label"="Add", "shape"="ellipse"]; ######### # 边缘 # ######### "Pow(Symbol(y), NegativeOne())(2,)" -> "Symbol(y)(2, 0)"; "Pow(Symbol(x), Integer(2))(1, 1)" -> "Symbol(x)(1, 1, 0)"; "Pow(Symbol(x), Integer(2))(1, 1)" -> "Integer(2)(1, 1, 1)"; "Pow(Symbol(y), NegativeOne())(2,)" -> "NegativeOne()(2, 1)"; "Mul(Symbol(x), Symbol(y))(0, 1, 0)" -> "Symbol(x)(0, 1, 0, 0)"; "Mul(Symbol(x), Symbol(y))(0, 1, 0)" -> "Symbol(y)(0, 1, 0, 1)"; "Mul(Half(), sin(Mul(Symbol(x), Symbol(y))))(0,)" -> "Half()(0, 0)"; "Mul(NegativeOne(), Pow(Symbol(x), Integer(2)))(1,)" -> "NegativeOne()(1, 0)"; "sin(Mul(Symbol(x), Symbol(y)))(0, 1)" -> "Mul(Symbol(x), Symbol(y))(0, 1, 0)"; "Mul(NegativeOne(), Pow(Symbol(x), Integer(2)))(1,)" -> "Pow(Symbol(x), Integer(2))(1, 1)"; "Mul(Half(), sin(Mul(Symbol(x), Symbol(y))))(0,)" -> "sin(Mul(Symbol(x), Symbol(y)))(0, 1)"; "Add(Mul(Half(), sin(Mul(Symbol(x), Symbol(y)))), Mul(NegativeOne(), Pow(Symbol(x), Integer(2))), Pow(Symbol(y), NegativeOne()))()" -> "Pow(Symbol(y), NegativeOne())(2,)"; "Add(Mul(Half(), sin(Mul(Symbol(x), Symbol(y)))), Mul(NegativeOne(), Pow(Symbol(x), Integer(2))), Pow(Symbol(y), NegativeOne()))()" -> "Mul(Half(), sin(Mul(Symbol(x), Symbol(y))))(0,)"; "Add(Mul(Half(), sin(Mul(Symbol(x), Symbol(y)))), Mul(NegativeOne(), Pow(Symbol(x), Integer(2))), Pow(Symbol(y), NegativeOne()))()" -> "Mul(NegativeOne(), Pow(Symbol(x), Integer(2)))_(1,)"; }

这个表达式揭示了一些关于 SymPy 表达树的有趣事情。让我们逐一了解它们。

• 让我们首先看看 x**2 项。正如我们预期的那样，我们看到 Pow(x, 2)。再上一层，我们有 Mul(-1, Pow(x, 2))。在 SymPy 中没有减法类。x - y 被表示为 x + -y，或者更完整地说，x + -1*y，即 Add(x, Mul(-1, y))。

>>> srepr(x - y)
"Add(Symbol('x'), Mul(Integer(-1), Symbol('y')))" 

digraph{ # 图形样式 "rankdir"="TD" ######### # 节点 # ######### "Symbol(x)(1,)" ["color"="black", "label"="Symbol('x')", "shape"="ellipse"]; "Symbol(y)(0, 1)" ["color"="black", "label"="Symbol('y')", "shape"="ellipse"]; "NegativeOne()(0, 0)" ["color"="black", "label"="Integer(-1)", "shape"="ellipse"]; "Mul(NegativeOne(), Symbol(y))(0,)" ["color"="black", "label"="Mul", "shape"="ellipse"]; "Add(Mul(NegativeOne(), Symbol(y)), Symbol(x))()" ["color"="black", "label"="Add", "shape"="ellipse"]; ######### # 边 # ######### "Mul(NegativeOne(), Symbol(y))(0,)" -> "Symbol(y)(0, 1)"; "Mul(NegativeOne(), Symbol(y))(0,)" -> "NegativeOne()(0, 0)"; "Add(Mul(NegativeOne(), Symbol(y)), Symbol(x))()" -> "Symbol(x)(1,)"; "Add(Mul(NegativeOne(), Symbol(y)), Symbol(x))()" -> "Mul(NegativeOne(), Symbol(y))_(0,)"; }

• 接下来，看看 1/y。我们可能期望看到类似 Div(1, y) 的东西，但类似于减法，在 SymPy 中没有除法类。相反，除法被表示为 -1 的幂。因此，我们有 Pow(y, -1)。如果我们将其他东西除以 y 而不是 1，例如 x/y，让我们看看。

>>> expr = x/y
>>> srepr(expr)
"Mul(Symbol('x'), Pow(Symbol('y'), Integer(-1)))" 

digraph{ # 图形样式 "rankdir"="TD" ######### # 节点 # ######### "Symbol(x)(0,)" ["color"="black", "label"="Symbol('x')", "shape"="ellipse"]; "Symbol(y)(1, 0)" ["color"="black", "label"="Symbol('y')", "shape"="ellipse"]; "NegativeOne()(1, 1)" ["color"="black", "label"="Integer(-1)", "shape"="ellipse"]; "Pow(Symbol(y), NegativeOne())(1,)" ["color"="black", "label"="Pow", "shape"="ellipse"]; "Mul(Symbol(x), Pow(Symbol(y), NegativeOne()))()" ["color"="black", "label"="Mul", "shape"="ellipse"]; ######### # 边 # ######### "Pow(Symbol(y), NegativeOne())(1,)" -> "Symbol(y)(1, 0)"; "Pow(Symbol(y), NegativeOne())(1,)" -> "NegativeOne()(1, 1)"; "Mul(Symbol(x), Pow(Symbol(y), NegativeOne()))()" -> "Symbol(x)(0,)"; "Mul(Symbol(x), Pow(Symbol(y), NegativeOne()))()" -> "Pow(Symbol(y), NegativeOne())_(1,)"; }

• 我们看到 x/y 被表示为 x*y**-1，即 Mul(x, Pow(y, -1))。

• 最后，让我们看看 sin(x*y)/2 项。按照前面示例的模式，我们可能期望看到 Mul(sin(x*y), Pow(Integer(2), -1))。但实际上，我们看到的是 Mul(Rational(1, 2), sin(x*y))。有理数总是组合成一个乘法项，因此当我们除以 2 时，表示为乘以 1/2。

• 最后，还有一点要注意。您可能已经注意到，我们输入表达式的顺序和从 srepr 或图形中得到的顺序不同。您可能也在本教程的早些时候注意到了这种现象。例如

>>> 1 + x
x + 1 

这是因为在 SymPy 中，交换操作 Add 和 Mul 的参数存储在任意（但一致！）的顺序中，这与输入的顺序无关（如果您担心非交换乘法，请放心。在 SymPy 中，您可以使用 Symbol('A', commutative=False) 创建非交换符号，并且非交换符号的乘法顺序与输入保持一致）。此外，正如我们将在下一节看到的那样，打印顺序和内部存储顺序也可能不同。

通常，在使用 SymPy 表达式树时需要记住的一件重要事情是：表达式的内部表示和打印方式可能不同。输入形式也是如此。如果某些表达式操作算法的工作方式与您预期的不同，很可能是对象的内部表示与您想象的不同。

透过表达式树进行递归

现在您知道了 SymPy 中表达式树的工作方式，让我们看看如何通过表达式树深入了解它。SymPy 中的每个对象都有两个非常重要的属性，func 和 args。

func

func 是对象的头部。例如，(x*y).func 是 Mul。通常它与对象的类相同（尽管有例外）。

有关 func 的两个注意事项。首先，对象的类不一定与用于创建它的类相同。例如

>>> expr = Add(x, x)
>>> expr.func
<class 'sympy.core.mul.Mul'> 

我们创建了 Add(x, x)，所以我们可能期望 expr.func 是 Add，但实际上我们得到的是 Mul。为什么呢？让我们仔细看一下 expr。

>>> expr
2*x 

Add(x, x)，即 x + x，自动转换为 Mul(2, x)，即 2*x，这是一个 Mul。SymPy 类大量使用 __new__ 类构造函数，与 __init__ 不同，它允许从构造函数返回不同的类。

其次，一些类是特例，通常出于效率原因[3]。

>>> Integer(2).func
<class 'sympy.core.numbers.Integer'>
>>> Integer(0).func
<class 'sympy.core.numbers.Zero'>
>>> Integer(-1).func
<class 'sympy.core.numbers.NegativeOne'> 

大多数情况下，这些问题不会困扰我们。特殊类 Zero、One、NegativeOne 等都是 Integer 的子类，因此只要使用 isinstance，这不会成为问题。

args

args 是对象的顶层参数。(x*y).args 将是 (x, y)。让我们看一些例子

>>> expr = 3*y**2*x
>>> expr.func
<class 'sympy.core.mul.Mul'>
>>> expr.args
(3, x, y**2) 

从这里，我们可以看到 expr == Mul(3, y**2, x)。事实上，我们可以完全通过其 func 和 args 重新构建 expr。

>>> expr.func(*expr.args)
3*x*y**2
>>> expr == expr.func(*expr.args)
True 

注意虽然我们输入了 3*y**2*x，但 args 是 (3, x, y**2)。在 Mul 中，有理数系数将首先出现在 args 中，但除此之外，其他所有顺序都没有特殊模式。但可以肯定的是，有一个顺序。

>>> expr = y**2*3*x
>>> expr.args
(3, x, y**2) 

Mul 的 args 是排序的，因此相同的 Mul 将具有相同的 args。但是排序是基于一些旨在使排序唯一和有效的标准，没有数学意义。

我们的 expr 的 srepr 形式是 Mul(3, x, Pow(y, 2))。如果我们想要获取 Pow(y, 2) 的 args，请注意 y**2 在 expr.args 的第三个位置，即 expr.args[2]。

>>> expr.args[2]
y**2 

因此，要获取这个的 args，我们调用 expr.args[2].args。

>>> expr.args[2].args
(y, 2) 

现在如果我们尝试更深入地查看。y 的参数是什么。或者 2 的。我们来看看。

>>> y.args
()
>>> Integer(2).args
() 

他们两者都具有空的 args。在 SymPy 中，空的 args 表示我们已经到达了表达式树的叶子。

因此，SymPy 表达式有两种可能性。要么它具有空的 args，在这种情况下，它是任何表达式树中的叶子节点，要么它具有 args，在这种情况下，它是任何表达式树中的分支节点。当它具有 args 时，可以完全从其 func 和 args 重建它。这体现了关键不变量。

（回想一下，在 Python 中，如果 a 是一个元组，那么 f(*a) 表示用元组 a 中的元素调用 f，例如，f(*(1, 2, 3)) 等同于 f(1, 2, 3)。）

这一关键不变量使我们能够编写简单的算法来遍历表达式树，修改它们，并将它们重建为新的表达式。

遍历树

有了这些知识，让我们看看如何通过表达式树进行递归。args 的嵌套特性非常适合递归函数。基本情况将是空的 args。让我们编写一个简单的函数，它可以遍历表达式并在每个级别打印所有的 args。

>>> def pre(expr):
...     print(expr)
...     for arg in expr.args:
...         pre(arg) 

看到 () 如何在表达式树中表示叶子节点，我们甚至不必为递归编写基本情况；它会被 for 循环自动处理。

让我们测试我们的函数。

>>> expr = x*y + 1
>>> pre(expr)
x*y + 1
1
x*y
x
y 

你能猜到我们为什么称呼我们的函数为 pre 吗？我们刚刚为我们的表达式树写了一个前序遍历函数。看看你能否编写一个后序遍历函数。

在 SymPy 中，这种遍历非常常见，提供了生成器函数 preorder_traversal 和 postorder_traversal 来简化这种遍历过程。我们也可以将我们的算法编写为

>>> for arg in preorder_traversal(expr):
...     print(arg)
x*y + 1
1
x*y
x
y 

防止表达式求值

通常有两种方法可以防止表达式求值，一种是在构建表达式时传递 evaluate=False 参数，另一种是通过将表达式包装在 UnevaluatedExpr 中创建一个停止求值。

例如：

>>> from sympy import Add
>>> from sympy.abc import x, y, z
>>> x + x
2*x
>>> Add(x, x)
2*x
>>> Add(x, x, evaluate=False)
x + x 

如果您不记得要构建的表达式对应的类（通常假设 evaluate=True），只需使用 sympify 并传递一个字符串：

>>> from sympy import sympify
>>> sympify("x + x", evaluate=False)
x + x 

注意，evaluate=False 不会防止在后续使用表达式时进行求值：

>>> expr = Add(x, x, evaluate=False)
>>> expr
x + x
>>> expr + x
3*x 

这就是为什么 UnevaluatedExpr 类很方便。UnevaluatedExpr 是 SymPy 提供的一种方法，允许用户保持表达式未求值。通过 未求值 意味着其中的值不会与外部表达式交互以提供简化的输出。例如：

>>> from sympy import UnevaluatedExpr
>>> expr = x + UnevaluatedExpr(x)
>>> expr
x + x
>>> x + expr
2*x + x 

保持独立的 x 是由 UnevaluatedExpr 包裹的 x。要释放它：

>>> (x + expr).doit()
3*x 

其他例子：

>>> from sympy import *
>>> from sympy.abc import x, y, z
>>> uexpr = UnevaluatedExpr(S.One*5/7)*UnevaluatedExpr(S.One*3/4)
>>> uexpr
(5/7)*(3/4)
>>> x*UnevaluatedExpr(1/x)
x*1/x 

值得注意的是，UnevaluatedExpr 无法阻止作为参数给出的表达式的评估。例如：

>>> expr1 = UnevaluatedExpr(x + x)
>>> expr1
2*x
>>> expr2 = sympify('x + x', evaluate=False)
>>> expr2
x + x 

记住，如果将 expr2 包含到另一个表达式中，它将被评估。结合这两种方法可以同时阻止内部和外部的评估：

>>> UnevaluatedExpr(sympify("x + x", evaluate=False)) + y
y + (x + x) 

UnevaluatedExpr 受 SymPy 打印机支持，并可用于以不同的输出形式打印结果。例如

>>> from sympy import latex
>>> uexpr = UnevaluatedExpr(S.One*5/7)*UnevaluatedExpr(S.One*3/4)
>>> print(latex(uexpr))
\frac{5}{7} \cdot \frac{3}{4} 

要释放表达式并获得评估后的 LaTeX 形式，只需使用 .doit()：

>>> print(latex(uexpr.doit()))
\frac{15}{28} 

脚注

接下来

原文：docs.sympy.org/latest/tutorials/intro-tutorial/next.html

恭喜您完成 SymPy 教程！

如果您是开发者，有兴趣在代码中使用 SymPy，请访问使用指南，这些指南讨论了关键的开发任务。

中级 SymPy 用户和开发者可能希望访问解释部分了解常见问题和高级主题。

SymPy API 参考详细描述了 SymPy API。

如果您有兴趣为 SymPy 做出贡献，请访问贡献指南。

物理教程

原文：docs.sympy.org/latest/tutorials/physics/index.html

物理教程旨在向尚未使用 SymPy 物理功能的用户介绍 SymPy。这里展示的功能比教程涵盖的要多得多。通过深入的示例和练习，教程演示了如何利用 SymPy 的各种物理子包，包括但不限于多体动力学、量子力学、光学和连续介质力学，来建模和模拟不同的物理系统及其行为。

对本教程或 SymPy 总体的反馈，我们随时欢迎。只需写信到我们的邮件列表。

内容

• 生物力学教程

  • 生物力学建模导论

  • 生物力学模型示例

生物力学教程

原文：docs.sympy.org/latest/tutorials/physics/biomechanics/index.html

这些教程提供了使用 SymPy 进行生物力学模拟和分析的全面指南。我们涵盖了各种模型，包括人臂移动杠杆、肌肉和肌腱使用希尔型肌肉模型产生的力。

在这些教程中，您可以期待：

• 模型描述：详细解释每个生物力学模型。

• 变量和运动学：定义用于建模的必要变量和运动学方程。

• 建模：逐步构建生物力学模型的过程。

• 运动方程：系统运动方程的推导和分析。

• 生物力学建模介绍

  • 载荷

  • 路径

  • 包裹几何体

  • 执行器

  • 激活动力学

  • 肌腱曲线

  • 肌腱动力学

  • 一个简单的肌腱模型

  • 参考文献

• 生物力学模型示例

  • 模型描述

  • 定义变量

  • 定义运动学

  • 定义惯性

  • 定义力

  • 运动方程

  • 肌肉激活微分方程

  • 评估系统微分方程

  • 模拟肌肉驱动运动

  • 结论

  • 参考文献

生物力学建模简介

原文：docs.sympy.org/latest/tutorials/physics/biomechanics/biomechanics.html

sympy.physics.biomechanics 提供了增强使用 sympy.physics.mechanics 创建的模型的功能，这些模型包括用于建模肌肉和肌腱的力元素。在本教程中，我们将介绍此模块的特性。

生物力学包的初始主要目的是介绍用于建模由Hill 型肌肉模型产生的力的工具。这些模型生成基于肌肉的收缩状态和肌腱的被动伸展作用于生物体的骨骼结构上的力。在本教程中，我们介绍组成肌腱模型的元素，然后展示其在具体实现中的运行，即MusculotendonDeGroote2016 模型。

载荷

sympy.physics.mechanics 包括两种类型的载荷：Force 和 Torque。力表示沿着作用线方向的绑定向量量，而力矩是未绑定的向量，表示一对力的结果力矩。

一个非常常见的力模型示例是平行线性弹簧和线性阻尼器。质量为 (m) 的粒子，其由广义坐标 (x(t)) 描述的一维运动，具有线性弹簧和阻尼器系数 (k) 和 (c)，其运动方程如下所示：

[m \ddot{x} = \sum F = -kx - c\dot{x}]

在 SymPy 中，我们可以这样表达作用在粒子 (P) 上的力，该粒子在参考系 (N) 中运动，并且相对于固定在 (N) 中的点 (O) 的位置如下：

>>> import sympy as sm
>>> import sympy.physics.mechanics as me 

>>> k, c = sm.symbols('k, c', real=True, nonnegative=True)
>>> x = me.dynamicsymbols('x', real=True) 

>>> N = me.ReferenceFrame('N')
>>> O, P = me.Point('O'), me.Point('P') 

>>> P.set_pos(O, x*N.x)
>>> P.set_vel(N, x.diff()*N.x) 

>>> force_on_P = me.Force(P, -k*P.pos_from(O) - c*P.vel(N))
>>> force_on_P
(P, (-c*Derivative(x(t), t) - k*x(t))*N.x) 

并且在 (O) 上会有一个等大而相反方向的力作用：

>>> force_on_O = me.Force(O, k*P.pos_from(O) + c*P.vel(N))
>>> force_on_O
(O, (c*Derivative(x(t), t) + k*x(t))*N.x) 

由单一肌肉和肌腱施加到一组刚体上的力将是进一步开发的肌肉肌腱模型的主要输出。

Pathways

肌肉及其相关的肌腱包裹在运动的骨骼系统周围，以及其他肌肉和器官。这带来了一个挑战，即确定肌肉和肌腱在接触到骨骼和器官时产生的力的作用线。我们引入了 pathway 模块来帮助管理几何关系到力的作用线的规范化。

上述弹簧阻尼器示例具有最简单的作用线定义，因此我们可以使用 LinearPathway 来建立该作用线。首先提供力将向两端点施加等大小但方向相反的应用点，路径的长度和两点之间的距离和两点之间的相对速度由 length 和 extension_velocity 路径计算。请注意，正速度意味着点正在远离彼此。还请注意，该公式处理了 (x) 是正数或负数的情况。

>>> lpathway = me.LinearPathway(O, P)
>>> lpathway
LinearPathway(O, P)
>>> lpathway.length
Abs(x(t))
>>> lpathway.extension_velocity
sign(x(t))*Derivative(x(t), t) 

然后，to_loads 方法会采用一种符号约定的力量大小，其中正的大小会推动两个点远离彼此，并返回作用于这两个点的所有力的列表。

>>> import pprint
>>> pprint.pprint(lpathway.to_loads(-k*x - k*x.diff()))
[Force(point=O, force=(k*x(t) + k*Derivative(x(t), t))*x(t)/Abs(x(t))*N.x),
 Force(point=P, force=(-k*x(t) - k*Derivative(x(t), t))*x(t)/Abs(x(t))*N.x)] 

可以用任意几何形状和任意数量的相互连接的粒子和刚体构造路径。例如，更复杂的路径是 ObstacleSetPathway。您可以指定两个路径端点之间的任意数量的中间点，力的驱动路径会沿着这些中间点连接 (N) 到 (O) 到 (Q) 到 (R) 到 (P)。这四个点各自会经历结果力。为简单起见，我们展示了仅弹簧力的效果。

>>> Q, R = me.Point('Q'), me.Point('R')
>>> Q.set_pos(O, 1*N.y)
>>> R.set_pos(O, 1*N.x + 1*N.y)
>>> opathway = me.ObstacleSetPathway(O, Q, R, P)
>>> opathway.length
sqrt((x(t) - 1)**2 + 1) + 2
>>> opathway.extension_velocity
(x(t) - 1)*Derivative(x(t), t)/sqrt((x(t) - 1)**2 + 1)
>>> pprint.pprint(opathway.to_loads(-k*opathway.length))
[Force(point=O, force=k*(sqrt((x(t) - 1)**2 + 1) + 2)*N.y),
 Force(point=Q, force=- k*(sqrt((x(t) - 1)**2 + 1) + 2)*N.y),
 Force(point=Q, force=k*(sqrt((x(t) - 1)**2 + 1) + 2)*N.x),
 Force(point=R, force=- k*(sqrt((x(t) - 1)**2 + 1) + 2)*N.x),
 Force(point=R, force=k*(sqrt((x(t) - 1)**2 + 1) + 2)*(x(t) - 1)/sqrt((x(t) - 1)**2 + 1)*N.x - k*(sqrt((x(t) - 1)**2 + 1) + 2)/sqrt((x(t) - 1)**2 + 1)*N.y),
 Force(point=P, force=- k*(sqrt((x(t) - 1)**2 + 1) + 2)*(x(t) - 1)/sqrt((x(t) - 1)**2 + 1)*N.x + k*(sqrt((x(t) - 1)**2 + 1) + 2)/sqrt((x(t) - 1)**2 + 1)*N.y)] 

如果设定 (x=1)，就更容易看到力的集合是正确的：

>>> for load in opathway.to_loads(-k*opathway.length):
...     pprint.pprint(me.Force(load[0], load[1].subs({x: 1})))
Force(point=O, force=3*k*N.y)
Force(point=Q, force=- 3*k*N.y)
Force(point=Q, force=3*k*N.x)
Force(point=R, force=- 3*k*N.x)
Force(point=R, force=- 3*k*N.y)
Force(point=P, force=3*k*N.y) 

你可以通过子类化PathwayBase来创建自己的路径。

包裹几何体

肌肉通常包裹在骨骼、组织或器官上。我们引入了包裹几何体和相关的包裹路径来帮助管理它们的复杂性。例如，如果两个路径端点位于圆柱体表面上，则力沿着连接端点的测地线的切线方向作用。WrappingCylinder对象计算了路径的复杂几何形状。然后，WrappingPathway利用这些几何形状构建力。沿着这条路径的弹簧力如下构建：

>>> r = sm.symbols('r', real=True, nonegative=True)
>>> theta = me.dynamicsymbols('theta', real=True)
>>> O, P, Q = sm.symbols('O, P, Q', cls=me.Point)
>>> A = me.ReferenceFrame('A') 

>>> A.orient_axis(N, theta, N.z) 

>>> P.set_pos(O, r*N.x)
>>> Q.set_pos(O, N.z + r*A.x) 

>>> cyl = me.WrappingCylinder(r, O, N.z)
>>> wpathway = me.WrappingPathway(P, Q, cyl)
>>> pprint.pprint(wpathway.to_loads(-k*wpathway.length))
[Force(point=P, force=- k*r*Abs(theta(t))*N.y - k*N.z),
 Force(point=Q, force=k*N.z + k*r*Abs(theta(t))*A.y),
 Force(point=O, force=k*r*Abs(theta(t))*N.y - k*r*Abs(theta(t))*A.y)] 

执行器

多体系统的模型通常具有时间变化的输入，如力或扭矩的大小。在许多情况下，这些指定的输入可以来自系统的状态甚至来自另一个动态系统的输出。actuator模块包含帮助管理这些力和扭矩输入模型的类。执行器旨在表示真实的物理组件。例如，上述弹簧-阻尼力可以通过子类化ActuatorBase并实现生成与该弹簧-阻尼执行器相关的负载的方法来创建。

>>> N = me.ReferenceFrame('N')
>>> O, P = me.Point('O'), me.Point('P')
>>> P.set_pos(O, x*N.x) 

>>> class SpringDamper(me.ActuatorBase):
...
...     # positive x spring is in tension
...     # negative x spring is in compression
...     def __init__(self, P1, P2, spring_constant, damper_constant):
...         self.P1 = P1
...         self.P2 = P2
...         self.k = spring_constant
...         self.c = damper_constant
...
...     def to_loads(self):
...         x = self.P2.pos_from(self.P1).magnitude()
...         v = x.diff(me.dynamicsymbols._t)
...         dir_vec = self.P2.pos_from(self.P1).normalize()
...         force_P1 = me.Force(self.P1,
...                             self.k*x*dir_vec + self.c*v*dir_vec)
...         force_P2 = me.Force(self.P2,
...                             -self.k*x*dir_vec - self.c*v*dir_vec)
...         return [force_P1, force_P2]
... 

>>> spring_damper = SpringDamper(O, P, k, c)
>>> pprint.pprint(spring_damper.to_loads())
[Force(point=O, force=(c*x(t)*sign(x(t))*Derivative(x(t), t)/Abs(x(t)) + k*x(t))*N.x),
 Force(point=P, force=(-c*x(t)*sign(x(t))*Derivative(x(t), t)/Abs(x(t)) - k*x(t))*N.x)] 

还有一个ForceActuator，可以与路径对象无缝集成。你只需在子类的初始化中设置.force属性即可。

>>> class SpringDamper(me.ForceActuator):
...
...     # positive x spring is in tension
...     # negative x spring is in compression
...     def __init__(self, pathway, spring_constant, damper_constant):
...         self.pathway = pathway
...         self.force = (-spring_constant*pathway.length -
...                       damper_constant*pathway.extension_velocity)
...
>>> spring_damper2 = SpringDamper(lpathway, k, c)
>>> pprint.pprint(spring_damper2.to_loads())
[Force(point=O, force=(c*sign(x(t))*Derivative(x(t), t) + k*Abs(x(t)))*x(t)/Abs(x(t))*N.x),
 Force(point=P, force=(-c*sign(x(t))*Derivative(x(t), t) - k*Abs(x(t)))*x(t)/Abs(x(t))*N.x)] 

这样就可以轻松地对其他路径应用弹簧-阻尼力，例如：

>>> spring_damper3 = SpringDamper(wpathway, k, c)
>>> pprint.pprint(spring_damper3.to_loads())
[Force(point=P, force=r*(-c*r**2*theta(t)*Derivative(theta(t), t)/sqrt(r**2*theta(t)**2 + 1) - k*sqrt(r**2*theta(t)**2 + 1))*Abs(theta(t))/sqrt(r**2*theta(t)**2 + 1)*N.y + (-c*r**2*theta(t)*Derivative(theta(t), t)/sqrt(r**2*theta(t)**2 + 1) - k*sqrt(r**2*theta(t)**2 + 1))/sqrt(r**2*theta(t)**2 + 1)*N.z),
 Force(point=Q, force=- (-c*r**2*theta(t)*Derivative(theta(t), t)/sqrt(r**2*theta(t)**2 + 1) - k*sqrt(r**2*theta(t)**2 + 1))/sqrt(r**2*theta(t)**2 + 1)*N.z - r*(-c*r**2*theta(t)*Derivative(theta(t), t)/sqrt(r**2*theta(t)**2 + 1) - k*sqrt(r**2*theta(t)**2 + 1))*Abs(theta(t))/sqrt(r**2*theta(t)**2 + 1)*A.y),
 Force(point=O, force=- r*(-c*r**2*theta(t)*Derivative(theta(t), t)/sqrt(r**2*theta(t)**2 + 1) - k*sqrt(r**2*theta(t)**2 + 1))*Abs(theta(t))/sqrt(r**2*theta(t)**2 + 1)*N.y + r*(-c*r**2*theta(t)*Derivative(theta(t), t)/sqrt(r**2*theta(t)**2 + 1) - k*sqrt(r**2*theta(t)**2 + 1))*Abs(theta(t))/sqrt(r**2*theta(t)**2 + 1)*A.y)] 

激活动力学

当肌腱模型被激活时，它们能够产生主动收缩力。在生物学上，当肌纤维中的(\textrm{Ca}^{2+})离子浓度足够高时，它们开始自发收缩。这种自发收缩状态称为“激活”。在生物力学模型中，通常用符号(a(t))表示，它被视为在([0, 1])范围内的归一化量。

生物体并不直接控制其肌肉中的(\textrm{Ca}^{2+})离子浓度，而是其由大脑控制的神经系统向肌肉发送电信号，导致(\textrm{Ca}^{2+})离子释放。这些离子扩散并在肌肉中增加浓度，导致激活。传输到肌肉以刺激收缩的电信号称为“兴奋”。在生物力学模型中，通常用符号(e(t))表示，也被视为在区间([0, 1])内的归一化量。

兴奋输入与激活状态之间的关系称为激活动力学。由于在生物力学模型中激活动力学非常普遍，SymPy 提供了activation模块，其中包含一些常见激活动力学模型的实现。这些模型包括基于[DeGroote2016]论文方程的零阶激活动力学和一阶激活动力学。接下来我们将手动实现这些模型，并展示它们与 SymPy 提供的类的关系。

零阶

激活动力学的最简单模型假设(\textrm{Ca}^{2+})离子的扩散是瞬时的。数学上，这给出了(a(t) = e(t))，一个零阶常微分方程。

>>> e = me.dynamicsymbols('e')
>>> e
e(t)
>>> a = e
>>> a
e(t) 

或者，您可以给(a(t))自己的dynamicsymbols并使用替换，在任何方程中用(e(t))替换它。

>>> a = me.dynamicsymbols('a')
>>> zeroth_order_activation = {a: e}
>>> a.subs(zeroth_order_activation)
e(t) 

SymPy 提供了ZerothOrderActivation类，在activation模块中。这个类必须用单个参数name实例化，将一个名称与实例关联起来。这个名称应该在每个实例中是唯一的。

>>> from sympy.physics.biomechanics import ZerothOrderActivation
>>> actz = ZerothOrderActivation('zeroth')
>>> actz
ZerothOrderActivation('zeroth') 

传递给name的参数试图确保在实例之间自动创建的dynamicsymbols对于(e(t))和(a(t))是唯一的。

>>> actz.excitation
e_zeroth(t)
>>> actz.activation
e_zeroth(t) 

ZerothOrderActivation 是 ActivationBase 的子类，为模型的所有具体激活动力学类提供一致的接口。这包括一个方法来检查与模型相关联的常微分方程(s)。由于零阶激活动力学对应于零阶常微分方程，因此返回一个空列矩阵。

>>> actz.rhs()
Matrix(0, 1, []) 

一阶

在实践中，(\textrm{Ca}^{2+})离子的扩散和浓度增加并非瞬时完成。在真实的生物肌肉中，兴奋度的步增将导致激活的平稳逐渐增加。[DeGroote2016] 使用一阶常微分方程模型来描述：

[\begin{split}\frac{da}{dt} &= \left( \frac{1}{\tau_a \left(1 + 3a(t)\right)} (1 + 2f) + \frac{1 + 3a(t)}{4\tau_d} (1 - 2f) \right) \left(e(t) - a(t) \right) \ f &= \frac{1}{2} \tanh{\left(b \left(e(t) -a(t)\right)\right)}\end{split}]

其中 (\tau_a) 是激活的时间常数，(\tau_d) 是去激活的时间常数，(b) 是平滑系数。

>>> tau_a, tau_d, b = sm.symbols('tau_a, tau_d, b')
>>> f = sm.tanh(b*(e - a))/2
>>> dadt = ((1/(tau_a*(1 + 3*a)))*(1 + 2*f) + ((1 + 3*a)/(4*tau_d))*(1 - 2*f))*(e - a) 

然后可以使用这个一阶常微分方程来在模拟中传播输入 (e(t)) 下的状态 (a(t))。

与之前一样，SymPy 在 activation 模块中提供了类 FirstOrderActivationDeGroote2016。这个类是 ActivationBase 的另一个子类，使用从 [DeGroote2016] 定义的一阶激活动力学模型。这个类必须用四个参数来实例化：名称和三个 sympifiable 对象，用来表示三个常数 (\tau_a)、(\tau_d) 和 (b)。

>>> from sympy.physics.biomechanics import FirstOrderActivationDeGroote2016
>>> actf = FirstOrderActivationDeGroote2016('first', tau_a, tau_d, b)
>>> actf.excitation
e_first(t)
>>> actf.activation
a_first(t) 

这一阶常微分方程可以像之前一样访问，但这次返回一个长度为 1 的列向量。

>>> actf.rhs()
Matrix([[((1/2 - tanh(b*(-a_first(t) + e_first(t)))/2)*(3*a_first(t)/2 + 1/2)/tau_d + (tanh(b*(-a_first(t) + e_first(t)))/2 + 1/2)/(tau_a*(3*a_first(t)/2 + 1/2)))*(-a_first(t) + e_first(t))]]) 

您还可以使用建议的每个常数的值来实例化该类。这些值为：(\tau_a = 0.015)，(\tau_d = 0.060)，和 (b = 10.0)。

>>> actf2 = FirstOrderActivationDeGroote2016.with_defaults('first')
>>> actf2.rhs()
Matrix([[((1/2 - tanh(10.0*a_first(t) - 10.0*e_first(t))/2)/(0.0225*a_first(t) + 0.0075) + 16.6666666666667*(3*a_first(t)/2 + 1/2)*(tanh(10.0*a_first(t) - 10.0*e_first(t))/2 + 1/2))*(-a_first(t) + e_first(t))]])
>>> constants = {tau_a: sm.Float('0.015'), tau_d: sm.Float('0.060'), b: sm.Float('10.0')}
>>> actf.rhs().subs(constants)
Matrix([[(66.6666666666667*(1/2 - tanh(10.0*a_first(t) - 10.0*e_first(t))/2)/(3*a_first(t)/2 + 1/2) + 16.6666666666667*(3*a_first(t)/2 + 1/2)*(tanh(10.0*a_first(t) - 10.0*e_first(t))/2 + 1/2))*(-a_first(t) + e_first(t))]]) 

自定义

要创建您自己的激活动力学模型，您可以子类化ActivationBase并重写抽象方法。具体的类将符合预期的 API，并且会自动集成到sympy.physics.mechanics和sympy.physics.biomechanics的其余部分中。

肌肉肌腱曲线

多年来，已发布了许多不同配置的希尔型肌肉模型，其中包含不同的串联和并联元素组合。我们将考虑模型的一个非常常见版本，将肌腱建模为串联元件，而肌肉纤维则建模为三个并联元件：弹性元件、收缩元件和阻尼器。

示意图显示了四元素希尔型肌肉模型。(SE)是代表肌腱的串联元件，(CE)是收缩元件，(EE)是代表肌肉纤维弹性的并联元件，(DE)是阻尼器。

每个组成部分通常具有描述其特性曲线。接下来的子章节将描述并实现文章中描述的特性曲线[DeGroote2016]。

肌腱力长度

通常将肌腱建模为刚性（不可伸长）和弹性元件。如果肌腱被视为刚性，肌腱长度不会改变，肌肉纤维的长度直接随着肌肉肌腱长度的变化而变化。刚性肌腱不具有相关的特性曲线；它本身没有任何产生力量的能力，只是直接传递由肌肉纤维产生的力量。

如果肌腱是弹性的，则通常被建模为非线性弹簧。因此，我们有了第一个特性曲线，即肌腱力长度曲线，它是标准化肌腱长度的函数：

[\tilde{l}^T = \frac{l^T}{l^T_{slack}}]

其中 (l^T) 是肌腱长度，(l^T_{slack}) 是“松弛肌腱长度”，一个代表在无力情况下的肌腱长度的常数。肌肉肌腱特性曲线是以“标准化”（或“无量纲”）量表述的，例如 (\tilde{l}^T)，因为这些曲线通用地适用于所有肌肉纤维和肌腱。通过选择不同的常数值，可以调整其属性以模拟特定的肌肉肌腱。在肌腱力长度特性中，通过调节 (l^T_{slack}) 来实现。较小的此常数值会导致肌腱更加僵硬。

肌腱力长度曲线(fl^T\left(\tilde{l}^T\right))的方程来自[DeGroote2016]：

[fl^T\left(\tilde{l}^T\right) = c_0 \exp{c_3 \left( \tilde{l}^T - c_1 \right)} - c_2]

要在 SymPy 中实现这一点，我们需要一个表示时间变化的动态符号，代表(\tilde{l}^T)，以及四个表示四个常数的符号。

>>> l_T_tilde = me.dynamicsymbols('l_T_tilde')
>>> c0, c1, c2, c3 = sm.symbols('c0, c1, c2, c3') 

>>> fl_T = c0*sm.exp(c3*(l_T_tilde - c1)) - c2
>>> fl_T
c0*exp(c3*(-c1 + l_T_tilde(t))) - c2 

或者，我们可以根据(l^T)和(l^T_{slack})来定义这个。

>>> l_T = me.dynamicsymbols('l_T')
>>> l_T_slack = sm.symbols('l_T_slack') 

>>> fl_T = c0*sm.exp(c3*(l_T/l_T_slack - c1)) - c2
>>> fl_T
c0*exp(c3*(-c1 + l_T(t)/l_T_slack)) - c2 

SymPy 的biomechanics模块提供了这条确切曲线的类，TendonForceLengthDeGroote2016。它可以用五个参数实例化。第一个参数是(\tilde{l}^T)，它不一定是一个符号；它可以是一个表达式。其余四个参数都是常数。预期这些将是常数或可以转换为数值的 SymPy 类型。

>>> from sympy.physics.biomechanics import TendonForceLengthDeGroote2016 

>>> fl_T2 = TendonForceLengthDeGroote2016(l_T/l_T_slack, c0, c1, c2, c3)
>>> fl_T2
TendonForceLengthDeGroote2016(l_T(t)/l_T_slack, c0, c1, c2, c3) 

这个类是Function的子类，因此实现了用于替换、评估、微分等的通常 SymPy 方法。doit方法允许访问曲线的方程。

>>> fl_T2.doit()
c0*exp(c3*(-c1 + l_T(t)/l_T_slack)) - c2 

该类提供了一个替代构造函数，允许它使用推荐的常数值在[DeGroote2016]中预填充。这需要一个参数，再次对应于(\tilde{l}^T)，它可以是符号或表达式。

>>> fl_T3 = TendonForceLengthDeGroote2016.with_defaults(l_T/l_T_slack)
>>> fl_T3
TendonForceLengthDeGroote2016(l_T(t)/l_T_slack, 0.2, 0.995, 0.25, 33.93669377311689) 

在上述内容中，常数已被替换为 SymPy 数值类型的实例，如Float。

TendonForceLengthDeGroote2016类还支持代码生成，因此与 SymPy 的代码打印机无缝集成。为了可视化这条曲线，我们可以使用lambdify在函数的实例上，它将创建一个可调用函数，用于对给定的(\tilde{l}^T)值进行评估。合理的(\tilde{l}^T)值落在范围([0.95, 1.05])内，我们将在下面进行绘制。

>>> import matplotlib.pyplot as plt
>>> import numpy as np
>>> from sympy.physics.biomechanics import TendonForceLengthDeGroote2016 

>>> l_T_tilde = me.dynamicsymbols('l_T_tilde')
>>> fl_T = TendonForceLengthDeGroote2016.with_defaults(l_T_tilde)
>>> fl_T_callable = sm.lambdify(l_T_tilde, fl_T)
>>> l_T_tilde_num = np.linspace(0.95, 1.05) 

>>> fig, ax = plt.subplots()
>>> _ = ax.plot(l_T_tilde_num, fl_T_callable(l_T_tilde_num))
>>> _ = ax.set_xlabel('Normalized tendon length')
>>> _ = ax.set_ylabel('Normalized tendon force-length') 

(png, hires.png, pdf)

推导弹性肌腱模型的肌腱肌腱动力学方程时，了解肌腱力长度特征曲线的反函数是很有用的。[DeGroote2016] 中定义的曲线是可解析可逆的，这意味着我们可以直接根据给定的 (fl^T\left(\tilde{l}^T\right)) 确定 (\tilde{l}^T = \left[fl^T\left(\tilde{l}^T\right)\right]^{-1})。

[\tilde{l}^T = \left[fl^T\left(\tilde{l}^T\right)\right]^{-1} = \frac{\log{\frac{fl^T + c_2}{c_0}}}{c_3} + c_1]

在 生物力学 中也有一个类 TendonForceLengthInverseDeGroote2016，与 TendonForceLengthDeGroote2016 行为相同。可以用五个参数实例化它，第一个参数是 (fl^T)，后面是四个常数，或者可以使用单参数的构造函数来实例化 (fl^T)。

>>> from sympy.physics.biomechanics import TendonForceLengthInverseDeGroote2016 

>>> fl_T_sym =me.dynamicsymbols('fl_T')
>>> fl_T_inv = TendonForceLengthInverseDeGroote2016(fl_T_sym, c0, c1, c2, c3)
>>> fl_T_inv
TendonForceLengthInverseDeGroote2016(fl_T(t), c0, c1, c2, c3) 

>>> fl_T_inv2 = TendonForceLengthInverseDeGroote2016.with_defaults(fl_T_sym)
>>> fl_T_inv2
TendonForceLengthInverseDeGroote2016(fl_T(t), 0.2, 0.995, 0.25, 33.93669377311689) 

纤维被动力长度

用于建模肌肉纤维的第一个元素是纤维被动力长度。这本质上是另一个非线性弹簧，代表肌肉纤维的弹性特性。描述这个元素的特征曲线是关于归一化肌肉纤维长度的函数：

[\tilde{l}^M = \frac{l^M}{l^M_{opt}}]

其中 (l^M) 是肌肉纤维长度，(l^M_{opt}) 是“最佳纤维长度”，一个表示肌肉纤维在不产生被动弹性力的长度（也是肌肉纤维能够产生最大主动力的长度）。像调整 (l^T_{slack}) 以通过肌腱力长度特征来改变建模肌腱的刚度属性一样，我们可以调整 (l^M_{opt}) 以改变肌肉纤维的被动特性；减小 (l^M_{opt}) 将使建模肌肉纤维变得更加坚硬。

来自 [DeGroote2016] 的肌肉纤维被动力长曲线 (fl^M_{pas}\left(\tilde{l}^M\right)) 的方程如下图所示：

[fl^M_{pas} = \frac{\frac{\exp{c_1 \left(\tilde{l^M} - 1\right)}}{c_0} - 1}{\exp{c_1} - 1}]

类似之前，在 SymPy 中实现这个过程需要一个表示时间变化的动态符号 (\tilde{l}^M) 和表示两个常数的符号。

>>> l_M_tilde = me.dynamicsymbols('l_M_tilde')
>>> c0, c1 = sm.symbols('c0, c1') 

>>> fl_M_pas = (sm.exp(c1*(l_M_tilde - 1)/c0) - 1)/(sm.exp(c1) - 1)
>>> fl_M_pas
(exp(c1*(l_M_tilde(t) - 1)/c0) - 1)/(exp(c1) - 1) 

或者，我们可以用(l^M)和(l^M_{opt})定义这个。

>>> l_M = me.dynamicsymbols('l_M')
>>> l_M_opt = sm.symbols('l_M_opt') 

>>> fl_M_pas2 = (sm.exp(c1*(l_M/l_M_opt - 1)/c0) - 1)/(sm.exp(c1) - 1)
>>> fl_M_pas2
(exp(c1*(-1 + l_M(t)/l_M_opt)/c0) - 1)/(exp(c1) - 1) 

另外，SymPy 中的biomechanics模块提供了一个专门处理此曲线的类，FiberForceLengthPassiveDeGroote2016。可以用三个参数实例化此类。第一个参数是(\tilde{l}^M)，不一定是一个符号，可以是一个表达式。另外两个参数都是常数。预期这些将是常数或可以符号化的数值。

>>> from sympy.physics.biomechanics import FiberForceLengthPassiveDeGroote2016 

>>> fl_M_pas2 = FiberForceLengthPassiveDeGroote2016(l_M/l_M_opt, c0, c1)
>>> fl_M_pas2
FiberForceLengthPassiveDeGroote2016(l_M(t)/l_M_opt, c0, c1)
>>> fl_M_pas2.doit()
(exp(c1*(-1 + l_M(t)/l_M_opt)/c0) - 1)/(exp(c1) - 1) 

使用单参数的替代构造函数，可以创建一个已经预设使用[DeGroote2016]推荐常数值的实例。

>>> fl_M_pas3 = FiberForceLengthPassiveDeGroote2016.with_defaults(l_M/l_M_opt)
>>> fl_M_pas3
FiberForceLengthPassiveDeGroote2016(l_M(t)/l_M_opt, 0.6, 4.0)
>>> fl_M_pas3.doit()
2.37439874427164e-5*exp(6.66666666666667*l_M(t)/l_M_opt) - 0.0186573603637741 

(\tilde{l}^M)的合理值位于范围([0.0, 2.0])，我们将在下面绘制。

>>> import matplotlib.pyplot as plt
>>> import numpy as np
>>> from sympy.physics.biomechanics import FiberForceLengthPassiveDeGroote2016 

>>> l_M_tilde = me.dynamicsymbols('l_M_tilde')
>>> fl_M_pas = FiberForceLengthPassiveDeGroote2016.with_defaults(l_M_tilde)
>>> fl_M_pas_callable = sm.lambdify(l_M_tilde, fl_M_pas)
>>> l_M_tilde_num = np.linspace(0.0, 2.0) 

>>> fig, ax = plt.subplots()
>>> _ = ax.plot(l_M_tilde_num, fl_M_pas_callable(l_M_tilde_num))
>>> _ = ax.set_xlabel('Normalized fiber length')
>>> _ = ax.set_ylabel('Normalized fiber passive force-length') 

(png, hires.png, pdf)

要在制定肌腱动力学时，有时需要反转纤维被动力长度特性曲线。这条曲线的方程来自[DeGroote2016]，可以再次进行解析反转。

[\tilde{l}^M = \left[fl^M_{pas}\right]^{-1} = \frac{c_0 \log{\left(\exp{c_1} - 1\right)fl^M_{pas} + 1}}{c_1} + 1]

在biomechanics中也有一个用于此目的的类，FiberForceLengthPassiveInverseDeGroote2016。可以用三个参数实例化此类，第一个是(fl^M)，后跟一对常数，或者使用单参数的替代构造函数来设置(\tilde{l}^M)。

>>> from sympy.physics.biomechanics import FiberForceLengthPassiveInverseDeGroote2016 

>>> fl_M_pas_sym =me.dynamicsymbols('fl_M_pas')
>>> fl_M_pas_inv = FiberForceLengthPassiveInverseDeGroote2016(fl_M_pas_sym, c0, c1)
>>> fl_M_pas_inv
FiberForceLengthPassiveInverseDeGroote2016(fl_M_pas(t), c0, c1) 

>>> fl_M_pas_inv2 = FiberForceLengthPassiveInverseDeGroote2016.with_defaults(fl_M_pas_sym)
>>> fl_M_pas_inv2
FiberForceLengthPassiveInverseDeGroote2016(fl_M_pas(t), 0.6, 4.0) 

纤维主动力长度

当肌肉被激活时，它收缩产生力。这种现象由肌腱模型中平行纤维组分的收缩元素建模。纤维可以产生的力量取决于纤维的即时长度。描述纤维活动力长度曲线的特征曲线再次由 (\tilde{l}^M) 参数化。这条曲线呈“钟形”。对于非常小和非常大的 (\tilde{l}^M) 值，活动纤维力长度趋于零。当 (\tilde{l}^M = l^M_{opt}) 时，活动纤维力长度达到峰值，其值为 (0.0)。

来自 [DeGroote2016] 的纤维活动力长度曲线 (fl^M_{act}\left(\tilde{l}^M\right)) 的方程为：

[fl^M_{act}\left(\tilde{l}^M\right) = c_0 \exp{-\frac{1}{2}\left(\frac{\tilde{l}^M - c_1}{\left(c_2 + c_3 \tilde{l}^M\right)}\right)²} + c_4 \exp{-\frac{1}{2}\left(\frac{\tilde{l}^M - c_5}{\left(c_6 + c_7 \tilde{l}^M\right)}\right)²} + c_8 \exp{-\frac{1}{2}\left(\frac{\tilde{l}^M - c_9}{\left(c_{10} + c_{11} \tilde{l}^M\right)}\right)²}]

在 SymPy 中实现这一点，我们需要一个表示时间变化的动态符号，表示 (\tilde{l}^M)，以及十二个表示十二个常数的符号。

>>> constants = sm.symbols('c0:12')
>>> c0, c1, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11 = constants 

>>> fl_M_act = (c0*sm.exp(-(((l_M_tilde - c1)/(c2 + c3*l_M_tilde))**2)/2) + c4*sm.exp(-(((l_M_tilde - c5)/(c6 + c7*l_M_tilde))**2)/2) + c8*sm.exp(-(((l_M_tilde - c9)/(c10 + c11*l_M_tilde))**2)/2))
>>> fl_M_act
c0*exp(-(-c1 + l_M_tilde(t))**2/(2*(c2 + c3*l_M_tilde(t))**2)) + c4*exp(-(-c5 + l_M_tilde(t))**2/(2*(c6 + c7*l_M_tilde(t))**2)) + c8*exp(-(-c9 + l_M_tilde(t))**2/(2*(c10 + c11*l_M_tilde(t))**2)) 

SymPy 提供的确切曲线类是 FiberForceLengthActiveDeGroote2016，它可以用十三个参数实例化。第一个参数是 (\tilde{l}^M)，不一定需要是一个符号，可以是一个表达式。接下来的十二个参数都是常数。预期这些将是常数或可 sympify 的数值。

>>> from sympy.physics.biomechanics import FiberForceLengthActiveDeGroote2016 

>>> fl_M_act2 = FiberForceLengthActiveDeGroote2016(l_M/l_M_opt, *constants)
>>> fl_M_act2
FiberForceLengthActiveDeGroote2016(l_M(t)/l_M_opt, c0, c1, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11)
>>> fl_M_act2.doit()
c0*exp(-(-c1 + l_M(t)/l_M_opt)**2/(2*(c2 + c3*l_M(t)/l_M_opt)**2)) + c4*exp(-(-c5 + l_M(t)/l_M_opt)**2/(2*(c6 + c7*l_M(t)/l_M_opt)**2)) + c8*exp(-(-c9 + l_M(t)/l_M_opt)**2/(2*(c10 + c11*l_M(t)/l_M_opt)**2)) 

使用单参数 (\tilde{l}^M) 的替代构造函数，我们可以创建一个实例，其常数值推荐在 [DeGroote2016] 中。

>>> fl_M_act3 = FiberForceLengthActiveDeGroote2016.with_defaults(l_M/l_M_opt)
>>> fl_M_act3
FiberForceLengthActiveDeGroote2016(l_M(t)/l_M_opt, 0.814, 1.06, 0.162, 0.0633, 0.433, 0.717, -0.0299, 0.2, 0.1, 1.0, 0.354, 0.0)
>>> fl_M_act3.doit()
0.1*exp(-3.98991349867535*(-1 + l_M(t)/l_M_opt)**2) + 0.433*exp(-12.5*(-0.717 + l_M(t)/l_M_opt)**2/(-0.1495 + l_M(t)/l_M_opt)**2) + 0.814*exp(-21.4067977442463*(-1 + 0.943396226415094*l_M(t)/l_M_opt)**2/(1 + 0.390740740740741*l_M(t)/l_M_opt)**2) 

合理的 (\tilde{l}^M) 值落在范围 ([0.0, 2.0]) 内，我们将在下面绘制。

>>> import matplotlib.pyplot as plt
>>> import numpy as np
>>> from sympy.physics.biomechanics import FiberForceLengthActiveDeGroote2016 

>>> l_M_tilde = me.dynamicsymbols('l_M_tilde')
>>> fl_M_act = FiberForceLengthActiveDeGroote2016.with_defaults(l_M_tilde)
>>> fl_M_act_callable = sm.lambdify(l_M_tilde, fl_M_act)
>>> l_M_tilde_num = np.linspace(0.0, 2.0) 

>>> fig, ax = plt.subplots()
>>> _ = ax.plot(l_M_tilde_num, fl_M_act_callable(l_M_tilde_num))
>>> _ = ax.set_xlabel('Normalized fiber length')
>>> _ = ax.set_ylabel('Normalized fiber active force-length') 

(png, hires.png, pdf)

对于纤维活动力长度特征曲线，不存在反曲线，因为每个 (fl^M_{act}) 的值对应多个 (\tilde{l}^M) 的值。

纤维力-速度

收缩元素产生的力量也是其长度变化速度的函数。描述收缩元素动态中与速度有关部分的特征曲线是肌肉纤维长度变化速度的函数：

[\tilde{v}^M = \frac{v^M}{v^M_{max}}]

这里的 (v^M) 是肌肉纤维的拉伸速度，(v^M_{max}) 是“最大纤维速度”，代表肌肉纤维在同心收缩时不能产生任何收缩力的速度常数。(v^M_{max}) 通常被赋予一个值 (10 l^M_{opt})。

来自[DeGroote2016]的纤维力-速度曲线 (fv^M\left(\tilde{v}^M\right)) 的方程是：

[fv^M\left(\tilde{v}^M\right) = c_0 \log{\left(c_1 \tilde{v}^M + c_2\right) + \sqrt{\left(c_1 \tilde{v}^M + c_2\right)² + 1}} + c_3]

类似之前，要在 SymPy 中实现这一点，我们需要一个表示 (\tilde{v}^M) 的时间变化动态符号，以及四个表示四个常数的符号。

>>> v_M_tilde = me.dynamicsymbols('v_M_tilde')
>>> c0, c1, c2, c3 = sm.symbols('c0, c1, c2, c3') 

>>> fv_M = c0*sm.log(c1*v_M_tilde + c2 + sm.sqrt((c1*v_M_tilde + c2)**2 + 1)) + c3
>>> fv_M
c0*log(c1*v_M_tilde(t) + c2 + sqrt((c1*v_M_tilde(t) + c2)**2 + 1)) + c3 

或者，我们可以根据 (v^M) 和 (v^M_{max}) 来定义这个曲线。

>>> v_M = me.dynamicsymbols('v_M')
>>> v_M_max = sm.symbols('v_M_max') 

>>> fv_M_pas2 = c0*sm.log(c1*v_M/v_M_max + c2 + sm.sqrt((c1*v_M/v_M_max + c2)**2 + 1)) + c3
>>> fv_M_pas2
c0*log(c1*v_M(t)/v_M_max + c2 + sqrt((c1*v_M(t)/v_M_max + c2)**2 + 1)) + c3 

SymPy 提供的精确曲线类是FiberForceVelocityDeGroote2016。它可以用五个参数实例化。第一个参数是 (\tilde{v}^M)，这不一定需要是一个符号，可以是一个表达式。接下来的四个参数都是常数。预期这些将是常数，或者是可以转换为数值的符号。

>>> from sympy.physics.biomechanics import FiberForceVelocityDeGroote2016 

>>> fv_M2 = FiberForceVelocityDeGroote2016(v_M/v_M_max, c0, c1, c2, c3)
>>> fv_M2
FiberForceVelocityDeGroote2016(v_M(t)/v_M_max, c0, c1, c2, c3)
>>> fv_M2.doit()
c0*log(c1*v_M(t)/v_M_max + c2 + sqrt((c1*v_M(t)/v_M_max + c2)**2 + 1)) + c3 

使用另一种构造函数，它接受一个参数 (\tilde{v}^M)，我们可以创建一个预先填充了推荐常数值的实例，这些值来自于[DeGroote2016]。

>>> fv_M3 = FiberForceVelocityDeGroote2016.with_defaults(v_M/v_M_max)
>>> fv_M3
FiberForceVelocityDeGroote2016(v_M(t)/v_M_max, -0.318, -8.149, -0.374, 0.886)
>>> fv_M3.doit()
0.886 - 0.318*log(8.149*sqrt((-0.0458952018652595 - v_M(t)/v_M_max)**2 + 0.0150588346410601) - 0.374 - 8.149*v_M(t)/v_M_max) 

合理的 (\tilde{v}^M) 值落在范围 ([-1.0, 1.0]) 内，我们将在下面进行绘制。

>>> import matplotlib.pyplot as plt
>>> import numpy as np
>>> from sympy.physics.biomechanics import FiberForceVelocityDeGroote2016 

>>> v_M_tilde = me.dynamicsymbols('v_M_tilde')
>>> fv_M = FiberForceVelocityDeGroote2016.with_defaults(v_M_tilde)
>>> fv_M_callable = sm.lambdify(v_M_tilde, fv_M)
>>> v_M_tilde_num = np.linspace(-1.0, 1.0) 

>>> fig, ax = plt.subplots()
>>> _ = ax.plot(l_M_tilde_num, fv_M_callable(v_M_tilde_num))
>>> _ = ax.set_xlabel('Normalized fiber velocity')
>>> _ = ax.set_ylabel('Normalized fiber force-velocity') 

(png, hires.png, pdf)

当制定肌腱动力学时，有时需要纤维力-速度特性曲线的反函数。这条曲线的方程来自于[DeGroote2016]，同样可以进行解析反演。

[\tilde{v}^M = \left[fv^M\right]^{-1} = \frac{\sinh{\frac{fv^M - c_3}{c_0}} - c_2}{c_1}]

这个公式也在biomechanics中有一个类，名为FiberForceVelocityInverseDeGroote2016。它可以用五个参数实例化，第一个是 (fv^M)，后面四个是常数；或者可以用另一种构造函数，只需要一个参数 (\tilde{v}^M)。

>>> from sympy.physics.biomechanics import FiberForceVelocityInverseDeGroote2016 

>>> fv_M_sym = me.dynamicsymbols('fv_M')
>>> fv_M_inv = FiberForceVelocityInverseDeGroote2016(fv_M_sym, c0, c1, c2, c3)
>>> fv_M_inv
FiberForceVelocityInverseDeGroote2016(fv_M(t), c0, c1, c2, c3) 

>>> fv_M_inv2 = FiberForceVelocityInverseDeGroote2016.with_defaults(fv_M_sym)
>>> fv_M_inv2
FiberForceVelocityInverseDeGroote2016(fv_M(t), -0.318, -8.149, -0.374, 0.886) 

纤维阻尼

在肌腱模型中，纤维阻尼可能是最简单的元素。它没有关联的特征曲线，通常被建模为简单的线性阻尼器。我们将使用 (\beta) 作为阻尼系数，阻尼力可以描述为：

[f_{damp} = \beta \tilde{v}^M]

[DeGroote2016] 建议使用 (\beta = 0.1)。然而，SymPy 默认使用 (\beta = 10)。在进行前向模拟或解决最优控制问题时，这种增加阻尼通常不会显著影响肌腱动力学，但已经经验性地发现可以显著改善方程的数值条件。

肌腱动力学

刚性肌腱动力学

由于不可伸展的肌腱允许直接用肌腱长度表示肌肉纤维的归一化长度，刚性肌腱肌腱动力学实现起来相当直接。当不可伸展的肌腱 (l^T = l^T_{slack})，归一化肌腱长度即为单位长度，(\tilde{l}^T = 1)。利用三角法，肌肉纤维长度可以表达为

[l^M = \sqrt{\left(l^{MT} - l^T\right)² + \left(l^M_{opt} \sin{\alpha_{opt}} \right)²}]

其中 (\alpha_{opt}) 是“最佳斜角”，是肌腱的另一个恒定特性，描述了肌肉纤维相对于与肌腱平行方向的角度（肌肉纤维的角度）。常见的简化假设是假设 (\alpha_{opt} = 0)，这样上述公式简化为

[l^M = \sqrt{\left(l^{MT} - l^T\right)² + \left(l^M_{opt}\right)²}]

当 (\tilde{l}^M = \frac{l^M}{l^M_{opt}}) 时，肌肉纤维速度可以表达为

[v^M = v^{MT} \frac{l^{MT} - l^T_{slack}}{l^M}]

如前所述，肌肉纤维可以归一化，(\tilde{v}^M = \frac{v^M}{v^M_{max}})。利用上述曲线，我们可以将归一化肌肉纤维力（(\tilde{F}^M)) 表达为归一化肌腱长度（(\tilde{l}^T))、归一化纤维长度（(\tilde{l}^M))、归一化纤维速度（(\tilde{v}^M)) 和激活（(a)) 的函数。

[\tilde{F}^M = a \cdot fl^M_{act}\left(\tilde{l}^M\right) \cdot fv^M\left(\tilde{v}^M\right) + fl^M_{pas}\left(\tilde{l}^M\right) + \beta \cdot \tilde{v}^M]

我们引入一个新的常数 (F^M_{max})，即“最大等长收缩力”，描述了在全激活和等长（(v^M = 0)）收缩状态下肌腱最多可以产生的力。考虑到斜角度，肌腱力（(F^T))，即施加在肌腱起点和插入点骨骼上的力，可以表示为：

[F^T = F^M_{max} \cdot F^M \cdot \sqrt{1 - \sin{\alpha_{opt}}²}]

我们可以使用上面介绍过的 SymPy 和肌腱曲线类来描述所有这些。我们将需要时间变化的动力学符号(l^{MT})、(v_{MT})和(a)。我们还需要常量符号(l^T_{slack})、(l^M_{opt})、(F^M_{max})、(v^M_{max})、(\alpha_{opt})和(\beta)。

>>> l_MT, v_MT, a = me.dynamicsymbols('l_MT, v_MT, a')
>>> l_T_slack, l_M_opt, F_M_max = sm.symbols('l_T_slack, l_M_opt, F_M_max')
>>> v_M_max, alpha_opt, beta = sm.symbols('v_M_max, alpha_opt, beta') 

>>> l_M = sm.sqrt((l_MT - l_T_slack)**2 + (l_M_opt*sm.sin(alpha_opt))**2)
>>> l_M
sqrt(l_M_opt**2*sin(alpha_opt)**2 + (-l_T_slack + l_MT(t))**2) 

>>> v_M = v_MT*(l_MT - l_T_slack)/l_M
>>> v_M
(-l_T_slack + l_MT(t))*v_MT(t)/sqrt(l_M_opt**2*sin(alpha_opt)**2 + (-l_T_slack + l_MT(t))**2) 

>>> fl_M_pas = FiberForceLengthPassiveDeGroote2016.with_defaults(l_M/l_M_opt)
>>> fl_M_pas
FiberForceLengthPassiveDeGroote2016(sqrt(l_M_opt**2*sin(alpha_opt)**2 + (-l_T_slack + l_MT(t))**2)/l_M_opt, 0.6, 4.0) 

>>> fl_M_act = FiberForceLengthActiveDeGroote2016.with_defaults(l_M/l_M_opt)
>>> fl_M_act
FiberForceLengthActiveDeGroote2016(sqrt(l_M_opt**2*sin(alpha_opt)**2 + (-l_T_slack + l_MT(t))**2)/l_M_opt, 0.814, 1.06, 0.162, 0.0633, 0.433, 0.717, -0.0299, 0.2, 0.1, 1.0, 0.354, 0.0) 

>>> fv_M = FiberForceVelocityDeGroote2016.with_defaults(v_M/v_M_max)
>>> fv_M
FiberForceVelocityDeGroote2016((-l_T_slack + l_MT(t))*v_MT(t)/(v_M_max*sqrt(l_M_opt**2*sin(alpha_opt)**2 + (-l_T_slack + l_MT(t))**2)), -0.318, -8.149, -0.374, 0.886) 

>>> F_M = a*fl_M_act*fv_M + fl_M_pas + beta*v_M/v_M_max
>>> F_M
beta*(-l_T_slack + l_MT(t))*v_MT(t)/(v_M_max*sqrt(l_M_opt**2*sin(alpha_opt)**2 + (-l_T_slack + l_MT(t))**2)) + a(t)*FiberForceLengthActiveDeGroote2016(sqrt(l_M_opt**2*sin(alpha_opt)**2 + (-l_T_slack + l_MT(t))**2)/l_M_opt, 0.814, 1.06, 0.162, 0.0633, 0.433, 0.717, -0.0299, 0.2, 0.1, 1.0, 0.354, 0.0)*FiberForceVelocityDeGroote2016((-l_T_slack + l_MT(t))*v_MT(t)/(v_M_max*sqrt(l_M_opt**2*sin(alpha_opt)**2 + (-l_T_slack + l_MT(t))**2)), -0.318, -8.149, -0.374, 0.886) + FiberForceLengthPassiveDeGroote2016(sqrt(l_M_opt**2*sin(alpha_opt)**2 + (-l_T_slack + l_MT(t))**2)/l_M_opt, 0.6, 4.0) 

>>> F_T = F_M_max*F_M*sm.sqrt(1 - sm.sin(alpha_opt)**2)
>>> F_T
F_M_max*sqrt(1 - sin(alpha_opt)**2)*(beta*(-l_T_slack + l_MT(t))*v_MT(t)/(v_M_max*sqrt(l_M_opt**2*sin(alpha_opt)**2 + (-l_T_slack + l_MT(t))**2)) + a(t)*FiberForceLengthActiveDeGroote2016(sqrt(l_M_opt**2*sin(alpha_opt)**2 + (-l_T_slack + l_MT(t))**2)/l_M_opt, 0.814, 1.06, 0.162, 0.0633, 0.433, 0.717, -0.0299, 0.2, 0.1, 1.0, 0.354, 0.0)*FiberForceVelocityDeGroote2016((-l_T_slack + l_MT(t))*v_MT(t)/(v_M_max*sqrt(l_M_opt**2*sin(alpha_opt)**2 + (-l_T_slack + l_MT(t))**2)), -0.318, -8.149, -0.374, 0.886) + FiberForceLengthPassiveDeGroote2016(sqrt(l_M_opt**2*sin(alpha_opt)**2 + (-l_T_slack + l_MT(t))**2)/l_M_opt, 0.6, 4.0)) 

SymPy 提供了MusculotendonDeGroote2016类的刚性肌腱动力学实现，我们将在下面构建一个完整的简单肌腱模型时展示。

弹性肌腱动力学

弹性肌腱动力学更为复杂，因为我们不能直接用肌腱长度表达肌肉肌腱长度。相反，我们必须将肌腱中经历的力与肌肉纤维产生的力联系起来，确保两者处于平衡状态。为了在肌腱肌腱动力学中引入额外的状态变量，我们需要引入一个额外的一阶常微分方程。我们可以对这种状态进行多种选择，但其中一个最直观的选择可能是使用(\tilde{l}^M)。因此，我们需要创建肌腱力((F^T))的表达式，以及对(\frac{d \tilde{l}^M}{dt})进行一阶常微分方程的处理。(l^M)、(l^T)和(\tilde{l}^T)的计算方式与刚性肌腱动力学类似，记住我们已经有(\tilde{l}^M)作为已知值，因为它是一个状态变量。

[\begin{split}l^M &= \tilde{l}^M \cdot l^M_{opt} \ l^T &= l^{MT} - \sqrt{\left(l^M\right)² - \left(l^M_{opt} \sin{\alpha_{opt}}\right)²} \ \tilde{l}^T &= \frac{l^T}{l^T_{slack}}\end{split}]

使用(\tilde{l}^T)和肌腱力长度曲线（(fl^T\left(\tilde{l}^T\right))），我们可以编写归一化和绝对肌腱力的方程：

[\begin{split}\tilde{F}^T &= fl^T\left(\tilde{l}^T\right) \ F^T &= F^M_{max} \cdot \tilde{F}^T\end{split}]

要表达(F^M)，我们需要知道肌腱角的余弦值（(\cos{\alpha})）。我们可以使用三角学来编写此方程：

[\begin{split}\cos{\alpha} &= \frac{l^{MT} - l^T}{l^M} \ F^M &= \frac{F^T}{\cos{\alpha}}\end{split}]

如果我们假设阻尼系数(\beta = 0)，我们可以重新排列肌纤维力方程：

[\tilde{F}^M = a \cdot fl^M_{act}\left(\tilde{l}^M\right) \cdot fv^M\left(\tilde{v}^M\right) + fl^M_{pas}\left(\tilde{l}^M\right) + \beta \cdot \tilde{v}^M]

给出(fv^M\left(\tilde{v}^M\right))：

[fv^M\left(\tilde{v}^M\right) = \frac{\tilde{F}^M - fl^M_{pas}\left(\tilde{l}^M\right)}{a \cdot fl^M_{act}\left(\tilde{l}^M\right)}]

使用反向纤维力-速度曲线，(\left[fv^M\left(\tilde{v}^M\right)\right]^{-1})，并对(\tilde{l}^M)关于时间求导，我们最终可以写出(\frac{d \tilde{l}^M}{dt})的方程：

[\frac{d \tilde{l}^M}{dt} = \frac{v^M_{max}}{l^M_{opt}} \tilde{v}^M]

要制定这些弹性肌腱肌肉肌腱动力学，我们不得不假设(\beta = 0)，这在前向模拟和最优控制问题中是次优的。可以制定具有阻尼的弹性肌腱肌肉肌腱动力学，但这需要一个更复杂的公式，除了额外的输入变量外还需要额外的状态变量，因此肌腱肌肉肌腱动力学必须作为代数微分方程而不是普通微分方程来实施。这些类型的具体公式不会在此处讨论，但感兴趣的读者可以参考MusculotendonDeGroote2016的文档字符串，其中它们已实施。

一个简单的肌腱模型

要展示肌肉对简单系统的影响，我们可以模拟一个质量为(m)的粒子，在重力的影响下，由肌肉抵抗重力拉动质量。质量(m)具有单一的广义坐标(q)和广义速度(u)，以描述其位置和运动。以下代码建立了运动学和重力作用以及相关的粒子：

>>> import pprint
>>> import sympy as sm
>>> import sympy.physics.mechanics as me 

>>> q, u = me.dynamicsymbols('q, u', real=True)
>>> m, g = sm.symbols('m, g', real=True, positive=True) 

>>> N = me.ReferenceFrame('N')
>>> O, P = sm.symbols('O, P', cls=me.Point) 

>>> P.set_pos(O, q*N.x)
>>> O.set_vel(N, 0)
>>> P.set_vel(N, u*N.x) 

>>> gravity = me.Force(P, m*g*N.x) 

>>> block = me.Particle('block', P, m) 

SymPy 生物力学包括肌腱肌肉肌腱执行器模型。在这里，我们将使用特定的肌腱肌肉肌腱模型实现。肌腱肌肉肌腱执行器使用两个输入组件实例化，即路径和激活动力学模型。执行器必须沿连接肌肉起始点和插入点的路径执行。我们的起始点将附着在固定点(O)上，插入到移动粒子(P)上。

>>> from sympy.physics.mechanics.pathway import LinearPathway 

>>> muscle_pathway = LinearPathway(O, P) 

一条路径具有附着点：

>>> muscle_pathway.attachments
(O, P) 

并知道末端附着点之间的长度以及两个附着点之间的相对速度：

>>> muscle_pathway.length
Abs(q(t))
>>> muscle_pathway.extension_velocity
sign(q(t))*Derivative(q(t), t) 

最后，路径可以确定作用在两个附着点上的力的大小：

>>> muscle_pathway.to_loads(m*g)
[(O, - g*m*q(t)/Abs(q(t))*N.x), (P, g*m*q(t)/Abs(q(t))*N.x)] 

激活动力学模型表示了一组代数或普通微分方程，这些方程将肌肉兴奋与肌肉激活联系起来。在我们的情况下，我们将使用一个一阶普通微分方程，该方程从兴奋(e(t))产生平滑但延迟的激活(a(t))。

>>> from sympy.physics.biomechanics import FirstOrderActivationDeGroote2016
>>> muscle_activation = FirstOrderActivationDeGroote2016.with_defaults('muscle') 

激活模型具有状态变量(\mathbf{x})，输入变量(\mathbf{r})，以及一些常数参数(\mathbf{p})：

>>> muscle_activation.x
Matrix([[a_muscle(t)]])
>>> muscle_activation.r
Matrix([[e_muscle(t)]])
>>> muscle_activation.p
Matrix(0, 1, []) 

请注意，常量参数的返回值为空。如果我们正常实例化 FirstOrderActivationDeGroote2016，那么我们将不得不提供 (\tau_{a})、(\tau_{d}) 和 (b) 的三个值。如果这些是 Symbol 对象，那么它们将显示在返回的 MutableDenseMatrix 中。

这些与其一阶微分方程 (\dot{a} = f(a, e, t)) 相关：

>>> muscle_activation.rhs()
Matrix([[((1/2 - tanh(10.0*a_muscle(t) - 10.0*e_muscle(t))/2)/(0.0225*a_muscle(t) + 0.0075) + 16.6666666666667*(3*a_muscle(t)/2 + 1/2)*(tanh(10.0*a_muscle(t) - 10.0*e_muscle(t))/2 + 1/2))*(-a_muscle(t) + e_muscle(t))]]) 

通过路径和激活动力学，使用它们创建的肌腱肌模型需要一些参数来定义肌肉和肌腱的特定属性。您需要指定肌腱松弛长度、峰值等长力、最佳纤维长度、最大纤维速度、最佳偏角和纤维阻尼系数。

>>> from sympy.physics.biomechanics import MusculotendonDeGroote2016 

>>> F_M_max, l_M_opt, l_T_slack = sm.symbols('F_M_max, l_M_opt, l_T_slack', real=True)
>>> v_M_max, alpha_opt, beta = sm.symbols('v_M_max, alpha_opt, beta', real=True) 

>>> muscle = MusculotendonDeGroote2016.with_defaults(
...     'muscle',
...     muscle_pathway,
...     muscle_activation,
...     tendon_slack_length=l_T_slack,
...     peak_isometric_force=F_M_max,
...     optimal_fiber_length=l_M_opt,
...     maximal_fiber_velocity=v_M_max,
...     optimal_pennation_angle=alpha_opt,
...     fiber_damping_coefficient=beta,
... )
... 

因为这种肌腱肌致动器具有刚性肌腱模型，它具有与激活模型相同的状态和常微分方程：

>>> muscle.musculotendon_dynamics
0
>>> muscle.x
Matrix([[a_muscle(t)]])
>>> muscle.r
Matrix([[e_muscle(t)]])
>>> muscle.p
Matrix([
[l_T_slack],
[  F_M_max],
[  l_M_opt],
[  v_M_max],
[alpha_opt],
[     beta]])
>>> muscle.rhs()
Matrix([[((1/2 - tanh(10.0*a_muscle(t) - 10.0*e_muscle(t))/2)/(0.0225*a_muscle(t) + 0.0075) + 16.6666666666667*(3*a_muscle(t)/2 + 1/2)*(tanh(10.0*a_muscle(t) - 10.0*e_muscle(t))/2 + 1/2))*(-a_muscle(t) + e_muscle(t))]]) 

肌腱肌提供额外的常微分方程以及应用于路径的肌肉特定力量：

>>> muscle_loads = muscle.to_loads()
>>> pprint.pprint(muscle_loads)
[Force(point=O, force=F_M_max*(beta*(-l_T_slack + Abs(q(t)))*sign(q(t))*Derivative(q(t), t)/(v_M_max*sqrt(l_M_opt**2*sin(alpha_opt)**2 + (-l_T_slack + Abs(q(t)))**2)) + a_muscle(t)*FiberForceLengthActiveDeGroote2016(sqrt(l_M_opt**2*sin(alpha_opt)**2 + (-l_T_slack + Abs(q(t)))**2)/l_M_opt, 0.814, 1.06, 0.162, 0.0633, 0.433, 0.717, -0.0299, 0.2, 0.1, 1.0, 0.354, 0.0)*FiberForceVelocityDeGroote2016((-l_T_slack + Abs(q(t)))*sign(q(t))*Derivative(q(t), t)/(v_M_max*sqrt(l_M_opt**2*sin(alpha_opt)**2 + (-l_T_slack + Abs(q(t)))**2)), -0.318, -8.149, -0.374, 0.886) + FiberForceLengthPassiveDeGroote2016(sqrt(l_M_opt**2*sin(alpha_opt)**2 + (-l_T_slack + Abs(q(t)))**2)/l_M_opt, 0.6, 4.0))*q(t)*cos(alpha_opt)/Abs(q(t))*N.x),
 Force(point=P, force=- F_M_max*(beta*(-l_T_slack + Abs(q(t)))*sign(q(t))*Derivative(q(t), t)/(v_M_max*sqrt(l_M_opt**2*sin(alpha_opt)**2 + (-l_T_slack + Abs(q(t)))**2)) + a_muscle(t)*FiberForceLengthActiveDeGroote2016(sqrt(l_M_opt**2*sin(alpha_opt)**2 + (-l_T_slack + Abs(q(t)))**2)/l_M_opt, 0.814, 1.06, 0.162, 0.0633, 0.433, 0.717, -0.0299, 0.2, 0.1, 1.0, 0.354, 0.0)*FiberForceVelocityDeGroote2016((-l_T_slack + Abs(q(t)))*sign(q(t))*Derivative(q(t), t)/(v_M_max*sqrt(l_M_opt**2*sin(alpha_opt)**2 + (-l_T_slack + Abs(q(t)))**2)), -0.318, -8.149, -0.374, 0.886) + FiberForceLengthPassiveDeGroote2016(sqrt(l_M_opt**2*sin(alpha_opt)**2 + (-l_T_slack + Abs(q(t)))**2)/l_M_opt, 0.6, 4.0))*q(t)*cos(alpha_opt)/Abs(q(t))*N.x)] 

这些负载由描述长度和速度关系与肌肉纤维力量的各种函数组成。

现在我们有了肌肉和肌腱产生的力，系统的运动方程可以通过例如 Kane's 方法形成：

>>> kane = me.KanesMethod(N, (q,), (u,), kd_eqs=(u - q.diff(),))
>>> Fr, Frs = kane.kanes_equations((block,), (muscle_loads + [gravity])) 

运动方程由运动学微分方程、动力学微分方程（牛顿第二定律）和肌肉激活微分方程组成。每个的显式形式可以如下形成：

>>> dqdt = u
>>> dudt = kane.forcing[0]/m
>>> dadt = muscle.rhs()[0] 

现在我们可以创建一个数值函数，根据状态、输入和常量参数评估运动方程。首先通过符号列出每个：

>>> a = muscle.a
>>> e = muscle.e
>>> state = [q, u, a]
>>> inputs = [e]
>>> constants = [m, g, F_M_max, l_M_opt, l_T_slack, v_M_max, alpha_opt, beta] 

然后，评估显式常微分方程右手边的数值函数是：

>>> eval_eom = sm.lambdify((state, inputs, constants), (dqdt, dudt, dadt)) 

另外，数值评估肌肉力量也会很有趣，因此也创建一个函数：

>>> force = muscle.force.xreplace({q.diff(): u})
>>> eval_force = sm.lambdify((state, constants), force) 

为了测试这些函数，我们需要一些合适的数值。这种肌肉能够产生最大 10N 的力量来提起一个质量为 0.5kg 的物体：

>>> import numpy as np
>>> p_vals = np.array([
...     0.5,  # m [kg]
...     9.81,  # g [m/s/s]
...     10.0,  # F_M_max [N]
...     0.18,  # l_M_opt [m]
...     0.17,  # l_T_slack [m]
...     10.0,  # v_M_max [m/s]
...     0.0,  # alpha_opt
...     0.1,  # beta
... ])
... 

我们的肌腱是刚性的，所以肌肉的长度将是 (q-l_{T_\textrm{slack}})，我们希望初始肌肉长度接近其产力峰值，因此我们选择 (q_0=l_{M_\textrm{opt}} + l_{T_\textrm{slack}})。还让肌肉具有小的初始激活以产生非零力量：

>>> x_vals = np.array([
...     p_vals[3] + p_vals[4],  # q [m]
...     0.0,  # u [m/s]
...     0.1,  # a [unitless]
... ])
... 

将兴奋设为 1.0 并测试数值函数：

>>> r_vals = np.array([
...     1.0,  # e
... ])
...
>>> eval_eom(x_vals, r_vals, p_vals)
(0.0, 7.817106179880225, 92.30769105034035)
>>> eval_force(x_vals, p_vals)
-0.9964469100598874 

这两个函数工作，因此我们现在可以模拟这个系统，看看肌肉如何提起负载：

>>> def eval_rhs(t, x):
...
...     r = np.array([1.0])
...
...     return eval_eom(x, r, p_vals)
...
>>> from scipy.integrate import solve_ivp
>>> t0, tf = 0.0, 6.0
>>> times = np.linspace(t0, tf, num=601)
>>> sol = solve_ivp(eval_rhs,
...                 (t0, tf),
...                 x_vals, t_eval=times)
...
>>> import matplotlib.pyplot as plt
>>> fig, axes = plt.subplots(4, 1, sharex=True)
>>> _ = axes[0].plot(sol.t, sol.y[0] - p_vals[4], label='length of muscle')
>>> _ = axes[0].set_ylabel('Distance [m]')
>>> _ = axes[1].plot(sol.t, sol.y[1], label=state[1])
>>> _ = axes[1].set_ylabel('Speed [m/s]')
>>> _ = axes[2].plot(sol.t, sol.y[2], label=state[2])
>>> _ = axes[2].set_ylabel('Activation')
>>> _ = axes[3].plot(sol.t, eval_force(sol.y, p_vals).T, label='force')
>>> _ = axes[3].set_ylabel('Force [N]')
>>> _ = axes[3].set_xlabel('Time [s]')
>>> _ = axes[0].legend(), axes[1].legend(), axes[2].legend(), axes[3].legend() 

(png, hires.png, pdf)

肌肉抵抗重力对质量施加力，使其达到平衡状态为 5 N。

参考文献

[DeGroote2016] (1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16)

De Groote, F., Kinney, A. L., Rao, A. V., & Fregly, B. J., 解决肌肉冗余问题的直接共轭最优控制问题公式评估, 生物医学工程学报, 44(10), (2016) pp. 2922-2936