企业官网建设流程全解析

深入 SymPy 级数展开模块：Gruntz 极限算法、series 展开、Order 阶项与留数计算

【免费下载链接】sympyA computer algebra system written in pure Python项目地址: https://gitcode.com/GitHub_Trending/sy/sympy

SymPy 的sympy.series模块是纯 Python 计算机代数系统（sympy/）中负责极限计算与级数展开的核心模块。本指南以官方文档 doc/src/modules/series/series.rst 为主线，完整讲解该模块的六大功能：limit()/Limit极限计算、Gruntz 算法及其理论、series()直观级数展开、Order大 O 阶项自动追踪、Richardson/Shanks 级数加速以及residue()留数计算。读完本文，你既能熟练使用这些 API 求解具体问题，也能理解其底层的增长阶比较（MRV 集）与渐近展开原理，并知道如何在源码 sympy/series/ 中定位实现与调试。

模块总览：以极限计算为核心

正如文档开头所述，本模块的主要目的是极限的计算。模块由sympy.series包内多个文件分工实现（sympy/series/）：

这一架构体现了"直观接口 + 底层严格算法"的设计：日常求极限走limit()，它会先用启发式方法快速处理简单情形，复杂情形才落到gruntz()。

极限计算：limit()与Limit类

函数入口limit(e, z, z0, dir="+")

文档给出的核心接口是 limits.py 中的limit()：

from sympy import limit, sin, oo from sympy.abc import x limit(sin(x)/x, x, 0) # 1 limit(1/x, x, 0) # oo（默认 dir='+'，从右侧） limit(1/x, x, 0, dir="-") # -oo（从左侧） limit(1/x, x, 0, dir='+-') # zoo（双侧极限，发散为复无穷） limit(1/x, x, oo) # 0

参数说明（与 limit 源码 docstring 一致）：

• e：要求极限的表达式；
• z：极限变量（其他符号视为常量；多元极限不支持）；
• z0：趋于的点，可以是任意表达式，包括oo与-oo；
• dir：方向，"+"表示右极限（z→z0+），"-"表示左极限（z→z0-），"+-"表示双侧。当z0为无穷时，dir由无穷的方向决定（即对oo取dir="-"）。

limit()的实现（limits.py）只是Limit(e, z, z0, dir).doit(deep=False)的语法糖，返回的是Limit对象求值后的结果。

未求值极限类Limit

Limit表示一个"尚未求值"的极限表达式，可以直接构造：

from sympy import Limit, sin from sympy.abc import x Limit(sin(x)/x, x, 0) # Limit(sin(x)/x, x, 0, dir='+') Limit(1/x, x, 0, dir="-") # Limit(1/x, x, 0, dir='-')

从 构造逻辑 可以看出几个重要约束：当z0本身含极限变量时抛出NotImplementedError（不支持趋于动点）；dir只能是'+'、'-'、'+-'三者之一。其doit()方法（limits.py）内部的分支顺序本身就揭示了求解路径：

1. 双侧极限分别求左右极限并比较；
2. 趋于无穷时做代换、规范化符号；
3. 对Float浮点做nsimplify精确化，避免舍入误差导致条件分支误判；
4. 若表达式在z0处亚纯（meromorphic），直接用leadterm提取首项判定；
5. 兜底调用gruntz()（见下节）；若 Gruntz 也失败（PoleError/ValueError），再退回heuristics()逐项求值。

这里的关键事实是文档明确指出的：极限计算的真正主力是gruntz()函数，它实现了 Gruntz 算法。limit()只是"先启发式、后 Gruntz"的调度层。相关的启发式逻辑可在 heuristics() 中查看——它逐项计算各参数的极限再合并结果，只对简单情形快而有效。

Gruntz 算法：极限计算背后的理论

函数增长速度的序关系

Gruntz 算法的第一步是建立函数之间的"增长序"。设f(x)、g(x)是两个实值函数，且当x→∞时都趋于+∞。文档定义了支配关系：

若对任意a, b ∈ ℝ>₀都有lim_{x→∞} f(x)^a / g(x)^b = ∞，则称f(x)支配g(x)，记作f(x) ≻ g(x)。若f(x)与g(x)互不支配，则称它们属于同一个可比较类（comparability class），记作f(x) ≍ g(x)。

该定义可通过f ~ 1/f、f ~ -f自然地推广到趋于 0 或±∞的函数。源码 gruntz.py 的模块 docstring 给出了等价的工程化判据：计算L = lim log|f(x)| / log|g(x)|（x→∞），L=±∞则f ≻ g，L=0则f ≺ g，L为有限非零值则同属一类。

文档给出的一组经典例子（在源码 docstring 中也重复出现）：

• e^x ≻ x^m
• e^{x²} ≻ e^{mx}
• e^{e^x} ≻ e^{x^m}
• x^m ≍ x^n（任何正幂同属一类）
• e^{x+1/x} ≍ e^{x+log x} ≍ e^x

也就是说：2 < x < exp(x) < exp(x**2) < exp(exp(x))，而2 ~ 3 ~ -5、x ~ x**2 ~ x**3 ~ 1/x、exp(x) ~ exp(-x) ~ exp(2x)。最低的可比较类被设定为2~3~-5（即常数类），log(x)低于该类，不参与比较。

支撑极限判定的定理

文档紧接着给出两条核心结论：

定理一（主项支配求和）：设ω, g₁, g₂, …是x的函数，lim ω = 0且ω ≻ gᵢ对所有i成立；设c₁ < c₂ < …为严格递增的实数。则

lim_{x→∞} Σᵢ gᵢ ω^{cᵢ} = lim_{x→∞} g₁ ω^{c₁}

即求和式的极限由指数最小的首项决定，其余项都更快地趋于零。

推论（对g·ω^c的三种情形）：

• c > 0时，lim g·ω^c = 0；
• c < 0时，lim g·ω^c = ±∞，符号由g的（最终）符号决定；
• c = 0时，lim g·ω^0 = lim g。

这两条正是 limitinf() 的实现依据：算出展开式的首项(c0, e0)后，依据e0的符号直接返回0、±∞或递归地求lim c0。

四步求解策略

据此，文档给出了计算lim_{x→∞} f(x)的标准流程：

1. 找出f(x)的 MRV 集（most rapidly varying subexpressions，最快变化子表达式集）：从f的所有子表达式中，选出在≻关系下极大的那些元素；
2. 选取ω：选择一个与 MRV 集元素同属一个可比较类、且lim ω = 0的函数（例如 MRV 集为{eˣ, e²ˣ}时，可取ω = e⁻ˣ）；
3. 按ω展开级数：把f写成f = c₀ω^{e₀} + c₁ω^{e₁} + … + O(ω^{eₙ})（e₀<e₁<…<eₙ，c₀≠0），使之满足定理一的前提；
4. 应用定理并递归：由首项判定极限，必要时对g₁(x)递归执行以上步骤。

源码 gruntz.py 的注释精确对应这四步：求 MRV 集用mrv()，选ω并改写用rewrite()，求首项用mrv_leadterm()/leadterm()，最终收敛在limitinf()中完成。

实现中的关键函数

文档 Reference 节逐一列出了 gruntz.py 的公开接口，对应实现如下：

• gruntz(e, z, z0, dir="+")：算法入口。z0可为任意表达式（含oo/-oo）；dir="+"求右极限、dir="-"求左极限（z0无穷时dir无意义）。内部先把一切极限通过代换（z0 + 1/z、z0 - 1/z、-z）归约为z→∞，再交给limitinf()，最后用rewrite('intractable')把结果美化回熟悉的形式。注意：gruntz()不接收非符号作为第二个参数（第 679 行）。
• compare(a, b, x)：计算L = lim log(a)/log(b)，返回"<"、"="或">"。对exp类表达式会直接取指数部分（log(exp(…))必须在此化简以保证算法终止）。
• mrv(e, x)：返回一个SubsSet，包含e的所有 MRV 子表达式及其重写信息。对Add/Mul递归拆分比较（mrv_max1/mrv_max3），对log、exp、一般函数、导数分别处理；未定义函数会抛ValueError，多变量函数（如BesselJ(x, x)）与导数暂未实现（NotImplementedError）。
• mrv_leadterm(e, x)：返回e的首项(c₀, e₀)。内部调用rewrite()后对w求leadterm，失败时逐步增大_eval_nseries的阶数重试；还用ilcm对分母取最小公倍数消除非整数幂展开的困难。
• limitinf(e, x)：计算x→∞的极限。先做rewrite('tractable')、powdenest等规范化，再依e₀符号返回0/±∞/递归求lim c₀。
• sign(e, x)：返回e在x→∞时的最终符号（1/0/-1）；对反复变号的函数（如sin(x)）结果未定义。实现先查假设（is_positive等），再对Mul、exp、Pow、log结构递归，最后兜底走mrv_leadterm。
• rewrite(e, Omega, x, wsym)：算法中最复杂的一步。docstring 明确说"当算法失败时，bug 通常在级数展开（即 SymPy 本身）或 rewrite 中"。它借助 build_expression_tree() 按依赖高度排序 MRV 表达式，把e改写为关于w（趋零变量）的级数形式，并同步计算出log(w)。
• SubsSet：MRV 表达式→哑变量的字典，附带rewrites映射。类 docstring 解释了为何不能用普通的subs()：它"太聪明"，会错误化简（文档中的-1/w + exp(exp(p)*exp(-exp(-p))/(1 - 1/p))例子说明普通subs无法正确处理嵌套 MRV 表达式）。
• mrv_max1 / mrv_max3：两组 MRV 集合取"更大"者；同可比较类时取并集（gruntz.py）。

文档未展开的四个细节与调试手段

文档"Notes"节诚实指出，前面的叙述略过了四个关键问题：

1. 如何判定f ≻ g、g ≻ f还是f ≍ g？（→ 由compare()实现）
2. 如何求 MRV 集？（→ 由mrv()实现）
3. 如何计算级数展开？（→ 依赖 SymPy 的leadterm/_eval_nseries）
4. 算法为何必然终止？（→ 由 Gruntz 博士论文证明）

gruntz.py 的模块 docstring 还提供了调试方法：Gruntz 算法高度递归，难以在调试器里追踪，可通过设置环境变量SYMPY_DEBUG=True开启漂亮的调试打印，例如：

SYMPY_DEBUG=True python -m isympy

然后执行limit(sin(x)/x, x, 0)，控制台会输出完整的递归树：

limitinf(_x*sin(1/_x), _x) = 1 +-mrv_leadterm(_x*sin(1/_x), _x) = (1, 0) | +-mrv(_x*sin(1/_x), _x) = set([_x]) | | +-mrv(_x, _x) = set([_x]) | | +-mrv(sin(1/_x), _x) = set([_x]) ... +-sign(0, _x) = 0 +-limitinf(1, _x) = 1

逐行核对哪一步出错，再回到对应源码函数定位问题。这得益于sympy.utilities.misc.debug_decorator（见 gruntz.py）对mrv、sign、limitinf、mrv_leadterm、rewrite等函数的装饰。此外，算法实现几乎是 Gruntz 博士论文中 Maple 代码的直接改写（gruntz.py 注释），论文包含完整描述与大量示例，是深入理解算法细节的首选资料。

更直观的级数展开：series()包装

文档指出x*cos(x)这种写法比(x*cos(x)).series(x)更直观，因此模块提供了一个围绕Basic.series()的包装函数 series()：

from sympy import Symbol, cos, series x = Symbol('x') series(cos(x), x) # 1 - x**2/2 + x**4/24 + O(x**6)

签名与参数（series.py）：

• expr：要展开的表达式；
• x：展开变量；
• x0：展开点，取值范围-oo到oo；
• n：展开到第n项（默认 6）；
• dir：方向，"+"为x→x0+，"−"为x→x0−；当x0为无穷时方向由无穷本身决定（对oo即dir="-"）。

例如在x=2处对tan(x)展开：

from sympy import series, tan from sympy.abc import x series(tan(x), x, 2, 6, "+") # tan(2) + (1 + tan(2)**2)*(x - 2) + ... + O((x - 2)**6, (x, 2)) series(tan(x), x, 2, 3, "-") # tan(2) + (2 - x)*(-tan(2)**2 - 1) + ... + O((x - 2)**3, (x, 2))

注意：n=oo会抛出TypeError（源码 docstring 中即有此示例），因为项数必须是整数。实现上 series() 只是sympify(expr)后调用expr.series(x, x0, n, dir)——完整语义请参见 Expr.series 的 docstring。

阶项追踪：Order类

级数展开不可避免要记录"截断误差"，本模块通过 Order（别名O）自动完成阶项的追踪。文档示例：

from sympy import Symbol, Order x = Symbol('x') Order(x) + x**2 # O(x) Order(x) + 1 # 1 + O(x)

Order用大 O 记号表征函数的极限行为：g(x) = O(f(x))（x→a）当且仅当存在δ>0与M>0，使得|x-a|<δ时|g(x)| ≤ M|f(x)|，等价于limsup |g(x)/f(x)| < ∞。以sin(x)在 0 处展开为例，sin(x) = x - x³/3! + O(x⁵)，此时O(x⁵) = x⁵/5! - x⁷/7! + …，且lim |O(x⁵)/x⁵| = 1/5!（order.py）。

Order 的构造逻辑 还包含若干自动规范化规则：

• 不传变量时默认使用表达式全部自由符号，展开点默认为 0；
• O(expr*f(x), x)规约为O(f(x), x)，O(expr, x)规约为O(1)，O(0, x)规约为0；
• O(f(x), x)会自动取f(x).as_leading_term(x)的首项；
• 支持多元O(f(x,y), x, y)，假定各变量取极限可交换；支持在oo、-oo等展开点（内部做1/Dummy()代换归化到 0）；
• 提供contains()方法判定阶项包含关系，如O(x) in O(1, x)为 True（order.py）：
  • O(1) in O(1, x)→ True，O(1, x) in O(1)→ False
  • O(x) in O(1, x)→ True，O(x**2) in O(x)→ True
• 算术规则：O(x)*x→O(x**2)，O(x) - O(x)→O(x)，O(cos(x))→O(1)，O(cos(x), (x, pi/2))→O(x - pi/2, (x, pi/2))（order.py 示例）。

Order的removeO()返回 0、getO()返回自身，是提取/丢弃阶项的常用工具；_eval_nseries返回自身，保证它不参与再次展开。

级数加速：Richardson 外推与 Shanks 变换

文档中 "Series Acceleration" 一节标注为 TODO，但 acceleration.py 中的两个函数实际上已完整实现，参考书目为 Bender & Orszag《Advanced Mathematical Methods for Scientists and Engineers》（Shanks 变换见 pp. 368-375，Richardson 外推见 pp. 375-377）。

richardson(A, k, n, N)：用项A(n), A(n+1), …, A(n+N+1)做 Richardson 外推，近似lim_{k→∞} A(k)。经验上取N ≈ 2n效果较好。经典例子是用极限定义算e：e = (1+1/n)**n收敛极慢，n=100时只有两位有效数字（2.7048138294），而用richardson(e, n, 10, 20)一次即可得到2.7182818285。对 ζ(2) 级数1/k²的部分和取前 100 项仅得1.6349839002，Richardson 外推则给出1.6449340668，与精确值π²/6一致。

shanks(A, k, n, m=1)：n 项 Shanks 变换S(A)(n)；m>1时做 m 重递归变换S(S(…S(A)…))(n)。它对在极点/奇点附近收敛慢的 Taylor 级数特别有效，例如交错调和级数求log(2)：前 100 项只有0.6881721793，shanks(A, n, 25)得到0.6931396564，shanks(A, n, 25, 5)达到0.6931471806，与log(2)的精确值0.6931471805599453…高度吻合。

两个函数的实现均基于文献中的递推/插值公式：Richardson 使用(n+j)^N与二项式系数的加权和（acceleration.py），Shanks 使用三对角(z*x - y²)/(z + x - 2y)递推表（acceleration.py）。这两个函数的正确性在 test_series.py 中有断言覆盖。

留数计算：residue()

文档 "Residues" 一节同样标注 TODO，但 residues.py 中的 residue(expr, x, x0) 已实现。留数定义为表达式在x=x0处幂级数展开中1/(x-x0)项的系数：

from sympy import Symbol, residue, sin x = Symbol("x") residue(1/x, x, 0) # 1 residue(1/x**2, x, 0) # 0 residue(2/sin(x), x, 0) # 2

实现原理（源码注释）：当前实现基于级数展开，先平移x→x+x0，用nseries在n ∈ (0,1,2,4,8,16,32)中逐步取到足够的展开阶，collect后逐项扫描，凡是1/x项就把系数累加进结果；若出现无法识别的项则抛NotImplementedError。对于一般函数可参考 Bronstein《Symbolic Integration I》5.6 节，纯有理函数则有基于结式（resultant）的简单算法（见 residues.py 的注释）。

test_residues.py 提供了丰富的验证用例，覆盖了各种情形：

residue(1/(x**2 + 1), x, I) # -I/2 residue(1/(x**2 + 1), x, -I) # I/2 residue(1/sin(x)**5, x, 0) # Rational(3, 8) residue(1/(x**4 + 1), x, exp(I*pi/4)).equals(-(Rational(1, 4) + I/4)/sqrt(2))

以及"展开点不是极点则留数为 0"的性质：residue(1/x, x, 1) == 0。留数定理（residues.py 引用）使其成为复分析积分与级数求和的核心工具。

从文档到实战：使用建议

综合文档与源码，可以总结出几条实用建议：

1. 日常求极限直接使用limit()。它先走启发式快速路径，失败才进入 Gruntz 算法，绝大多数情形开箱即用。
2. 需要研究算法行为时，把dir方向与Limit类配合使用：Limit让你保留"未求值"形式的极限（便于符号操作），.doit()触发求值；dir='+-'会校验左右极限是否一致，不一致时抛出带左右值的ValueError。
3. 级数展开优先用series()函数式写法，并善用Order做截断管理；需要丢弃阶项时用.removeO()。
4. gruntz()直接调用的适用前提：第二个参数必须是Symbol（gruntz.py），文档也提示它"只在更快的limit()失败时才常被用到"。
5. 当 Gruntz 结果异常时，用SYMPY_DEBUG=True开启递归打印，从输出树中定位出错的mrv/rewrite环节，再回到 gruntz.py 对应函数排查。
6. 数值上收敛很慢的极限/级数，可尝试richardson/shanks加速；符号留数计算则用residue()，注意其实现基于nseries展开，复杂奇点（如本性奇点）可能需要更高阶展开。

如需验证行为或深入学习，可查阅 sympy/series/tests/ 下的测试文件：Gruntz 算法与极限的断言集中在 test_gruntz.py 与 test_limits.py，级数展开见 test_series.py 与 test_nseries.py，阶项语义见 test_order.py，留数见 test_residues.py。

【免费下载链接】sympyA computer algebra system written in pure Python项目地址: https://gitcode.com/GitHub_Trending/sy/sympy