SymPy-1-13-中文文档-一-

SymPy 1.13 中文文档（一）

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

安装

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

SymPy CAS 可以安装在几乎任何安装了 Python 的计算机上。SymPy 需要首先安装 mpmath Python 库。建议的安装方法是通过 Anaconda，它包括 mpmath 和其他几个有用的库。或者，一些 Linux 发行版提供了 SymPy 的软件包。

SymPy 官方支持 Python 3.8、3.9、3.10 和 PyPy。

Anaconda

Anaconda 是由 Continuum Analytics 提供的免费 Python 发行版，包含 SymPy、Matplotlib、IPython、NumPy 等许多科学计算的有用包。推荐使用它，因为只有在安装了特定的库时，SymPy 才能启用许多优秀的功能。例如，如果没有安装 Matplotlib，则只能启用简单的基于文本的绘图。通过 IPython 笔记本或 qtconsole，你可以通过运行 init_printing() 获得更漂亮的 (\mathrm{\LaTeX}) 打印输出。

如果你已经安装了 Anaconda 并希望将 SymPy 更新到最新版本，请使用：

conda update sympy 

Git

如果你希望为 SymPy 做贡献或者想要获取最新更新，请从 git 安装 SymPy。要下载代码库，请在命令行中执行以下操作：

git clone https://github.com/sympy/sympy.git 

要更新到最新版本，请进入你的代码库并执行：

git pull origin master 

如果你想安装 SymPy，但仍想使用 git 版本，请从你的代码库运行：

python -m pip install -e . 

这将导致安装的版本始终指向 git 目录中的版本。

其他方法

你也可以使用 pip 或从源代码安装 SymPy。此外，大多数 Linux 发行版和 Python 发行版都可以通过其包管理器安装 SymPy 的某个版本。以下是几个这样的 Python 发行版列表：

• Anaconda

• Enthought Deployment Manager

• ActivePython

• Spack

运行 SymPy

安装完成后，最好验证你刚刚安装的 SymPy 是否正常工作。为此，请启动 Python 并导入 SymPy 库：

$ python
>>> from sympy import * 

从这里开始执行一些简单的 SymPy 语句，比如以下示例：

>>> x = Symbol('x')
>>> limit(sin(x)/x, x, 0)
1
>>> integrate(1/x, x)
log(x) 

要了解如何有效使用 SymPy 的入门指南，请参阅 Introductory Tutorial。

mpmath

SymPy 1.0 之前的版本包含了 mpmath，但现在它作为一个外部依赖项存在。如果你使用 Anaconda 安装了 SymPy，它将已经包含 mpmath。使用：

conda install mpmath 

确保已安装。

如果你不想使用 Anaconda，你可以使用 pip install mpmath。

如果你的代码中使用 sympy.mpmath 调用 mpmath，你需要修改为直接调用 mpmath。如果你依赖的代码不容易更改，可以通过以下方式解决：

import sys
import mpmath
sys.modules['sympy.mpmath'] = mpmath 

在导入sympy.mpmath的代码之前。建议尽可能修改使用sympy.mpmath的代码直接使用mpmath。

问题

如果您对安装或 SymPy 总体有任何问题，请随时访问我们在Gitter上的聊天。此外，我们的邮件列表是社区支持的绝佳来源。

如果您认为存在 bug 或者希望请求功能，请开启一个问题票。

教程

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

对于任何初学 SymPy 或其特定功能的人来说，教程是开始的最佳地方。

入门教程

如果你是 SymPy 的新手，从这里开始。

物理教程

对于 SymPy 中的物理功能，从这里开始。

入门教程

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

本教程旨在为以前未使用过该库的人提供 SymPy 的介绍。本教程将介绍 SymPy 的许多功能，但不会详尽。事实上，本教程中显示的每个功能几乎都有更多选项或功能。SymPy 文档的其余部分作为 API 文档，详细列出了每个函数的所有功能和选项。

这些是本教程的目标：

• 给出一份指南，适合从未使用过 SymPy 的人（但已使用 Python 并了解必要的数学）。

• 以叙述格式编写，既易于理解又有趣。应该读起来像一本书。

• 提供富有洞见的例子和练习，帮助读者学习并使其有趣。

• 以逻辑顺序介绍概念。

• 使用良好的实践和习惯用法，避免反模式。避免导致反模式的函数或方法。不显示仅对高级用户有用的功能。

• 保持一致。如果有多种方法可以做到，只显示最佳方法。

• 为了避免不必要的重复，假设已经阅读了教程的前几部分。

欢迎随时在此教程或 SymPy 中提供反馈意见。只需写信到我们的邮件列表。

内容

• 准备工作

  • 安装

  • 练习

• 介绍

  • 符号计算是什么？

  • 更有趣的例子

  • 符号计算的力量

  • 为什么选择 SymPy？

• 注意事项

  • 注意事项

  • 等号

  • 两个最后注意事项：^ 和 /

  • 进一步阅读

• SymPy 特性

  • 基本操作

  • 打印

  • 简化

  • 微积分

  • 求解器

  • 矩阵

  • 高级表达式操作

• 接下来是什么

准备工作

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

这篇教程假设读者已经掌握了 Python 编程语言的基础知识。如果没有，可以参考官方 Python 教程，这个教程非常优秀。

这篇教程假设读者具备相当的数学背景。大多数示例需要低于微积分水平的知识，有些示例需要微积分水平的知识。一些高级特性需要更多的数学知识。如果遇到使用你不熟悉的数学函数的部分，你可以跳过它，或者用你更熟悉的类似函数替换它。或者在维基百科上查找这个函数并学习新知识。一些重要的数学概念将会在需要时介绍。

安装

首先，您需要安装 SymPy。请参阅安装指南。

练习

这篇教程是基于 2013 年在德克萨斯州奥斯汀举行的 SciPy 会议上进行的一次教程。该教程的网站在这里。网站上有视频链接、材料和 IPython 笔记本练习。特别推荐大家完成这篇教程时参与 IPython 笔记本练习。

简介

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

什么是符号计算？

符号计算处理数学对象的符号计算。这意味着数学对象被精确地表示，而不是近似表示，并且带有未评估变量的数学表达式以符号形式保留。

让我们举个例子。假设我们想要使用内置的 Python 函数来计算平方根。我们可能会做这样的事情

>>> import math
>>> math.sqrt(9)
3.0 

9 是一个完全平方数，所以我们得到了精确答案，3。但假设我们计算的是不是完全平方数的平方根

>>> math.sqrt(8)
2.82842712475 

这里我们得到了一个近似结果。2.82842712475 并不是 8 的精确平方根（事实上，8 的实际平方根不能用有限小数表示，因为它是一个无理数）。如果我们只关心 8 的小数形式的平方根，我们就结束了。

但假设我们想进一步。回想一下 (\sqrt{8} = \sqrt{4\cdot 2} = 2\sqrt{2})。我们很难从上面的结果推断出这一点。这就是符号计算发挥作用的地方。使用类似 SymPy 的符号计算系统，不是完全平方数的数字的平方根默认保留为未评估状态

>>> import sympy
>>> sympy.sqrt(3)
sqrt(3) 

此外——这是我们开始看到符号计算真正力量的地方——符号结果可以被符号化简。

>>> sympy.sqrt(8)
2*sqrt(2) 

更有趣的例子

上面的例子开始展示了如何使用 SymPy 精确地操作无理数。但它比那更强大。符号计算系统（顺便说一下，它们通常也被称为计算机代数系统，或者只是 CAS）如 SymPy 能够计算带有变量的符号表达式。

正如我们后面将看到的，在 SymPy 中，变量是用 symbols 定义的。与许多符号操作系统不同，SymPy 中的变量必须在使用之前定义（这个原因将在下一节中讨论）。

让我们定义一个符号表达式，表示数学表达式 (x + 2y)。

>>> from sympy import symbols
>>> x, y = symbols('x y')
>>> expr = x + 2*y
>>> expr
x + 2*y 

注意，我们写的 x + 2*y 就像我们如果 x 和 y 是普通的 Python 变量一样。但在这种情况下，而不是评估为某些东西，表达式仍然只是 x + 2*y。现在让我们玩弄一下它：

>>> expr + 1
x + 2*y + 1
>>> expr - x
2*y 

注意上面示例中的一些情况。当我们输入 expr - x 时，我们没有得到 x + 2*y - x，而只是得到 2*y。x 和 -x 自动抵消了。这类似于上面 sqrt(8) 自动变成 2*sqrt(2) 的情况。然而在 SymPy 中，并非总是这样：

>>> x*expr
x*(x + 2*y) 

在这里，我们可能希望 (x(x + 2y)) 转换成 (x² + 2xy)，但实际上我们看到表达式保持不变。这是 SymPy 中的一个常见主题。除了像 (x - x = 0) 和 (\sqrt{8} = 2\sqrt{2}) 这样的显而易见的简化外，大多数简化不会自动执行。这是因为我们可能更喜欢因式分解形式 (x(x + 2y))，或者更喜欢展开形式 (x² + 2xy)。在不同的情况下，两种形式都是有用的。在 SymPy 中，有函数可以在这两种形式之间转换。

>>> from sympy import expand, factor
>>> expanded_expr = expand(x*expr)
>>> expanded_expr
x**2 + 2*x*y
>>> factor(expanded_expr)
x*(x + 2*y) 

符号计算的威力

一个符号计算系统（如 SymPy）真正的力量在于能够进行各种各样的符号计算。SymPy 能够简化表达式、计算导数、积分和极限、解方程、处理矩阵等等，并且全部都是符号计算。它包括绘图、打印（例如数学公式的二维漂亮打印输出，或者 (\mathrm{\LaTeX})）、代码生成、物理学、统计学、组合数学、数论、几何学、逻辑学等模块。这里是 SymPy 能够实现的一小部分符号计算能力，以激发您的兴趣。

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

这将使得所有后续的示例都能使用 unicode 字符进行漂亮打印。

>>> init_printing(use_unicode=True) 

对 (\sin{(x)}e^x) 求导数。

>>> diff(sin(x)*exp(x), x)
 x           x
ℯ ⋅sin(x) + ℯ ⋅cos(x) 

计算 (\int(e^x\sin{(x)} + e^x\cos{(x)}),dx)。

>>> integrate(exp(x)*sin(x) + exp(x)*cos(x), x)
 x
ℯ ⋅sin(x) 

计算 (\int_{-\infty}^\infty \sin{(x²)},dx)。

>>> integrate(sin(x**2), (x, -oo, oo))
√2⋅√π
─────
 2 

求 (\lim_{x\to 0}\frac{\sin{(x)}}{x})。

>>> limit(sin(x)/x, x, 0)
1 

解 (x² - 2 = 0)。

>>> solve(x**2 - 2, x)
[-√2, √2] 

解微分方程 (y'' - y = e^t)。

>>> y = Function('y')
>>> dsolve(Eq(y(t).diff(t, t) - y(t), exp(t)), y(t))
 -t   ⎛     t⎞  t
y(t) = C₂⋅ℯ   + ⎜C₁ + ─⎟⋅ℯ
 ⎝     2⎠ 

求 (\left[\begin{smallmatrix}1 & 2\2 & 2\end{smallmatrix}\right]) 的特征值。

>>> Matrix([[1, 2], [2, 2]]).eigenvals()
⎧3   √17     3   √17   ⎫
⎨─ - ───: 1, ─ + ───: 1⎬
⎩2    2      2    2    ⎭ 

将贝塞尔函数 (J_{\nu}\left(z\right)) 重写为球贝塞尔函数 (j_\nu(z)) 的形式。

>>> besselj(nu, z).rewrite(jn)
√2⋅√z⋅jn(ν - 1/2, z)
────────────────────
 √π 

使用 (\mathrm{\LaTeX})，打印 (\int_{0}^{\pi} \cos^{2}{\left (x \right )}, dx)。

>>> latex(Integral(cos(x)**2, (x, 0, pi)))
\int\limits_{0}^{\pi} \cos^{2}{\left(x \right)}\, dx 

为什么选择 SymPy？

有许多计算代数系统可供选择。这篇文章列出了许多。SymPy 为何比其他选择更好呢？

首先，SymPy 是完全免费的。它是开源的，并且在宽松的 BSD 许可下发布，因此您可以修改源代码，甚至出售它。这与像 Maple 或 Mathematica 这样的流行商业系统形成对比，后者需要花费数百美元购买许可。

其次，SymPy 使用 Python。大多数计算代数系统都会发明自己的语言。但 SymPy 不同。SymPy 完全用 Python 编写，并且完全在 Python 中执行。这意味着如果您已经了解 Python，那么使用 SymPy 就更容易上手，因为您已经熟悉语法（而如果您尚未了解 Python，则学习起来也非常简单）。我们已经知道 Python 是一种经过良好设计和经受过考验的语言。SymPy 的开发者对自己编写数学软件的能力很有信心，但编程语言设计是完全不同的事情。通过重复使用现有语言，我们能够专注于那些真正重要的事情：数学。

另一个计算机代数系统 Sage 也使用 Python 作为其语言。但是 Sage 很大，下载超过一千兆字节。SymPy 的一个优点是它很轻量级。除了 Python 之外，它没有其他依赖项，因此几乎可以在任何地方轻松使用。此外，Sage 和 SymPy 的目标不同。Sage 的目标是成为一个完整的数学系统，并通过将所有主要的开源数学系统编译到一起来实现这一目标。当您在 Sage 中调用某些函数（如integrate）时，它会调用其中包含的一个开源包。实际上，SymPy 已包含在 Sage 中。另一方面，SymPy 的目标是成为一个独立的系统，所有功能都在 SymPy 本身实现。

SymPy 的另一个重要特点是可以作为一个库来使用。许多计算机代数系统侧重于在交互式环境中的可用性，但如果您希望自动化或扩展它们，这将变得困难。使用 SymPy，您既可以在交互式 Python 环境中轻松使用它，也可以在您自己的 Python 应用程序中导入它。SymPy 还提供了 API 来方便您使用自定义函数扩展它。

需要注意的地方

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

要开始，我们应该澄清一些关于 SymPy 的事情。SymPy 只是一个 Python 库，就像NumPy、Django或者 Python 标准库中的模块sys或re一样。这意味着 SymPy 并没有为 Python 语言添加任何东西。Python 语言固有的限制在 SymPy 中也同样存在。这也意味着 SymPy 尽可能使用 Python 惯用法，使得那些已经熟悉 Python 编程的人可以轻松地使用 SymPy。作为一个简单的例子，SymPy 使用 Python 语法来构建表达式。在 Python 中不允许隐式乘法（如3x或3 x），因此在 SymPy 中也不允许。要将3和x相乘，必须输入3*x，带上*。

Symbols

这个事实的一个后果是，SymPy 可以在任何可用 Python 的环境中使用。我们只需像导入任何其他库一样导入它：

>>> from sympy import * 

这将把 SymPy 中的所有函数和类导入到我们的交互式 Python 会话中。现在，假设我们开始进行计算。

>>> x + 1
Traceback (most recent call last):
...
NameError: name 'x' is not defined 

糟糕！这里发生了什么？我们试图使用变量x，但它告诉我们x未定义。在 Python 中，变量在定义之前没有意义。SymPy 也不例外。与你可能使用过的许多符号操作系统不同，在 SymPy 中，变量不会自动定义。要定义变量，我们必须使用symbols。

>>> x = symbols('x')
>>> x + 1
x + 1 

symbols接受一个由空格或逗号分隔的变量名字符串，并从中创建 Symbols。然后我们可以将这些符号赋予变量名。稍后，我们将研究一些方便的方法来解决这个问题。现在，让我们先定义最常见的变量名x、y和z，在本节的其余部分使用。

>>> x, y, z = symbols('x y z') 

最后需要注意的是，Symbol 的名称与其分配给的变量名称之间不一定有任何关系。

>>> a, b = symbols('b a')
>>> a
b
>>> b
a 

在这里，我们做了一件非常令人困惑的事情，将一个名为a的 Symbol 赋给变量b，将一个名为b的 Symbol 赋给变量a。现在名为a的 Python 变量指向名为b的 SymPy Symbol，反之亦然。多么令人困惑。我们也可以做类似的事情：

>>> crazy = symbols('unrelated')
>>> crazy + 1
unrelated + 1 

这也表明 Symbols 的名称可以比一个字符长，如果我们希望的话。

通常，最佳实践是将 Symbol 分配给同名的 Python 变量，尽管也有例外：Symbol 名称可以包含 Python 变量名称中不允许的字符，或者可能只是想通过将长名称的 Symbols 分配给单个字母 Python 变量来避免输入长名称。

为避免混淆，在本教程中，Symbol 名称和 Python 变量名称将始终相同。此外，“Symbol”一词将指代 SymPy Symbol，“variable”一词将指代 Python 变量。

最后，让我们确保我们理解 SymPy 符号和 Python 变量之间的区别。考虑以下内容：

x = symbols('x')
expr = x + 1
x = 2
print(expr) 

你认为这段代码的输出会是什么？如果你认为是3，你错了。让我们看看实际发生了什么

>>> x = symbols('x')
>>> expr = x + 1
>>> x = 2
>>> print(expr)
x + 1 

将x更改为2对expr没有影响。这是因为x = 2将 Python 变量x更改为2，但对创建expr时使用的 SymPy 符号x没有影响。当我们创建expr时，Python 变量x是一个 Symbol。在创建后，我们将 Python 变量x更改为 2。但expr保持不变。这种行为不是 SymPy 特有的。所有 Python 程序都是这样工作的：如果一个变量被更改，已经使用该变量创建的表达式不会自动更改。例如

>>> x = 'abc'
>>> expr = x + 'def'
>>> expr
'abcdef'
>>> x = 'ABC'
>>> expr
'abcdef' 

在这个例子中，如果我们想知道expr在新值x下是什么，我们需要重新评估创建expr的代码，即expr = x + 1。如果有几行创建了expr，这可能会很复杂。使用像 SymPy 这样的符号计算系统的一个优点是我们可以为expr构建一个符号表示，然后用值替换x。在 SymPy 中正确的方法是使用subs，稍后将更详细讨论。

>>> x = symbols('x')
>>> expr = x + 1
>>> expr.subs(x, 2)
3 
```  ## 等号

另一个非常重要的结果是，SymPy 不扩展 Python 语法的事实是，`=`在 SymPy 中不表示相等。而是 Python 变量赋值。这是硬编码到 Python 语言中的，SymPy 不会尝试改变这一点。

然而，你可能认为在 Python 中用于相等测试的`==`也用于 SymPy 作为相等。这也不完全正确。让我们看看当我们使用`==`时会发生什么。

```py
>>> x + 1 == 4
False 

代替将x + 1 == 4象征性地处理，我们只是得到了False。在 SymPy 中，==表示精确的结构相等测试。这意味着a == b意味着我们在询问是否 (a = b)。我们始终得到bool作为==的结果。有一个单独的对象，称为Eq，可以用来创建符号相等性。

>>> Eq(x + 1, 4)
Eq(x + 1, 4) 

关于==还有一个额外的警告。假设我们想知道是否 ((x + 1)² = x² + 2x + 1)。我们可能会尝试这样做：

>>> (x + 1)**2 == x**2 + 2*x + 1
False 

我们再次得到了False。然而，((x + 1)²)确实等于(x² + 2x + 1)。这里发生了什么？我们在 SymPy 中找到了一个错误吗，还是它只是不能识别这个基本的代数事实？

从上面回顾，==表示精确的结构相等测试。“精确”在这里意味着只有两个表达式在结构上完全相等时才会用==比较相等。在这里，((x + 1)²)和(x² + 2x + 1)在结构上不相同。一个是两项相加的幂，另一个是三项相加。

结果表明，在作为库使用 SymPy 时，将==用于精确的结构相等性比将其用于表示符号相等性或进行数学相等性的检测更加有用。然而，作为新用户，您可能更关心后两者。我们已经看到了表示等式的另一种选择，Eq。要测试两个事物是否相等，最好记住一个基本事实，即如果(a = b)，那么(a - b = 0)。因此，检查(a = b)的最佳方法是取(a - b)并简化它，看看它是否变为 0。我们将在后面学习到，执行此操作的函数称为simplify。这种方法并非万无一失——事实上，可以从理论上证明无法确定一般情况下两个符号表达式是否完全相等——但对于大多数常见的表达式，它效果非常好。

>>> a = (x + 1)**2
>>> b = x**2 + 2*x + 1
>>> simplify(a - b)
0
>>> c = x**2 - 2*x + 1
>>> simplify(a - c)
4*x 

还有一种叫做equals的方法，通过在随机点数值上评估它们来测试两个表达式是否相等。

>>> a = cos(x)**2 - sin(x)**2
>>> b = cos(2*x)
>>> a.equals(b)
True 
```  ## 两个最后注意事项：`^` 和 `/`

您可能已经注意到，我们一直在使用`**`来表示乘方，而不是标准的`^`。这是因为 SymPy 遵循 Python 的惯例。在 Python 中，`^` 表示逻辑异或。SymPy 也遵循了这一惯例：

```py
>>> True ^ False
True
>>> True ^ True
False
>>> Xor(x, y)
x ^ y 

最后，需要对 SymPy 的工作原理进行一点技术性讨论。当您键入类似x + 1的表达式时，SymPy 的符号x会与 Python 的整数1相加。然后 Python 的操作规则允许 SymPy 告诉 Python SymPy 对象知道如何与 Python 整数相加，因此1会自动转换为 SymPy 的整数对象。

这种运算符的魔术是在幕后自动发生的，你很少需要知道它正在发生。然而，也有一个例外。每当您结合一个 SymPy 对象和一个 SymPy 对象，或者一个 SymPy 对象和一个 Python 对象时，您会得到一个 SymPy 对象，但是每当您结合两个 Python 对象时，SymPy 从不参与其中，因此您会得到一个 Python 对象。

>>> type(Integer(1) + 1)
<class 'sympy.core.numbers.Integer'>
>>> type(1 + 1)
<... 'int'> 

这通常不是什么大问题。Python 的整数与 SymPy 的整数工作方式基本相同，但有一个重要的例外：除法。在 SymPy 中，两个整数的除法会得到一个有理数：

>>> Integer(1)/Integer(3)
1/3
>>> type(Integer(1)/Integer(3))
<class 'sympy.core.numbers.Rational'> 

但是在 Python 中，/表示整数除法或浮点数除法，具体取决于您使用的是 Python 2 还是 Python 3，以及是否在 Python 2 中运行了from __future__ import division，这在 SymPy 1.5.1 以上的版本中不再支持：

>>> from __future__ import division
>>> 1/2
0.5 

为了避免这种情况，我们可以显式地构造有理数对象。

>>> Rational(1, 2)
1/2 

当我们在一个较大的符号表达式中遇到带有int/int的情况时，也会出现这个问题。例如：

>>> x + 1/2
x + 0.5 

这是因为 Python 首先将1/2计算为0.5，然后在将其与x相加时将其转换为 SymPy 类型。同样，我们可以通过显式创建有理数来避免这种情况：

>>> x + Rational(1, 2)
x + 1/2 

在陷阱与风险文档中有几个避免这种情况的建议。

进一步阅读

关于本节涵盖的主题的更多讨论，请参阅陷阱与风险。

SymPy 特性

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

本节讨论常见和高级的 SymPy 操作和特性。

内容

• 基本操作

  • 替换

  • 将字符串转换为 SymPy 表达式

  • evalf

  • lambdify

• 打印

  • 打印机

  • 设置漂亮打印

  • 打印函数

• 简化

  • simplify

  • 多项式/有理函数简化

  • 三角函数简化

  • 幂运算

  • 指数和对数

  • 特殊函数

  • 示例：连分数

• 微积分

  • 导数

  • 积分

  • 数值积分

  • 极限

  • 级数展开

  • 有限差分

• 求解器

  • 关于方程的注释

  • 代数方程的求解

  • 微分方程的求解

• 矩阵

  • 基本操作

  • 基本方法

  • 矩阵构造器

  • 高级方法

  • 可能的问题

• 高级表达式操作

  • 理解表达式树

  • 递归遍历表达式树

  • 防止表达式求值

基本操作

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

在 SymPy 中，我们将讨论表达式操作所需的一些最基本操作。稍后将在高级表达式操作部分讨论一些更高级的操作。

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

替换

数学表达式中最常见的操作之一是替换。替换将表达式中的某些内容替换为其他内容。可以使用subs方法来完成。例如

>>> expr = cos(x) + 1
>>> expr.subs(x, y)
cos(y) + 1 

替换通常出于两个原因之一进行：

1. 在某一点评估一个表达式。例如，如果我们的表达式是cos(x) + 1，我们想在点x = 0处评估它，这样我们得到cos(0) + 1，结果是 2。

   >>> expr.subs(x, 0)
   2 
   

2. 将一个子表达式替换为另一个子表达式。我们可能有两个原因想要这样做。第一个是如果我们试图构建具有某种对称性的表达式，比如 (x^{x{x^x}})。为了构建这个表达式，我们可以从x**y开始，然后用x**y替换y。然后我们会得到x**(x**y)。如果我们在这个新表达式中用x**x替换y，我们会得到x**(x**(x**x))，这就是我们想要的表达式。

   >>> expr = x**y
   >>> expr
   x**y
   >>> expr = expr.subs(y, x**y)
   >>> expr
   x**(x**y)
   >>> expr = expr.subs(y, x**x)
   >>> expr
   x**(x**(x**x)) 
   

   第二个原因是，如果我们想进行非常精确的简化，或者可能是 SymPy 无法做的简化。例如，假设我们有 (\sin(2x) + \cos(2x))，我们想要用 (\sin(2x)) 替换为 (2\sin(x)\cos(x))。正如我们稍后将学习的，函数 expand_trig 就是这么做的。然而，这个函数也会展开 (\cos(2x))，这可能不是我们想要的。虽然有方法可以进行如此精确的简化，我们将在高级表达式操作部分学习其中一些，一个简单的方法就是直接用 (2\sin(x)\cos(x)) 替换 (\sin(2x))。

   >>> expr = sin(2*x) + cos(2*x)
   >>> expand_trig(expr)
   2*sin(x)*cos(x) + 2*cos(x)**2 - 1
   >>> expr.subs(sin(2*x), 2*sin(x)*cos(x))
   2*sin(x)*cos(x) + cos(2*x) 
   

关于subs有两个重要的事项需要注意。首先，它返回一个新的表达式。SymPy 对象是不可变的。这意味着subs不会就地修改它。例如

>>> expr = cos(x)
>>> expr.subs(x, 0)
1
>>> expr
cos(x)
>>> x
x 

在这里，我们看到执行expr.subs(x, 0)不会改变expr。实际上，由于 SymPy 表达式是不可变的，没有函数会就地修改它们。所有函数都会返回新的表达式。

要一次执行多个替换，请将 (old, new) 对的列表传递给 subs。

>>> expr = x**3 + 4*x*y - z
>>> expr.subs([(x, 2), (y, 4), (z, 0)])
40 

将这个与列表推导结合起来，可以做大量类似的替换。例如，假设我们有 (x⁴ - 4x³ + 4x² - 2x + 3)，我们想要替换所有偶次幂为 (y) 的 (x) 实例，得到 (y⁴ - 4x³ + 4y² - 2x + 3)。

>>> expr = x**4 - 4*x**3 + 4*x**2 - 2*x + 3
>>> replacements = [(x**i, y**i) for i in range(5) if i % 2 == 0]
>>> expr.subs(replacements)
-4*x**3 - 2*x + y**4 + 4*y**2 + 3 

将字符串转换为 SymPy 表达式

函数 sympify（注意是 sympify，不要与 simplify 混淆）可用于将字符串转换为 SymPy 表达式。

例如

>>> str_expr = "x**2 + 3*x - 1/2"
>>> expr = sympify(str_expr)
>>> expr
x**2 + 3*x - 1/2
>>> expr.subs(x, 2)
19/2 

警告

sympify 使用 eval。不要对未经过滤的输入使用它。

evalf

要将数值表达式评估为浮点数，请使用 evalf。

>>> expr = sqrt(8)
>>> expr.evalf()
2.82842712474619 

SymPy 可以将浮点表达式计算到任意精度。默认使用 15 位数字精度，但您可以将任何数字作为 evalf 的参数传递。让我们计算 (\pi) 的前 100 位小数。

>>> pi.evalf(100)
3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068 

要在点处使用符号评估表达式的数值表达式，我们可以使用 subs，然后使用 evalf，但将替换传递给 evalf 使用 subs 标志更有效和数值稳定，该标志接受 Symbol: point 对的字典。

>>> expr = cos(2*x)
>>> expr.evalf(subs={x: 2.4})
0.0874989834394464 

有时，在求值表达式后，可能会保留比所需精度小的舍入误差。可以通过将 chop 标志设置为 True，由用户自行决定是否移除这些数字。

>>> one = cos(1)**2 + sin(1)**2
>>> (one - 1).evalf()
-0.e-124
>>> (one - 1).evalf(chop=True)
0 

lambdify

如果您打算在许多点评估表达式，则 subs 和 evalf 很好用。但如果您打算在一千个点评估表达式，则使用 SymPy 比必要的更慢，特别是如果您只关心机器精度。相反，您应该使用像 NumPy 和 SciPy 这样的库。

将 SymPy 表达式转换为可以进行数值评估的表达式的最简单方法是使用 lambdify 函数。lambdify 的功能类似于 lambda 函数，但它将 SymPy 名称转换为给定数值库（通常是 NumPy）的名称。例如

>>> import numpy 
>>> a = numpy.arange(10) 
>>> expr = sin(x)
>>> f = lambdify(x, expr, "numpy") 
>>> f(a) 
[ 0\.          0.84147098  0.90929743  0.14112001 -0.7568025  -0.95892427
 -0.2794155   0.6569866   0.98935825  0.41211849] 

警告

lambdify 使用 eval。不要在未经过消毒处理的输入上使用它。

您可以使用除 NumPy 之外的其他库。例如，要使用标准库 math 模块，请使用 "math"。

>>> f = lambdify(x, expr, "math")
>>> f(0.1)
0.0998334166468 

要与 lambdify 不了解的数值库一起使用 lambdify，请传递一个 sympy_name:numerical_function 对的字典。例如

>>> def mysin(x):
...  """
...     My sine. Note that this is only accurate for small x.
...     """
...     return x
>>> f = lambdify(x, expr, {"sin":mysin})
>>> f(0.1)
0.1 

打印（Printing）

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

正如我们已经看到的，SymPy 可以使用 Unicode 字符对其输出进行漂亮打印。 这是关于 SymPy 中最常见的打印选项的简短介绍。

打印机

SymPy 提供了几种打印方式。其中最常见的是

• 字符串（str）

• srepr

• ASCII 漂亮打印器

• Unicode 漂亮打印器

• LaTeX

• MathML

• 点（Dot）

除了这些外，还有可以将 SymPy 对象输出为代码的“打印机”，例如 C、Fortran、Javascript、Theano 和 Python。 这些在本教程中不讨论。

设置漂亮打印

如果您只希望得到最佳的漂亮打印效果，请使用 init_printing() 函数。 这将自动启用您环境中可用的最佳打印机。

>>> from sympy import init_printing
>>> init_printing() 

如果您计划在交互式计算器类型的会话中工作，则 init_session() 函数将自动导入 SymPy 中的所有内容，创建一些常见符号，设置绘图，并运行 init_printing()。

>>> from sympy import init_session

>>> init_session() 

Python console for SymPy 0.7.3 (Python 2.7.5-64-bit) (ground types: gmpy)



These commands were executed:

>>> from __future__ import division

>>> from sympy import *

>>> x, y, z, t = symbols('x y z t')

>>> k, m, n = symbols('k m n', integer=True)

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

>>> init_printing() # doctest: +SKIP



Documentation can be found at https://www.sympy.org/ 

>>> 

无论如何，都会发生以下情况：

• 在 IPython QTConsole 中，如果安装了 (\mathrm{\LaTeX})，它将启用使用 (\mathrm{\LaTeX}) 的打印机。

  

  如果未安装 (\mathrm{\LaTeX})，但安装了 Matplotlib，则将使用 Matplotlib 渲染引擎。 如果未安装 Matplotlib，则使用 Unicode 漂亮打印器。

• 在 IPython 笔记本中，它将使用 MathJax 渲染 (\mathrm{\LaTeX})。

  

• 在 IPython 控制台会话或常规 Python 会话中，如果终端支持 Unicode，则将使用 Unicode 漂亮打印器。

  

• 在不支持 Unicode 的终端上，将使用 ASCII 漂亮打印器。

  

若要显式禁用 (\mathrm{\LaTeX})，请向 init_printing() 或 init_session() 传递 use_latex=False。 若要显式禁用 Unicode，请传递 use_unicode=False。

打印函数（Printing Functions）

除了自动打印外，您还可以通过调用相应的函数显式地使用任何一个打印机。

字符串（str）

要获得表达式的字符串形式，请使用 str(expr)。 这也是由 print(expr) 生成的形式。 字符串形式设计为易于阅读，但形式上正确，以便可以复制和粘贴。 表达式的 str() 形式通常看起来与输入它时的表达式完全相同。

>>> from sympy import *
>>> x, y, z = symbols('x y z')
>>> str(Integral(sqrt(1/x), x))
'Integral(sqrt(1/x), x)'
>>> print(Integral(sqrt(1/x), x))
Integral(sqrt(1/x), x) 

srepr

表达式的 srepr 形式旨在显示表达式的确切形式。将在高级表达式操作部分进一步讨论。要获取它，请使用 srepr() [1]。

>>> srepr(Integral(sqrt(1/x), x))
"Integral(Pow(Pow(Symbol('x'), Integer(-1)), Rational(1, 2)), Tuple(Symbol('x')))" 

srepr 形式在大多数情况下对理解表达式内部结构很有用。

ASCII 漂亮打印机

ASCII 漂亮打印机可以从 pprint() 访问。如果终端不支持 Unicode，则默认使用 ASCII 打印机。否则，必须传递 use_unicode=False。

>>> pprint(Integral(sqrt(1/x), x), use_unicode=False)
 /
 |
 |     ___
 |    / 1
 |   /  -  dx
 | \/   x
 |
/ 

pprint() 将输出打印到屏幕。如果需要字符串形式，请使用 pretty()。

>>> pretty(Integral(sqrt(1/x), x), use_unicode=False)
'  /          \n |           \n |     ___   \n |    / 1    \n |   /  -  dx\n | \\/   x    \n |           \n/            '
>>> print(pretty(Integral(sqrt(1/x), x), use_unicode=False))
 /
 |
 |     ___
 |    / 1
 |   /  -  dx
 | \/   x
 |
/ 

Unicode 漂亮打印机

Unicode 漂亮打印机也可以从 pprint() 和 pretty() 访问。如果终端支持 Unicode，则会自动使用。如果 pprint() 无法检测到终端支持 Unicode，则可以传递 use_unicode=True 强制使用 Unicode。

>>> pprint(Integral(sqrt(1/x), x), use_unicode=True)
⌠
⎮     ___
⎮    ╱ 1
⎮   ╱  ─  dx
⎮ ╲╱   x
⌡ 

(\mathrm{\LaTeX})

要获得表达式的 (\mathrm{\LaTeX}) 形式，请使用 latex()。

>>> print(latex(Integral(sqrt(1/x), x)))
\int \sqrt{\frac{1}{x}}\, dx 

latex() 函数有许多选项，可以改变不同事物的格式。详见其文档，获取更多细节。

MathML

还有一个名为 print_mathml() 的 MathML 打印机。它必须从 sympy.printing.mathml 导入。

>>> from sympy.printing.mathml import print_mathml
>>> print_mathml(Integral(sqrt(1/x), x))
<apply>
 <int/>
 <bvar>
 <ci>x</ci>
 </bvar>
 <apply>
 <root/>
 <apply>
 <power/>
 <ci>x</ci>
 <cn>-1</cn>
 </apply>
 </apply>
</apply> 

print_mathml() 将输出打印出来。如果需要字符串，请使用函数 mathml()。

点

sympy.printing.dot 中的 dotprint() 函数将输出打印到 dot 格式，可以使用 Graphviz 渲染。参见高级表达式操作部分，了解一些此打印机输出的示例。

这里是 dotprint() 函数的原始输出示例。

>>> from sympy.printing.dot import dotprint
>>> from sympy.abc import x
>>> print(dotprint(x+2))
digraph{

# Graph style
"ordering"="out"
"rankdir"="TD"

#########
# Nodes #
#########

"Add(Integer(2), Symbol('x'))_()" ["color"="black", "label"="Add", "shape"="ellipse"];
"Integer(2)_(0,)" ["color"="black", "label"="2", "shape"="ellipse"];
"Symbol('x')_(1,)" ["color"="black", "label"="x", "shape"="ellipse"];

#########
# Edges #
#########

"Add(Integer(2), Symbol('x'))_()" -> "Integer(2)_(0,)";
"Add(Integer(2), Symbol('x'))_()" -> "Symbol('x')_(1,)";
} 

脚注

简化

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

为了使这份文档更易读，我们将启用漂亮的打印输出。

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

simplify

现在让我们跳进来做一些有趣的数学。符号操作系统最有用的功能之一是简化数学表达式的能力。SymPy 有几十个函数来执行各种简化。还有一个名为简化()的通用函数，它尝试以智能方式应用所有这些函数，以得到表达式的最简形式。以下是一些示例

>>> simplify(sin(x)**2 + cos(x)**2)
1
>>> simplify((x**3 + x**2 - x - 1)/(x**2 + 2*x + 1))
x - 1
>>> simplify(gamma(x)/gamma(x - 2))
(x - 2)⋅(x - 1) 

在这里，gamma(x)是(\Gamma(x))，即伽玛函数。我们看到简化()能够处理大量的表达式类。

但简化()存在一个陷阱。它只是应用 SymPy 中的所有主要简化操作，并使用启发式方法来确定最简结果。但是，“最简”并不是一个明确定义的术语。例如，假设我们想将(x² + 2x + 1)“简化”为((x + 1)²)：

>>> simplify(x**2 + 2*x + 1)
 2
x  + 2⋅x + 1 

我们没有得到我们想要的。有一个执行这种简化的函数，称为factor()，将在下面讨论。

另一个简化()的陷阱是它可能会不必要地慢，因为它在选择最佳选项之前尝试了许多种简化方法。如果您已经确切地知道需要哪种简化，最好应用特定的简化函数来执行这些简化。

使用特定的简化函数而不是简化()还有一个优点，即特定函数对其输出的形式有一定的保证。这些将在每个函数下面讨论。例如，当在具有有理系数的多项式上调用factor()时，保证将多项式因式分解为不可约因子。简化()没有任何保证。它完全是启发式的，并且如上所示，甚至可能会忽略 SymPy 能够执行的一种可能简化类型。

简化()在交互式使用时效果最佳，当您只想将表达式简化为更简单的形式时。然后，一旦看到简化()的返回结果，您可以选择应用特定函数以获得更精确的结果。当您不知道表达式将采用什么形式时，这也是一个有用的通用函数来简化它。

多项式/有理函数简化

expand

expand()是 SymPy 中最常见的简化函数之一。虽然它有很多用途，但现在我们将考虑它在展开多项式表达式中的功能。例如：

>>> expand((x + 1)**2)
 2
x  + 2⋅x + 1
>>> expand((x + 2)*(x - 3))
 2
x  - x - 6 

给定多项式，expand()将其放入一组单项式的标准形式。

expand() 可能听起来不像一个简化函数。毕竟，从其名称来看，它使表达式变大而不是变小。通常情况下确实如此，但经常调用 expand() 后表达式会因为约简而变小。

>>> expand((x + 1)*(x - 2) - (x - 1)*x)
-2 

factor

factor() 接受一个多项式并在有理数域上将其因式分解为不可约因子。例如：

>>> factor(x**3 - x**2 + x - 1)
 ⎛ 2    ⎞
(x - 1)⋅⎝x  + 1⎠
>>> factor(x**2*z + 4*x*y*z + 4*y**2*z)
 2
z⋅(x + 2⋅y) 

对于多项式，factor() 是 expand() 的反函数。factor() 使用完整的多变量因式分解算法在有理数域上运行，这意味着 factor() 返回的每个因子都保证是不可约的。

如果你对因子本身感兴趣，factor_list 返回一个更结构化的输出。

>>> factor_list(x**2*z + 4*x*y*z + 4*y**2*z)
(1, [(z, 1), (x + 2⋅y, 2)]) 

请注意，factor 和 expand 的输入不一定严格是多项式。它们会智能地因式分解或展开任何类型的表达式（尽管如果输入不再是有理数域上的多项式，则因子可能不是不可约的）。

>>> expand((cos(x) + sin(x))**2)
 2                           2
sin (x) + 2⋅sin(x)⋅cos(x) + cos (x)
>>> factor(cos(x)**2 + 2*cos(x)*sin(x) + sin(x)**2)
 2
(sin(x) + cos(x)) 

collect

collect() 在表达式中收集一个项的公共幂次。例如

>>> expr = x*y + x - 3 + 2*x**2 - z*x**2 + x**3
>>> expr
 3    2        2
x  - x ⋅z + 2⋅x  + x⋅y + x - 3
>>> collected_expr = collect(expr, x)
>>> collected_expr
 3    2
x  + x ⋅(2 - z) + x⋅(y + 1) - 3 

collect() 特别在与 .coeff() 方法一起使用时非常有用。expr.coeff(x, n) 给出 expr 中 x**n 的系数：

>>> collected_expr.coeff(x, 2)
2 - z 

cancel

cancel() 将任何有理函数放入标准的规范形式 (\frac{p}{q})，其中 (p) 和 (q) 是没有公因式的展开多项式，并且 (p) 和 (q) 的首项系数没有分母（即为整数）。

>>> cancel((x**2 + 2*x + 1)/(x**2 + x))
x + 1
─────
 x 

>>> expr = 1/x + (3*x/2 - 2)/(x - 4)
>>> expr
3⋅x
─── - 2
 2        1
─────── + ─
 x - 4    x
>>> cancel(expr)
 2
3⋅x  - 2⋅x - 8
──────────────
 2
 2⋅x  - 8⋅x 

>>> expr = (x*y**2 - 2*x*y*z + x*z**2 + y**2 - 2*y*z + z**2)/(x**2 - 1)
>>> expr
 2                2    2            2
x⋅y  - 2⋅x⋅y⋅z + x⋅z  + y  - 2⋅y⋅z + z
───────────────────────────────────────
 2
 x  - 1
>>> cancel(expr)
 2            2
y  - 2⋅y⋅z + z
───────────────
 x - 1 

请注意，由于 factor() 将完全因式分解表达式的分子和分母，因此它也可以用来做同样的事情：

>>> factor(expr)
 2
(y - z)
────────
 x - 1 

然而，如果你只关心表达式是否处于约简形式，cancel() 比 factor() 更有效。

apart

apart() 在有理函数上执行偏分数分解。

>>> expr = (4*x**3 + 21*x**2 + 10*x + 12)/(x**4 + 5*x**3 + 5*x**2 + 4*x)
>>> expr
 3       2
4⋅x  + 21⋅x  + 10⋅x + 12
────────────────────────
 4      3      2
 x  + 5⋅x  + 5⋅x  + 4⋅x
>>> apart(expr)
 2⋅x - 1       1     3
────────── - ───── + ─
 2           x + 4   x
x  + x + 1 

三角函数简化

注意

SymPy 遵循 Python 的反三角函数命名约定，即在函数名前加上 a。例如，反余弦或弧余弦称为 acos()。

>>> acos(x)
acos(x)
>>> cos(acos(x))
x
>>> asin(1)
π
─
2 

trigsimp

要使用三角函数恒等式简化表达式，请使用 trigsimp()。

>>> trigsimp(sin(x)**2 + cos(x)**2)
1
>>> trigsimp(sin(x)**4 - 2*cos(x)**2*sin(x)**2 + cos(x)**4)
cos(4⋅x)   1
──────── + ─
 2       2
>>> trigsimp(sin(x)*tan(x)/sec(x))
 2
sin (x) 

trigsimp() 也适用于双曲三角函数。

>>> trigsimp(cosh(x)**2 + sinh(x)**2)
cosh(2⋅x)
>>> trigsimp(sinh(x)/tanh(x))
cosh(x) 

与 simplify() 类似，trigsimp() 对输入表达式应用各种三角函数恒等式，然后使用启发式算法返回“最佳”表达式。

expand_trig

要展开三角函数，即应用和角或双角恒等式，请使用 expand_trig()。

>>> expand_trig(sin(x + y))
sin(x)⋅cos(y) + sin(y)⋅cos(x)
>>> expand_trig(tan(2*x))
 2⋅tan(x)
───────────
 2
1 - tan (x) 

因为 expand_trig() 倾向于使三角函数表达式变大，而 trigsimp() 倾向于使其变小，所以可以用 trigsimp() 反向应用这些恒等式。

>>> trigsimp(sin(x)*cos(y) + sin(y)*cos(x))
sin(x + y) 

幂次

在我们介绍幂简化函数之前，有必要就指数满足的三种恒等式进行数学讨论。

1. (x^axb = x^{a + b})

2. (x^aya = (xy)^a)

3. ((x^a)b = x^{ab})

标识 1 总是成立。

标识 2 并非总是成立。例如，如果 (x = y = -1) 且 (a = \frac{1}{2})，则 (x^aya = \sqrt{-1}\sqrt{-1} = i\cdot i = -1)，而 ((xy)^a = \sqrt{-1\cdot-1} = \sqrt{1} = 1)。然而，标识 2 至少在 (x) 和 (y) 非负且 (a) 是实数时成立（它可能在其他条件下也成立）。标识 2 失败的常见后果是 (\sqrt{x}\sqrt{y} \neq \sqrt{xy})。

标识 3 并非总是成立。例如，如果 (x = -1)，(a = 2)，(b = \frac{1}{2})，则 ((x^a)b = {\left((-1)²\right)}^{1/2} = \sqrt{1} = 1)，而 (x^{ab} = (-1)^{2\cdot1/2} = (-1)¹ = -1)。然而，当 (b) 是整数时，标识 3 成立（同样，在其他情况下也可能成立）。标识 3 失败的两个常见后果是 (\sqrt{x²}\neq x) 和 (\sqrt{\frac{1}{x}} \neq \frac{1}{\sqrt{x}})。

总结

标识	满足条件的充分条件	条件不满足时的反例	重要后果

|

1. (x^axb = x^{a + b})

总是成立	无	无

|

1. (x^aya = (xy)^a)

(x, y \geq 0) 和 (a \in \mathbb{R})	((-1){1/2}(-1) \neq (-1\cdot-1)^{1/2})	(\sqrt{x}\sqrt{y} \neq \sqrt{xy}) 在一般情况下

|

1. ((x^a)b = x^{ab})

(b \in \mathbb{Z})	({\left((-1)²\right)}^{1/2} \neq (-1)^{2\cdot1/2})	(\sqrt{x²}\neq x) 和 (\sqrt{\frac{1}{x}}\neq\frac{1}{\sqrt{x}}) 在一般情况下

这一点很重要，因为默认情况下，如果在一般情况下简化不成立，SymPy 将不会执行这些简化。

为了使 SymPy 进行涉及仅在某些假设下才成立的简化，我们需要在符号上加上假设。稍后我们将对假设系统进行全面讨论，但现在我们只需知道以下内容。

• 默认情况下，SymPy 符号被假定为复数（属于 (\mathbb{C}) 的元素）。也就是说，除非对所有复数成立，否则不会对具有给定符号的表达式进行简化。

• 通过将假设传递给 symbols()，可以为符号提供不同的假设。在本节的其余部分中，我们假设 x 和 y 是正数，a 和 b 是实数。我们将保留 z、t 和 c 作为任意复杂符号，以演示在这种情况下会发生什么。

  >>> x, y = symbols('x y', positive=True)
  >>> a, b = symbols('a b', real=True)
  >>> z, t, c = symbols('z t c') 
  

注意

在 SymPy 中，sqrt(x) 只是 x**Rational(1, 2) 的一个快捷方式。它们是完全相同的对象。

>>> sqrt(x) == x**Rational(1, 2)
True 

powsimp

powsimp() 应用上述标识的 1 和 2 号，从左到右。

>>> powsimp(x**a*x**b)
 a + b
 x
>>> powsimp(x**a*y**a)
 a
(x⋅y) 

注意，如果简化不成立，powsimp() 将拒绝执行该简化。

>>> powsimp(t**c*z**c)
 c  c
t ⋅z 

如果您知道要应用这种简化，但不想涉及假设，可以传递 force=True 标志。这将强制进行简化，而不考虑假设。

>>> powsimp(t**c*z**c, force=True)
 c
(t⋅z) 

请注意，在某些情况下，特别是当指数是整数或有理数，并且恒等式 2 成立时，它将自动应用。

>>> (z*t)**2
 2  2
 t ⋅z
>>> sqrt(x*y)
 √x⋅√y 

这意味着使用 powsimp() 将无法撤销此恒等式，因为即使 powsimp() 将基数放在一起，它们也会自动再次分开。

>>> powsimp(z**2*t**2)
 2  2
 t ⋅z
>>> powsimp(sqrt(x)*sqrt(y))
 √x⋅√y 

expand_power_exp / expand_power_base

expand_power_exp() 和 expand_power_base() 分别从右到左应用恒等式 1 和 2。

>>> expand_power_exp(x**(a + b))
 a  b
x ⋅x 

>>> expand_power_base((x*y)**a)
 a  a
x ⋅y 

与 powsimp() 类似，如果不成立，则不会应用恒等式 2。

>>> expand_power_base((z*t)**c)
 c
(t⋅z) 

与 powsimp() 类似，您可以通过使用 force=True 来强制发生扩展，而无需操作假设。

>>> expand_power_base((z*t)**c, force=True)
 c  c
 t ⋅z 

与恒等式 2 一样，如果幂是一个数，则恒等式 1 会自动应用，因此无法通过expand_power_exp()来撤销。

>>> x**2*x**3
 5
 x
>>> expand_power_exp(x**5)
 5
 x 

powdenest

powdenest()从左到右应用恒等式 3。

>>> powdenest((x**a)**b)
 a⋅b
x 

与之前一样，如果在给定的假设下恒等式不成立，则不会应用该恒等式。

>>> powdenest((z**a)**b)
 b
⎛ a⎞
⎝z ⎠ 

同样，可以通过force=True手动覆盖这一点。

>>> powdenest((z**a)**b, force=True)
 a⋅b
z 

指数和对数

注意

在 SymPy 中，就像在 Python 和大多数编程语言中一样，log 是自然对数，也称为 ln。SymPy 自动提供 ln = log 的别名以防您忘记这一点。

>>> ln(x)
log(x) 

对数与幂具有类似的问题。主要有两个恒等式

1. (\log{(xy)} = \log{(x)} + \log{(y)})

2. (\log{(x^n)} = n\log{(x)})

由于复数对数的分支切割，恒等式对任意复数 (x) 和 (y) 都不成立。但是，如果 (x) 和 (y) 是正数，并且 (n) 是实数，则恒等式成立的充分条件。

>>> x, y = symbols('x y', positive=True)
>>> n = symbols('n', real=True) 

与之前一样，z 和 t 将是没有额外假设的符号。

请注意，恒等式 (\log{\left(\frac{x}{y}\right)} = \log(x) - \log(y)) 是恒等式 1 和 2 的一个特殊情况，由 (\log{\left(\frac{x}{y}\right)} =) (\log{\left(x\cdot\frac{1}{y}\right)} =) (\log(x) + \log{\left( y^{-1}\right)} =) (\log(x) - \log(y)) 得出，因此如果 (x) 和 (y) 是正数，则它也成立，但不一定在一般情况下成立。

我们还看到 (\log{\left( e^x \right)} = x) 来自 (\log{\left( e^x \right)} = x\log(e) = x)，因此在 (x) 是实数时成立（并且可以验证它对于任意复数 (x) 并不总是成立，例如，(\log{\left(e^{x + 2\pi i}\right)} = \log{\left(e^x\right)} = x \neq x + 2\pi i))。

expand_log

要从左到右应用恒等式 1 和 2，请使用 expand_log()。除非它们有效，否则恒等式不会应用。

>>> expand_log(log(x*y))
log(x) + log(y)
>>> expand_log(log(x/y))
log(x) - log(y)
>>> expand_log(log(x**2))
2⋅log(x)
>>> expand_log(log(x**n))
n⋅log(x)
>>> expand_log(log(z*t))
log(t⋅z) 

与 powsimp() 和 powdenest() 一样，expand_log() 也有一个 force 选项，可以用于忽略假设。

>>> expand_log(log(z**2))
 ⎛ 2⎞
log⎝z ⎠
>>> expand_log(log(z**2), force=True)
2⋅log(z) 

logcombine

要从右到左应用恒等式 1 和 2，请使用 logcombine()。

>>> logcombine(log(x) + log(y))
log(x⋅y)
>>> logcombine(n*log(x))
 ⎛ n⎞
log⎝x ⎠
>>> logcombine(n*log(z))
n⋅log(z) 

logcombine()还有一个force选项，可用于忽略假设。

>>> logcombine(n*log(z), force=True)
 ⎛ n⎞
log⎝z ⎠ 

特殊函数

SymPy 实现了数十个特殊函数，涵盖从组合数学到数学物理的各种函数。

SymPy 包含的特殊函数及其文档的详细列表位于函数模块页面。

为了本教程的目的，让我们介绍 SymPy 中的一些特殊函数。

让我们将x、y和z定义为常规复数符号，去除我们在前一节中对它们的任何假设。我们还将定义k、m和n。

>>> x, y, z = symbols('x y z')
>>> k, m, n = symbols('k m n') 

阶乘函数是factorial。factorial(n)代表(n!= 1\cdot2\cdots(n - 1)\cdot n)。(n!)表示(n)个不同项目的排列数。

>>> factorial(n)
n! 

二项式系数函数是binomial。binomial(n, k)代表(\binom{n}{k})，即从(n)个不同项目中选择(k)个项目的方法数。它通常写作(nCk)，发音为“(n) choose (k)”。

>>> binomial(n, k)
⎛n⎞
⎜ ⎟
⎝k⎠ 

阶乘函数与伽玛函数密切相关，gamma。gamma(z)表示(\Gamma(z) = \int_0^\infty t^{z - 1}e^{-t},dt)，对于正整数(z)，与((z - 1)!)相同。

>>> gamma(z)
Γ(z) 

广义超几何函数是hyper。hyper([a_1, ..., a_p], [b_1, ..., b_q], z)表示({}_pF_q\left(\begin{matrix} a_1, \cdots, a_p \ b_1, \cdots, b_q \end{matrix} \middle| z \right))。最常见的情况是({}_2F_1)，通常称为普通超几何函数。

>>> hyper([1, 2], [3], z)
 ┌─  ⎛1, 2 │  ⎞
 ├─  ⎜     │ z⎟
2╵ 1 ⎝ 3   │  ⎠ 

重写

处理特殊函数的常见方法是将它们重写为彼此的函数。这适用于 SymPy 中的任何函数，而不仅仅是特殊函数。要将表达式重写为函数形式，请使用expr.rewrite(function)。例如，

>>> tan(x).rewrite(cos)
 ⎛    π⎞
cos⎜x - ─⎟
 ⎝    2⎠
──────────
 cos(x)
>>> factorial(x).rewrite(gamma)
Γ(x + 1) 

关于应用更有针对性的重写的一些提示，请参阅高级表达式操作部分。

expand_func

要根据一些恒等式扩展特殊函数，请使用expand_func()。例如

>>> expand_func(gamma(x + 3))
x⋅(x + 1)⋅(x + 2)⋅Γ(x) 

hyperexpand

要用更标准的函数重写hyper，请使用hyperexpand()。

>>> hyperexpand(hyper([1, 1], [2], z))
-log(1 - z)
────────────
 z 

hyperexpand()还适用于更一般的 Meijer G 函数（有关更多信息，请参阅其文档）。

>>> expr = meijerg([[1],[1]], [[1],[]], -z)
>>> expr
╭─╮1, 1 ⎛1  1 │   ⎞
│╶┐     ⎜     │ -z⎟
╰─╯2, 1 ⎝1    │   ⎠
>>> hyperexpand(expr)
 1
 ─
 z
ℯ 

combsimp

要简化组合表达式，使用combsimp()。

>>> n, k = symbols('n k', integer = True)
>>> combsimp(factorial(n)/factorial(n - 3))
n⋅(n - 2)⋅(n - 1)
>>> combsimp(binomial(n+1, k+1)/binomial(n, k))
n + 1
─────
k + 1 

gammasimp

要简化带有伽玛函数或非整数参数的组合函数的表达式，请使用gammasimp()。

>>> gammasimp(gamma(x)*gamma(1 - x))
 π
────────
sin(π⋅x) 

示例：连分数

让我们使用 SymPy 探索连分数。连分数是形式为

[a_0 + \cfrac{1}{a_1 + \cfrac{1}{a_2 + \cfrac{1}{ \ddots + \cfrac{1}{a_n} }}}]

其中 (a_0, \ldots, a_n) 是整数，(a_1, \ldots, a_n) 是正数。连分数也可以是无限的，但无限对象在计算机中更难表示，因此我们这里只讨论有限情况。

上述形式的连分数通常表示为一个列表 ([a_0; a_1, \ldots, a_n])。让我们写一个简单的函数，将这样的列表转换为其连分数形式。从列表构造连分数的最简单方法是从后向前工作。请注意，尽管定义表面上对称，但第一个元素 (a_0) 通常必须与其他元素不同地处理。

>>> def list_to_frac(l):
...     expr = Integer(0)
...     for i in reversed(l[1:]):
...         expr += i
...         expr = 1/expr
...     return l[0] + expr
>>> list_to_frac([x, y, z])
 1
x + ─────
 1
 y + ─
 z 

我们在 list_to_frac 中使用 Integer(0)，这样即使我们只传入 Python 整数，结果也将始终是一个 SymPy 对象。

>>> list_to_frac([1, 2, 3, 4])
43
──
30 

每个有限的连分数都是有理数，但在这里我们对符号感兴趣，因此让我们创建一个符号连分数。我们一直在使用的 symbols() 函数有一个快捷方式来创建带编号的符号。symbols('a0:5') 将创建符号 a0、a1，直到 a4。

>>> syms = symbols('a0:5')
>>> syms
(a₀, a₁, a₂, a₃, a₄)
>>> a0, a1, a2, a3, a4 = syms
>>> frac = list_to_frac(syms)
>>> frac
 1
a₀ + ─────────────────
 1
 a₁ + ────────────
 1
 a₂ + ───────
 1
 a₃ + ──
 a₄ 

这种形式对于理解连分数很有用，但让我们用 cancel() 将其放入标准有理函数形式中。

>>> frac = cancel(frac)
>>> frac
a₀⋅a₁⋅a₂⋅a₃⋅a₄ + a₀⋅a₁⋅a₂ + a₀⋅a₁⋅a₄ + a₀⋅a₃⋅a₄ + a₀ + a₂⋅a₃⋅a₄ + a₂ + a₄
─────────────────────────────────────────────────────────────────────────
 a₁⋅a₂⋅a₃⋅a₄ + a₁⋅a₂ + a₁⋅a₄ + a₃⋅a₄ + 1 

现在假设我们在上述取消的形式中给出了 frac。事实上，我们可能会以任何形式得到分数，但我们总是可以用 cancel() 将其转换为上述的标准形式。假设我们知道它可以被重写为一个连分数。我们可以如何使用 SymPy 做到这一点？一个连分数递归地是 (c + \frac{1}{f})，其中 (c) 是一个整数，(f) 是一个（较小的）连分数。如果我们能以这种形式写出表达式，我们就可以递归地提取每个 (c) 并添加到一个列表中。然后我们可以用我们的 list_to_frac() 函数得到一个连分数。

这里的关键观察是，我们可以通过对 (c) 进行部分分解来将一个表达式转换为 (c + \frac{1}{f}) 的形式。这是因为 (f) 不包含 (c)。这意味着我们需要使用 apart() 函数。我们使用 apart() 来分离项，然后从表达式中减去它，并取倒数来得到 (f) 部分。

>>> l = []
>>> frac = apart(frac, a0)
>>> frac
 a₂⋅a₃⋅a₄ + a₂ + a₄
a₀ + ───────────────────────────────────────
 a₁⋅a₂⋅a₃⋅a₄ + a₁⋅a₂ + a₁⋅a₄ + a₃⋅a₄ + 1
>>> l.append(a0)
>>> frac = 1/(frac - a0)
>>> frac
a₁⋅a₂⋅a₃⋅a₄ + a₁⋅a₂ + a₁⋅a₄ + a₃⋅a₄ + 1
───────────────────────────────────────
 a₂⋅a₃⋅a₄ + a₂ + a₄ 

现在我们重复这个过程。

>>> frac = apart(frac, a1)
>>> frac
 a₃⋅a₄ + 1
a₁ + ──────────────────
 a₂⋅a₃⋅a₄ + a₂ + a₄
>>> l.append(a1)
>>> frac = 1/(frac - a1)
>>> frac = apart(frac, a2)
>>> frac
 a₄
a₂ + ─────────
 a₃⋅a₄ + 1
>>> l.append(a2)
>>> frac = 1/(frac - a2)
>>> frac = apart(frac, a3)
>>> frac
 1
a₃ + ──
 a₄
>>> l.append(a3)
>>> frac = 1/(frac - a3)
>>> frac = apart(frac, a4)
>>> frac
a₄
>>> l.append(a4)
>>> list_to_frac(l)
 1
a₀ + ─────────────────
 1
 a₁ + ────────────
 1
 a₂ + ───────
 1
 a₃ + ──
 a₄ 

当然，这个练习似乎毫无意义，因为我们已经知道我们的 frac 是 list_to_frac([a0, a1, a2, a3, a4])。所以试试以下练习。取一个符号列表并将它们随机化，然后创建被取消的连分数，并看看能否复制原始列表。例如

>>> import random
>>> l = list(symbols('a0:5'))
>>> random.shuffle(l)
>>> orig_frac = frac = cancel(list_to_frac(l))
>>> del l 

在 SymPy 中，对上面的例子，尝试从 frac 复制出 l。我已经删除了末尾的 l 以消除偷看的诱惑（你可以在最后调用 cancel(list_to_frac(l)) 来检查你的答案，并与 orig_frac 进行比较。

看看你能否想出在每个阶段传递给 apart() 的符号的方法（提示：想想在公式 (a_0 + \frac{1}{a_1 + \cdots}) 中 (a_0) 发生了什么）。

微积分

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

本节介绍了如何在 SymPy 中执行基本的微积分任务，例如导数、积分、极限和级数展开。如果您对本节的任何数学内容不熟悉，可以放心跳过。

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

导数

要进行导数计算，请使用diff函数。

>>> diff(cos(x), x)
-sin(x)
>>> diff(exp(x**2), x)
 ⎛ 2⎞
 ⎝x ⎠
2⋅x⋅ℯ 

diff可以一次进行多个导数的计算。要进行多次导数计算，传递变量多少次即可，或在变量后传递一个数字。例如，以下两个示例都找到了(x⁴)的三阶导数。

>>> diff(x**4, x, x, x)
24⋅x
>>> diff(x**4, x, 3)
24⋅x 

您还可以一次性相对多个变量进行导数计算。只需按顺序传递每个导数，使用与单变量导数相同的语法。例如，以下每个示例都将计算(\frac{\partial⁷}{\partial x\partial y²\partial z⁴} e^{x y z})。

>>> expr = exp(x*y*z)
>>> diff(expr, x, y, y, z, z, z, z)
 3  2 ⎛ 3  3  3       2  2  2                ⎞  x⋅y⋅z
x ⋅y ⋅⎝x ⋅y ⋅z  + 14⋅x ⋅y ⋅z  + 52⋅x⋅y⋅z + 48⎠⋅ℯ
>>> diff(expr, x, y, 2, z, 4)
 3  2 ⎛ 3  3  3       2  2  2                ⎞  x⋅y⋅z
x ⋅y ⋅⎝x ⋅y ⋅z  + 14⋅x ⋅y ⋅z  + 52⋅x⋅y⋅z + 48⎠⋅ℯ
>>> diff(expr, x, y, y, z, 4)
 3  2 ⎛ 3  3  3       2  2  2                ⎞  x⋅y⋅z
x ⋅y ⋅⎝x ⋅y ⋅z  + 14⋅x ⋅y ⋅z  + 52⋅x⋅y⋅z + 48⎠⋅ℯ 

diff也可以作为方法调用。调用diff的两种方式完全相同，仅为方便起见。

>>> expr.diff(x, y, y, z, 4)
 3  2 ⎛ 3  3  3       2  2  2                ⎞  x⋅y⋅z
x ⋅y ⋅⎝x ⋅y ⋅z  + 14⋅x ⋅y ⋅z  + 52⋅x⋅y⋅z + 48⎠⋅ℯ 

要创建一个未计算的导数，请使用Derivative类。它具有与diff相同的语法。

>>> deriv = Derivative(expr, x, y, y, z, 4)
>>> deriv
 7
 ∂     ⎛ x⋅y⋅z⎞
──────────⎝ℯ     ⎠
 4   2
∂z  ∂y  ∂x 

要计算未计算的导数，使用doit方法。

>>> deriv.doit()
 3  2 ⎛ 3  3  3       2  2  2                ⎞  x⋅y⋅z
x ⋅y ⋅⎝x ⋅y ⋅z  + 14⋅x ⋅y ⋅z  + 52⋅x⋅y⋅z + 48⎠⋅ℯ 

这些未计算的对象对于延迟导数的计算或打印目的非常有用。当 SymPy 不知道如何计算表达式的导数时（例如，如果它包含未定义的函数，这些函数在 Solving Differential Equations 部分中描述），它们也会被使用。

可以使用元组(x, n)创建未指定顺序的导数，其中n是相对于x的导数的阶数。

>>> m, n, a, b = symbols('m n a b')
>>> expr = (a*x + b)**m
>>> expr.diff((x, n))
 n
 ∂ ⎛         m⎞
───⎝(a⋅x + b) ⎠
 n
∂x 

积分

要计算积分，请使用integrate函数。有两种积分类型，即定积分和不定积分。要计算不定积分，即反导数或原函数，只需在表达式后传递变量。

>>> integrate(cos(x), x)
sin(x) 

请注意，SymPy 不包括积分常数。如果需要，您可以自行添加，或将问题重新表述为微分方程并使用dsolve来解决，后者会添加常数（请参阅 Solving Differential Equations）。

要计算定积分，请传递参数(integration_variable, lower_limit, upper_limit)。例如，要计算

[\int_0^\infty e^{-x},dx,]

我们将执行

>>> integrate(exp(-x), (x, 0, oo))
1 

与不定积分一样，您可以传递多个限制元组以执行多重积分。例如，要计算

[\int_{-\infty}^{{\infty}\int_{-\infty}} e^{- x^{2} - y^{2}}, dx, dy,]

原文：

>>> integrate(exp(-x**2 - y**2), (x, -oo, oo), (y, -oo, oo))
π 

如果integrate无法计算积分，它会返回一个未计算的Integral对象。

>>> expr = integrate(x**x, x)
>>> print(expr)
Integral(x**x, x)
>>> expr
⌠
⎮  x
⎮ x  dx
⌡ 

与Derivative一样，您可以使用Integral创建一个未计算的积分。要稍后评估此积分，请调用doit。

>>> expr = Integral(log(x)**2, x)
>>> expr
⌠
⎮    2
⎮ log (x) dx
⌡
>>> expr.doit()
 2
x⋅log (x) - 2⋅x⋅log(x) + 2⋅x 

integrate使用强大的算法来计算定积分和不定积分，包括启发式模式匹配类型算法、Risch 算法的部分实现以及使用Meijer G 函数的算法，后者特别适用于以特殊函数形式计算积分，尤其是定积分。以下是integrate的一些强大功能示例。

>>> integ = Integral((x**4 + x**2*exp(x) - x**2 - 2*x*exp(x) - 2*x -
...     exp(x))*exp(x)/((x - 1)**2*(x + 1)**2*(exp(x) + 1)), x)
>>> integ
⌠
⎮ ⎛ 4    2  x    2        x          x⎞  x
⎮ ⎝x  + x ⋅ℯ  - x  - 2⋅x⋅ℯ  - 2⋅x - ℯ ⎠⋅ℯ
⎮ ──────────────────────────────────────── dx
⎮               2        2 ⎛ x    ⎞
⎮        (x - 1) ⋅(x + 1) ⋅⎝ℯ  + 1⎠
⌡
>>> integ.doit()
 x
 ⎛ x    ⎞     ℯ
log⎝ℯ  + 1⎠ + ──────
 2
 x  - 1 

>>> integ = Integral(sin(x**2), x)
>>> integ
⌠
⎮    ⎛ 2⎞
⎮ sin⎝x ⎠ dx
⌡
>>> integ.doit()
 ⎛√2⋅x⎞
3⋅√2⋅√π⋅S⎜────⎟⋅Γ(3/4)
 ⎝ √π ⎠
──────────────────────
 8⋅Γ(7/4) 

>>> integ = Integral(x**y*exp(-x), (x, 0, oo))
>>> integ
∞
⌠
⎮  y  -x
⎮ x ⋅ℯ   dx
⌡
0
>>> integ.doit()
⎧ Γ(y + 1)    for re(y) > -1
⎪
⎪∞
⎪⌠
⎨⎮  y  -x
⎪⎮ x ⋅ℯ   dx    otherwise
⎪⌡
⎪0
⎩ 

这个最后的例子返回了一个Piecewise表达式，因为积分只有在(\Re(y) > -1.)时才收敛。

数值积分

数值积分是数学分析中用来估计函数在简化范围内定积分的方法。SymPy 不仅支持符号积分，还支持数值积分。它利用mpmath库的精度功能来增强数值积分计算的准确性。

>>> from sympy import Integral, Symbol, sqrt
>>> x = Symbol('x')
>>> integral = Integral(sqrt(2)*x, (x, 0, 1))
>>> integral
1
⌠
⎮ √2⋅x dx
⌡
0
>>> integral.evalf()
0.707106781186548 

要计算具有指定精度的积分：

>>> integral.evalf(50)
0.70710678118654752440084436210484903928483593768847 

数值积分在符号积分不可行或不可能的情况下成为一种可行的方法。这种方法允许通过数值技术计算积分，即使处理无限区间或被积函数时也是如此：

>>> Integral(exp(-(x ** 2)), (x, -oo, oo)).evalf()
1.77245385090552 

>>> Integral(1 / sqrt(x), (x, 0, 1)).evalf()
2.00000000000000 

极限

SymPy 可以使用limit函数计算符号极限。计算

[\lim_{x\to x_0} f(x)]

是limit(f(x), x, x0)。

>>> limit(sin(x)/x, x, 0)
1 

当评估点是奇点时，应该使用limit而不是subs。尽管 SymPy 有表示(\infty)的对象，但在评估时不可靠，因为它们不会跟踪增长速度等信息。此外，诸如(\infty - \infty)和(\frac{\infty}{\infty})会返回(\mathrm{nan})（非数字）。例如

>>> expr = x**2/exp(x)
>>> expr.subs(x, oo)
nan
>>> limit(expr, x, oo)
0 

与Derivative和Integral类似，limit也有一个未求值的对应物Limit。要对其求值，请使用doit。

>>> expr = Limit((cos(x) - 1)/x, x, 0)
>>> expr
 ⎛cos(x) - 1⎞
 lim ⎜──────────⎟
x─→0⁺⎝    x     ⎠
>>> expr.doit()
0 

要在单侧计算极限，将'+'或'-'作为第四个参数传递给limit函数。例如，要计算

[\lim_{x\to 0^+}\frac{1}{x},]

做

>>> limit(1/x, x, 0, '+')
∞ 

与

>>> limit(1/x, x, 0, '-')
-∞ 

级数展开

SymPy 可以计算函数在某点周围的渐近级数展开。要计算(f(x))在点(x = x_0)处的(x^n)阶展开，可以使用f(x).series(x, x0, n)。可以省略x0和n，此时将使用默认值x0=0和n=6。

>>> expr = exp(sin(x))
>>> expr.series(x, 0, 4)
 2
 x     ⎛ 4⎞
1 + x + ── + O⎝x ⎠
 2 

结尾处的(O\left(x⁴\right))项表示在(x=0)处的 Landau 阶次项（不要与计算机科学中使用的大 O 符号混淆，后者通常表示(x \rightarrow \infty)时的 Landau 阶次项）。它意味着所有幂次大于或等于(x⁴)的 x 项都被省略了。阶次项可以在series之外创建和操作。它们会自动吸收更高阶次的项。

>>> x + x**3 + x**6 + O(x**4)
 3    ⎛ 4⎞
x + x  + O⎝x ⎠
>>> x*O(1)
O(x) 

如果不想要顺序项，请使用removeO方法。

>>> expr.series(x, 0, 4).removeO()
 2
x
── + x + 1
2 

O符号支持任意的极限点（非 0）：

>>> exp(x - 6).series(x, x0=6)
 2          3          4          5
 (x - 6)    (x - 6)    (x - 6)    (x - 6)         ⎛       6       ⎞
-5 + ──────── + ──────── + ──────── + ──────── + x + O⎝(x - 6) ; x → 6⎠
 2          6          24        120 

有限差分

到目前为止，我们分别查看了具有解析导数和原始函数的表达式。但是，如果我们想要估计缺乏闭合形式表示的曲线的导数，或者我们尚不知道其功能值，该怎么办呢？一种方法是使用有限差分方法。

使用有限差分最简单的方法是使用differentiate_finite函数：

>>> f, g = symbols('f g', cls=Function)
>>> differentiate_finite(f(x)*g(x))
-f(x - 1/2)⋅g(x - 1/2) + f(x + 1/2)⋅g(x + 1/2) 

如果您已经有一个Derivative实例，可以使用as_finite_difference方法生成任意阶导数的近似值：

>>> f = Function('f')
>>> dfdx = f(x).diff(x)
>>> dfdx.as_finite_difference()
-f(x - 1/2) + f(x + 1/2) 

这里对 x 周围的一阶导数使用了最少数量的点（一阶导数为 2 个点），等间距地使用步长为 1 进行评估。我们可以使用任意步长（可能包含符号表达式）：

>>> f = Function('f')
>>> d2fdx2 = f(x).diff(x, 2)
>>> h = Symbol('h')
>>> d2fdx2.as_finite_difference([-3*h,-h,2*h])
f(-3⋅h)   f(-h)   2⋅f(2⋅h)
─────── - ───── + ────────
 2        2        2
 5⋅h      3⋅h     15⋅h 

如果您只对评估权重感兴趣，可以手动执行：

>>> finite_diff_weights(2, [-3, -1, 2], 0)[-1][-1]
[1/5, -1/3, 2/15] 

注意我们只需要从finite_diff_weights返回的最后一个子列表的最后一个元素。原因是该函数还为更低阶导数生成权重，并且使用更少的点（详见finite_diff_weights的文档以获取更多详情）。

如果直接使用finite_diff_weights看起来复杂，并且Derivative实例的as_finite_difference方法不够灵活，你可以使用apply_finite_diff，它接受order、x_list、y_list和x0作为参数：

>>> x_list = [-3, 1, 2]
>>> y_list = symbols('a b c')
>>> apply_finite_diff(1, x_list, y_list, 0)
 3⋅a   b   2⋅c
- ─── - ─ + ───
 20   4    5 

解算器

原文：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}})离子释放。这些离子扩散并在肌肉中增加浓度，导致激活。传输到肌肉以刺激收缩的电信号称为“兴奋”。在生物力学模型中，通常用符号(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}{lT_{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)和(lT_{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}{lM_{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)和(lM_{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]} = \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}{vM_{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] = \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}{lM_{opt}}) 时，肌肉纤维速度可以表达为

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

如前所述，肌肉纤维可以归一化，(\tilde{v}^M = \frac{v^M}{vM_{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}{lT_{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}{lM} \ 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}}{lM_{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

生物力学模型示例

原文：docs.sympy.org/latest/tutorials/physics/biomechanics/biomechanical-model-example.html

sympy.physics.biomechanics 提供了增强使用 sympy.physics.mechanics 创建的模型的功能，包括模拟肌肉和肌腱的力产生元素。在本教程中，我们将通过向简单的人体手臂模型添加肌肉来介绍该包的功能。

模型描述

示意图显示了杠杆 (A)、上臂 (C) 和下臂 (D)。

杠杆 (A) 可以围绕 (\hat{n}_z) 轴旋转角度 (q_1)。其质心位于旋转轴上。肩位于 (P_2)，上臂 (C) 可以围绕 (\hat{n}_y) 轴伸展角度 (q_2) 并围绕 (\hat{b}_z) 轴旋转角度 (q_3)。肘位于点 (P_3)。下臂可以围绕 (\hat{c}_y) 轴弯曲角度 (q_4)。手位于点 (P_4)。手将通过强制条件 (\mathbf{r}^{P_4/O} = \mathbf{r}^{P_1/O}) 限制在杠杆上。杠杆、上臂和下臂将被建模为惯性简单的细圆柱体。

我们将引入两种肌肉肌腱模型，分别代表肱二头肌和肱三头肌。两个肌肉附着点 (C_m) 和 (D_m) 固定在上臂和下臂上。肱二头肌在从 (C_m) 到 (D_m) 的线性路径上发挥作用，收缩时导致肘部弯曲。定义一个半径为 (r) 的圆弧，其中心在 (P_3) 处，法线方向为 (\hat{c}_y)。肱三头肌将围绕圆弧绕过，并与肱二头肌相同的点连接，收缩时导致肘部伸展。

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

定义变量

引入四个坐标 (\mathbf{q} = [q_1, q_2, q_3, q_4]^T)，分别表示杠杆角度、肩部伸展、肩部旋转和肘部弯曲。我们还需要广义速度 (\mathbf{u} = [u_1,u_2,u_3,u_4]^T)，定义为 (\mathbf{u} = \dot{\mathbf{q}})。

>>> q1, q2, q3, q4 = me.dynamicsymbols('q1, q2, q3, q4', real=True)
>>> u1, u2, u3, u4 = me.dynamicsymbols('u1, u2, u3, u4', real=True) 

机械系统所需的常数参数包括：

• (d_x, l_A)：从 (O) 沿 (\hat{n}_x) 和 (\hat{a}_y) 方向定位 (P_1) 的位置

• (d_y, d_z)：从 (O) 沿 (N) 单位向量方向定位 (P_2) 的位置

• (l_C,l_D)：上臂和下臂的长度

• (m_A,m_C,m_D)：杠杆、上臂和下臂的质量

• (g)：重力加速度

• (k)：杠杆的线性转动弹簧系数

• (c)：杠杆的线性转动阻尼系数

>>> dx, dy, dz = sm.symbols('dx, dy, dz', real=True, nonnegative=True)
>>> lA, lC, lD = sm.symbols('lA, lC, lD', real=True, positive=True)
>>> mA, mC, mD = sm.symbols('mA, mC, mD', real=True, positive=True)
>>> g, k, c, r = sm.symbols('g, k, c, r', real=True, positive=True) 

定义运动学

定义 Schematic showing the lever A and the upper C and lower D arm.中显示的所有参考框架和点。(C_o)和(D_o)分别是上臂和下臂的质心。

>>> N, A, B, C, D = sm.symbols('N, A, B, C, D', cls=me.ReferenceFrame)
>>> O, P1, P2, P3, P4 = sm.symbols('O, P1, P2, P3, P4 ', cls=me.Point)
>>> Ao, Co, Cm, Dm, Do = sm.symbols('Ao, Co, Cm, Dm, Do', cls=me.Point) 

参考框架的方向和角速度如下：

>>> A.orient_axis(N, q1, N.z)
>>> B.orient_axis(N, q2, N.y)
>>> C.orient_axis(B, q3, B.z)
>>> D.orient_axis(C, q4, C.y)
>>> A.set_ang_vel(N, u1*N.z)
>>> B.set_ang_vel(N, u2*N.y)
>>> C.set_ang_vel(B, u3*B.z)
>>> D.set_ang_vel(C, u4*C.y) 

所有点的位置和速度如下：

>>> Ao.set_pos(O, dx*N.x)
>>> P1.set_pos(Ao, lA*A.y)
>>> P2.set_pos(O, dy*N.y + dz*N.z)
>>> Co.set_pos(P2, lC/2*C.z)
>>> Cm.set_pos(P2, 1*lC/3*C.z)
>>> P3.set_pos(P2, lC*C.z)
>>> Dm.set_pos(P3, 1*lD/3*D.z)
>>> Do.set_pos(P3, lD/2*D.z)
>>> P4.set_pos(P3, lD*D.z) 

>>> O.set_vel(N, 0)
>>> Ao.set_vel(N, 0)
>>> P1.v2pt_theory(Ao, N, A)
- lA*u1(t)*A.x
>>> P2.set_vel(N, 0)
>>> Co.v2pt_theory(P2, N, C)
lC*u2(t)*cos(q3(t))/2*C.x - lC*u2(t)*sin(q3(t))/2*C.y
>>> Cm.v2pt_theory(P2, N, C)
lC*u2(t)*cos(q3(t))/3*C.x - lC*u2(t)*sin(q3(t))/3*C.y
>>> P3.v2pt_theory(P2, N, C)
lC*u2(t)*cos(q3(t))*C.x - lC*u2(t)*sin(q3(t))*C.y
>>> Dm.v2pt_theory(P3, N, D)
lC*u2(t)*cos(q3(t))*C.x - lC*u2(t)*sin(q3(t))*C.y + lD*(u2(t)*cos(q3(t)) + u4(t))/3*D.x - lD*(u2(t)*sin(q3(t))*cos(q4(t)) - u3(t)*sin(q4(t)))/3*D.y
>>> Do.v2pt_theory(P3, N, D)
lC*u2(t)*cos(q3(t))*C.x - lC*u2(t)*sin(q3(t))*C.y + lD*(u2(t)*cos(q3(t)) + u4(t))/2*D.x - lD*(u2(t)*sin(q3(t))*cos(q4(t)) - u3(t)*sin(q4(t)))/2*D.y
>>> P4.v2pt_theory(P3, N, D)
lC*u2(t)*cos(q3(t))*C.x - lC*u2(t)*sin(q3(t))*C.y + lD*(u2(t)*cos(q3(t)) + u4(t))*D.x - lD*(u2(t)*sin(q3(t))*cos(q4(t)) - u3(t)*sin(q4(t)))*D.y 

有三个保角约束方程式需要保持手(P_4)在杠杆(P_1)上：

>>> holonomic = (P4.pos_from(O) - P1.pos_from(O)).to_matrix(N) 

定义惯性

惯性二阶张量然后可以假设杠杆、上臂和下臂是细长圆柱体：

>>> IA = me.Inertia(me.inertia(A, mA/12*lA**2, mA/2*lA**2, mA/12*lA**2), Ao)
>>> IC = me.Inertia(me.inertia(C, mC/12*lC**2, mC/12*lC**2, mC/2*lC**2), Co)
>>> ID = me.Inertia(me.inertia(D, mD/12*lD**2, mD/12*lD**2, mD/2*lD**2), Do) 

>>> lever = me.RigidBody('lever', masscenter=Ao, frame=A, mass=mA, inertia=IA)
>>> u_arm = me.RigidBody('upper arm', masscenter=Co, frame=C, mass=mC, inertia=IC)
>>> l_arm = me.RigidBody('lower arm', masscenter=Do, frame=D, mass=mD, inertia=ID) 

定义力

我们将在地球的重力场中模拟这个系统：

>>> gravC = me.Force(u_arm, mC*g*N.z)
>>> gravD = me.Force(l_arm, mD*g*N.z) 

杠杆具有惯性，但我们还将添加一个线性扭转弹簧和阻尼器，以提供更多的抵抗力供手臂推压和拉动：

>>> lever_resistance = me.Torque(A, (-k*q1 - c*u1)*N.z) 

肱二头肌

我们将模拟肱二头肌作为一个在两个肌肉附着点(C_m)和(D_m)之间收缩的致动器。给定指定的激励输入，这种肌肉可以收缩，而且我们将假设相关的腱是刚性的。肌肉腱致动器模型由两个组成部分组成：一个用于作用的路径和定义如何传播激励输入以激活肌肉的激活动力学。肱二头肌将沿着LinearPathway作用，并将使用从[DeGroote2016]中推导的特定肌肉动力学实现。

首先创建线性路径：

>>> biceps_pathway = me.LinearPathway(Cm, Dm) 

您可以创建一个激活模型，该模型完全符号化，或者使用从[DeGroote2016]中获取的特定调整的数值参数创建它（推荐）：

>>> biceps_activation = bm.FirstOrderActivationDeGroote2016.with_defaults('biceps') 

然后，将完整的肌肉腱致动器模型命名并构建与匹配的类：

>>> biceps = bm.MusculotendonDeGroote2016.with_defaults('biceps', biceps_pathway, biceps_activation) 

肱三头肌

肱三头肌致动器模型将需要一个定制的路径来管理肌肉和腱在半径为(r)的圆弧周围包裹的性质。这个路径由两个不改变长度的线性段和一个圆弧组成，随着肘部的伸展和屈曲而改变长度。作用于上臂和下臂的力可以被建模为分别作用于点(C_m)和(D_m)上的力，这些力始终与线性段平行，并且在圆弧端点上作用的相等且相反的力在(P_3)处产生一个合力。

为了开发这条路径，我们需要子类化 PathwayBase 并创建计算路径长度、路径扩展速度以及作用于涉及身体的负载的方法。我们将开发一个类，假设两个刚体之间有个销钉连接，两个肌肉附着点分别固定在每个身体上，并且销钉连接点和两个附着点位于与销钉轴线垂直的同一平面上。我们还假设销钉连接坐标以 (q_4) 为测量，其在 显示杠杆 A 及上 C 和下 D 臂的示意图。 中，并且 (0 \le q_4 \le \pi)。圆弧具有半径 (r)。在这些假设条件下，我们可以使用 __init__() 方法收集剩余方法中所需的必要信息。在 __init__() 中，我们还将计算一些在多个重载方法中需要的量。路径的长度是两条线性段和随销钉连接坐标变化的圆弧长度的总和。扩展速度简单地是弧长随时间的变化。负载由三个力组成：两个在路径线性部分的起始点和插入点上推拉的力以及对肘部的力量推拉在圆弧端点上的结果。

>>> class ExtensorPathway(me.PathwayBase):
...
...     def __init__(self, origin, insertion, axis_point, axis, parent_axis,
...                  child_axis, radius, coordinate):
...  """A custom pathway that wraps a circular arc around a pin joint.
...
...         This is intended to be used for extensor muscles. For example, a
...         triceps wrapping around the elbow joint to extend the upper arm at
...         the elbow.
...
...         Parameters
...         ==========
...         origin : Point
...             Muscle origin point fixed on the parent body (A).
...         insertion : Point
...             Muscle insertion point fixed on the child body (B).
...         axis_point : Point
...             Pin joint location fixed in both the parent and child.
...         axis : Vector
...             Pin joint rotation axis.
...         parent_axis : Vector
...             Axis fixed in the parent frame (A) that is directed from the pin
...             joint point to the muscle origin point.
...         child_axis : Vector
...             Axis fixed in the child frame (B) that is directed from the pin
...             joint point to the muscle insertion point.
...         radius : sympyfiable
...             Radius of the arc that the muscle wraps around.
...         coordinate : sympfiable function of time
...             Joint angle, zero when parent and child frames align. Positive
...             rotation about the pin joint axis, B with respect to A.
...
...         Notes
...         =====
...
...         Only valid for coordinate >= 0.
...
...         """
...         super().__init__(origin, insertion)
...
...         self.origin = origin
...         self.insertion = insertion
...         self.axis_point = axis_point
...         self.axis = axis.normalize()
...         self.parent_axis = parent_axis.normalize()
...         self.child_axis = child_axis.normalize()
...         self.radius = radius
...         self.coordinate = coordinate
...
...         self.origin_distance = axis_point.pos_from(origin).magnitude()
...         self.insertion_distance = axis_point.pos_from(insertion).magnitude()
...         self.origin_angle = sm.asin(self.radius/self.origin_distance)
...         self.insertion_angle = sm.asin(self.radius/self.insertion_distance)
...
...     @property
...     def length(self):
...  """Length of the pathway.
...
...         Length of two fixed length line segments and a changing arc length
...         of a circle.
...
...         """
...
...         angle = self.origin_angle + self.coordinate + self.insertion_angle
...         arc_length = self.radius*angle
...
...         origin_segment_length = self.origin_distance*sm.cos(self.origin_angle)
...         insertion_segment_length = self.insertion_distance*sm.cos(self.insertion_angle)
...
...         return origin_segment_length + arc_length + insertion_segment_length
...
...     @property
...     def extension_velocity(self):
...  """Extension velocity of the pathway.
...
...         Arc length of circle is the only thing that changes when the elbow
...         flexes and extends.
...
...         """
...         return self.radius*self.coordinate.diff(me.dynamicsymbols._t)
...
...     def to_loads(self, force_magnitude):
...  """Loads in the correct format to be supplied to `KanesMethod`.
...
...         Forces applied to origin, insertion, and P from the muscle wrapped
...         over circular arc of radius r.
...
...         """
...
...         parent_tangency_point = me.Point('Aw')  # fixed in parent
...         child_tangency_point = me.Point('Bw')  # fixed in child
...
...         parent_tangency_point.set_pos(
...             self.axis_point,
...             -self.radius*sm.cos(self.origin_angle)*self.parent_axis.cross(self.axis)
...             + self.radius*sm.sin(self.origin_angle)*self.parent_axis,
...         )
...         child_tangency_point.set_pos(
...             self.axis_point,
...             self.radius*sm.cos(self.insertion_angle)*self.child_axis.cross(self.axis)
...             + self.radius*sm.sin(self.insertion_angle)*self.child_axis),
...
...         parent_force_direction_vector = self.origin.pos_from(parent_tangency_point)
...         child_force_direction_vector = self.insertion.pos_from(child_tangency_point)
...         force_on_parent = force_magnitude*parent_force_direction_vector.normalize()
...         force_on_child = force_magnitude*child_force_direction_vector.normalize()
...         loads = [
...             me.Force(self.origin, force_on_parent),
...             me.Force(self.axis_point, -(force_on_parent + force_on_child)),
...             me.Force(self.insertion, force_on_child),
...         ]
...         return loads
... 

现在我们已经定义了一个自定义路径，可以以与肱二头肌相同的方式创建肌腱执行器模型：

>>> triceps_pathway = ExtensorPathway(Cm, Dm, P3, B.y, -C.z, D.z, r, q4)
>>> triceps_activation = bm.FirstOrderActivationDeGroote2016.with_defaults('triceps')
>>> triceps = bm.MusculotendonDeGroote2016.with_defaults('triceps', triceps_pathway, triceps_activation) 

最后，所有负载可以汇总到一个列表中：

>>> loads = biceps.to_loads() + triceps.to_loads() + [lever_resistance, gravC, gravD] 

运动方程

现在定义了所有负载，系统的运动方程可以生成。我们有三个完整约束，因此系统仅有一个自由度。

>>> kane = me.KanesMethod(
...     N,
...     (q1,),
...     (u1,),
...     kd_eqs=(
...         u1 - q1.diff(),
...         u2 - q2.diff(),
...         u3 - q3.diff(),
...         u4 - q4.diff(),
...     ),
...     q_dependent=(q2, q3, q4),
...     configuration_constraints=holonomic,
...     velocity_constraints=holonomic.diff(me.dynamicsymbols._t),
...     u_dependent=(u2, u3, u4),
... )
...
>>> Fr, Frs = kane.kanes_equations((lever, u_arm, l_arm), loads) 

不在 (\dot{\mathbf{u}}) 中的项包含与肌肉力相关的函数，这些函数除了坐标和广义速度外，还是激活状态变量的函数。

>>> me.find_dynamicsymbols(kane.forcing)
{a_biceps(t), a_triceps(t), q1(t), q2(t), q3(t), q4(t), u1(t), u2(t), u3(t), u4(t)} 

它们还包含与肌肉模型相关的新常数参数：

>>> kane.forcing.free_symbols
{F_M_max_biceps, F_M_max_triceps, c, g, k, lA, lC, lD, l_M_opt_biceps, l_M_opt_triceps, l_T_slack_biceps, l_T_slack_triceps, mC, mD, r, t} 

肌肉激活微分方程

每个肌肉的激活状态是与两个新的一阶微分方程相关联的新状态变量。这些微分方程可以从肌肉执行器模型中获取：

>>> biceps.rhs()
Matrix([[((1/2 - tanh(10.0*a_biceps(t) - 10.0*e_biceps(t))/2)/(0.0225*a_biceps(t) + 0.0075) + 16.6666666666667*(3*a_biceps(t)/2 + 1/2)*(tanh(10.0*a_biceps(t) - 10.0*e_biceps(t))/2 + 1/2))*(-a_biceps(t) + e_biceps(t))]]) 

>>> triceps.rhs()
Matrix([[((1/2 - tanh(10.0*a_triceps(t) - 10.0*e_triceps(t))/2)/(0.0225*a_triceps(t) + 0.0075) + 16.6666666666667*(3*a_triceps(t)/2 + 1/2)*(tanh(10.0*a_triceps(t) - 10.0*e_triceps(t))/2 + 1/2))*(-a_triceps(t) + e_triceps(t))]]) 

将所有肌肉激活微分方程存储在一个矩阵中：

>>> dadt = biceps.rhs().col_join(triceps.rhs()) 

评估系统微分方程

此系统的完整微分方程集的形式如下：

[\begin{split}\begin{bmatrix} \mathbf{I} & \mathbf{0} & \mathbf{0} \ \mathbf{0} & \mathbf{M}_d & \mathbf{0} \ \mathbf{0} & \mathbf{0} & \mathbf{I} \end{bmatrix} \begin{bmatrix} \dot{\mathbf{q}} \ \dot{\mathbf{u}} \ \dot{\mathbf{a}} \end{bmatrix} = \begin{bmatrix} \mathbf{u} \ \mathbf{g}_d(\mathbf{q}, \mathbf{u}, \mathbf{a}) \ \mathbf{g}_a(\mathbf{a}, \mathbf{e}) \end{bmatrix}\end{split}]

在这种情况下，只需求解线性系统以将动力学微分方程组转化为明确形式即可。

要评估系统的方程，我们首先需要收集所有状态、输入和常数变量，以便与lambdify一起使用。状态向量由坐标、广义速度和两个肌肉的激活状态组成：(\mathbf{x}=\begin{bmatrix}\mathbf{q}\\mathbf{u}\\mathbf{a}\end{bmatrix})。

>>> q, u = kane.q, kane.u
>>> a = biceps.x.col_join(triceps.x)
>>> x = q.col_join(u).col_join(a)
>>> x
Matrix([
[       q1(t)],
[       q2(t)],
[       q3(t)],
[       q4(t)],
[       u1(t)],
[       u2(t)],
[       u3(t)],
[       u4(t)],
[ a_biceps(t)],
[a_triceps(t)]]) 

指定输入只有两个肌肉的激活：

>>> e = biceps.r.col_join(triceps.r)
>>> e
Matrix([
[ e_biceps(t)],
[e_triceps(t)]]) 

常数由几何形状、质量、局部重力常数、杠杆的刚度和阻尼系数以及肌肉的各种参数组成。

>>> p = sm.Matrix([
...     dx,
...     dy,
...     dz,
...     lA,
...     lC,
...     lD,
...     mA,
...     mC,
...     mD,
...     g,
...     k,
...     c,
...     r,
...     biceps.F_M_max,
...     biceps.l_M_opt,
...     biceps.l_T_slack,
...     triceps.F_M_max,
...     triceps.l_M_opt,
...     triceps.l_T_slack,
... ])
...
>>> p
Matrix([
[               dx],
[               dy],
[               dz],
[               lA],
[               lC],
[               lD],
[               mA],
[               mC],
[               mD],
[                g],
[                k],
[                c],
[                r],
[   F_M_max_biceps],
[   l_M_opt_biceps],
[ l_T_slack_biceps],
[  F_M_max_triceps],
[  l_M_opt_triceps],
[l_T_slack_triceps]]) 

现在我们有所有符号组件来生成数值函数，以评估(\mathbf{M}_d,\mathbf{g}_d)和(\mathbf{g}_a)的时间导数。通过这些，我们可以计算状态的时间导数。我们还需要一个数值函数来评估全约束条件，以确保配置处于有效状态。

>>> eval_diffeq = sm.lambdify((q, u, a, e, p),
...                           (kane.mass_matrix, kane.forcing, dadt), cse=True)
>>> eval_holonomic = sm.lambdify((q, p), holonomic, cse=True) 

我们需要一些合理的数值常数值：

>>> import numpy as np 

>>> p_vals = np.array([
...     0.31,  # dx [m]
...     0.15,  # dy [m]
...     -0.31,  # dz [m]
...     0.2,   # lA [m]
...     0.3,  # lC [m]
...     0.3,  # lD [m]
...     1.0,  # mA [kg]
...     2.3,  # mC [kg]
...     1.7,  # mD [kg]
...     9.81,  # g [m/s/s]
...     5.0,  # k [Nm/rad]
...     0.5,  # c [Nms/rad]
...     0.03,  # r [m]
...     500.0,  # biceps F_M_max [?]
...     0.6*0.3,  # biceps l_M_opt [?]
...     0.55*0.3,  # biceps l_T_slack [?]
...     500.0,  # triceps F_M_max [?]
...     0.6*0.3,  # triceps l_M_opt [?]
...     0.65*0.3,  # triceps l_T_slack [?]
... ])
... 

由于三个全约束条件，其中三个坐标是剩余一个的函数。我们可以选择杠杆角度(q_1)作为独立坐标，并给出其余坐标的解，假设它们的值已知。

>>> from scipy.optimize import fsolve 

>>> q_vals = np.array([
...     np.deg2rad(5.0),  # q1 [rad]
...     np.deg2rad(-10.0),  # q2 [rad]
...     np.deg2rad(0.0),  # q3 [rad]
...     np.deg2rad(75.0),  # q4 [rad]
... ])
... 

>>> def eval_holo_fsolve(x):
...     q1 = q_vals[0]  # specified
...     q2, q3, q4 = x
...     return eval_holonomic((q1, q2, q3, q4), p_vals).squeeze()
... 

>>> q_vals[1:] = fsolve(eval_holo_fsolve, q_vals[1:]) 

>>> np.rad2deg(q_vals)
[ 5\.         -0.60986636  9.44918589 88.68812842] 

假设系统处于初始静止状态：

>>> u_vals = np.array([
...     0.0,  # u1, [rad/s]
...     0.0,  # u2, [rad/s]
...     0.0,  # u3, [rad/s]
...     0.0,  # u4, [rad/s]
... ])
... 

>>> a_vals = np.array([
...     0.0,  # a_bicep, nondimensional
...     0.0,  # a_tricep, nondimensional
... ]) 

肌肉兴奋刺激也将最初被停用：

>>> e_vals = np.array([
...     0.0,
...     0.0,
... ]) 

现在可以对系统方程进行数值评估：

>>> eval_diffeq(q_vals, u_vals, a_vals, e_vals, p_vals)
([[ 0.00333333 -0.15174161 -0.00109772 -0.00152436]
 [ 0.19923894  0.31       -0.04923615  0.00996712]
 [ 0.01743115  0\.          0.29585191  0.0011276 ]
 [ 0\.         -0.29256885 -0.0005241  -0.29983226]], [[-0.9121071]
 [ 0\.       ]
 [-0\.       ]
 [ 0\.       ]], [[0.]
 [0.]]) 

模拟肌肉驱动的运动

现在，系统方程可以在给定状态和常数值的情况下进行评估，我们可以模拟带有两个肌肉激活的手臂和杠杆运动。如果提供一个函数以明确形式评估它们，例如(\dot{\mathbf{x}}=)，SciPy 的solve_ivp()可以积分这些微分方程。我们将包含一个激活肌肉的函数，但对于第一次模拟设置为零。

>>> def eval_r(t):
...  """Returns the muscles' excitation as a function of time."""
...     e = np.array([0.0, 0.0])
...     return e
...
>>> def eval_rhs(t, x, r, p):
...  """Returns the time derivative of the state.
...
...     Parameters
...     ==========
...     t : float
...         Time in seconds.
...     x : array_like, shape(10,)
...         State vector.
...     r : function
...         Function f(t) that evaluates e.
...     p : array_like, shape(?, )
...         Parameter vector.
...
...     Returns
...     =======
...     dxdt : ndarray, shape(10,)
...       Time derivative of the state.
...
...     """
...
...     q = x[0:4]
...     u = x[4:8]
...     a = x[8:10]
...
...     e = r(t)
...
...     qd = u
...     m, f, ad = eval_diffeq(q, u, a, e, p)
...     ud = np.linalg.solve(m, f).squeeze()
...
...     return np.hstack((qd, ud, ad.squeeze()))
... 

现在，我们可以模拟系统在 3 秒内的运动，提供初始状态(\mathbf{x}_0)和上述函数，使用 SciPy 的solve_ivp()。

>>> from scipy.integrate import solve_ivp 

>>> t0, tf = 0.0, 3.0
>>> ts = np.linspace(t0, tf, num=301)
>>> x0 = np.hstack((q_vals, u_vals, a_vals))
>>> sol = solve_ivp(lambda t, x: eval_rhs(t, x, eval_r, p_vals),
...                 (t0, tf), x0, t_eval=ts) 

可以通过绘制随时间变化的状态轨迹来可视化运动：

>>> import matplotlib.pyplot as plt 

>>> def plot_traj(t, x, syms):
...  """Simple plot of state trajectories.
...
...     Parameters
...     ==========
...     t : array_like, shape(n,)
...         Time values.
...     x : array_like, shape(n, m)
...         State values at each time value.
...     syms : sequence of Symbol, len(m)
...         SymPy symbols associated with state.
...
...     """
...
...     fig, axes = plt.subplots(5, 2, sharex=True)
...
...     for ax, traj, sym in zip(axes.T.flatten(), x.T, syms):
...         if not sym.name.startswith('a'):
...             traj = np.rad2deg(traj)
...         ax.plot(t, traj)
...         ax.set_ylabel(sm.latex(sym, mode='inline'))
...
...     for ax in axes[-1, :]:
...         ax.set_xlabel('Time [s]')
...
...     fig.tight_layout()
...
...     return axes
... 

>>> plot_traj(ts, sol.y.T, x)
[[<Axes: ylabel='$q_{1}{\\left(t \\right)}$'>
 <Axes: ylabel='$u_{2}{\\left(t \\right)}$'>]
 [<Axes: ylabel='$q_{2}{\\left(t \\right)}$'>
 <Axes: ylabel='$u_{3}{\\left(t \\right)}$'>]
 [<Axes: ylabel='$q_{3}{\\left(t \\right)}$'>
 <Axes: ylabel='$u_{4}{\\left(t \\right)}$'>]
 [<Axes: ylabel='$q_{4}{\\left(t \\right)}$'>
 <Axes: ylabel='$a_{biceps}{\\left(t \\right)}$'>]
 [<Axes: xlabel='Time [s]', ylabel='$u_{1}{\\left(t \\right)}$'>
 <Axes: xlabel='Time [s]', ylabel='$a_{triceps}{\\left(t \\right)}$'>]] 

（png, hires.png, pdf)

在没有肌肉激活的情况下，手臂初始杠杆角度为 5 度，达到平衡位置的模拟。

模拟显示，手臂通过重力、杠杆阻力和肌腱模型的被动部分达到平衡。现在我们激活肱二头肌，持续 1 秒，激活率为 80%，以观察其对运动的影响：

>>> def eval_r(t):
...     if t < 0.5 or t > 1.5:
...         e = np.array([0.0, 0.0])
...     else:
...         e = np.array([0.8, 0.0])
...     return e
... 

>>> sol = solve_ivp(lambda t, x: eval_rhs(t, x, eval_r, p_vals), (t0, tf), x0, t_eval=ts) 

>>> plot_traj(ts, sol.y.T, x)
[[<Axes: ylabel='$q_{1}{\\left(t \\right)}$'>
 <Axes: ylabel='$u_{2}{\\left(t \\right)}$'>]
 [<Axes: ylabel='$q_{2}{\\left(t \\right)}$'>
 <Axes: ylabel='$u_{3}{\\left(t \\right)}$'>]
 [<Axes: ylabel='$q_{3}{\\left(t \\right)}$'>
 <Axes: ylabel='$u_{4}{\\left(t \\right)}$'>]
 [<Axes: ylabel='$q_{4}{\\left(t \\right)}$'>
 <Axes: ylabel='$a_{biceps}{\\left(t \\right)}$'>]
 [<Axes: xlabel='Time [s]', ylabel='$u_{1}{\\left(t \\right)}$'>
 <Axes: xlabel='Time [s]', ylabel='$a_{triceps}{\\left(t \\right)}$'>]] 

(png, hires.png, pdf)

在肱二头肌持续激活 1 秒的情况下，手臂初始杠杆角度为 5 度的模拟。

我们首先看到，手臂试图像之前那样达到平衡，但随后激活的肱二头肌将杠杆拉回肩部，导致手臂与被动运动相抗衡。一旦肌肉被停止激活，手臂便会如之前那样平衡。

结论

在这里，我们展示了如何通过构建简单和定制的肌肉-肌腱作用路径创建代表肌肉骨骼系统的数学模型。模型的运动可以通过激活肌肉来控制，模拟显示了预期的行为。

参考文献

[DeGroote2016] (1,2)

De Groote, F., Kinney, A. L., Rao, A. V., & Fregly, B. J.，《评估用于解决肌肉冗余问题的直接协同最优控制问题形式》，生物医学工程学年刊，44(10)，（2016）第 2922-2936 页