SymPy 微积分实战

参考链接

Calculus — SymPy documentation

环境准备

开始计算前，先导入必要的模块并定义符号变量：

from sympy import symbols, diff, integrate, limit, series, sin, exp, ln, oo

# 定义符号
x, y, z = symbols('x y z')

微分

使用 diff 进行符号求导：

# 对 f(x) = x^2 求导
f_prime = diff(x**2, x)  # 结果: 2*x

# 对 f(x, y) = x^2 + xy 求偏导
partial_f_x = diff(x**2 + x*y, x)  # 结果: 2*x + y
partial_f_y = diff(x**2 + x*y, y)  # 结果: x

积分

支持不定积分和定积分计算：

# 对 f(x) = x^2 求不定积分
integral_f = integrate(x**2, x)  # 结果: x**3/3

# 计算 f(x) = x^2 在 [0, 1] 区间上的定积分
definite_integral_f = integrate(x**2, (x, 0, 1))  # 结果: 1/3

极限

计算变量趋近特定值时的函数极限：

# x -> 0 时 sin(x)/x 的极限
lim_sin_x = limit(sin(x)/x, x, 0)  # 结果: 1

# x -> 无穷大时 (1+1/x)^x 的极限
lim_exp = limit((1+1/x)**x, x, oo)  # 结果: E

级数展开

将函数在指定点展开为泰勒级数：

# 保留低于 x**5 的幂；O(x**5) 表示余项阶数
series_exp = series(exp(x), x, 0, 5)
# 结果: 1 + x + x**2/2 + x**3/6 + x**4/24 + O(x**5)

# 保留低于 x**4 的幂；O(x**4) 表示余项阶数
series_ln = series(ln(1+x), x, 0, 4)
# 结果: x - x**2/2 + x**3/3 + O(x**4)

微分方程求解

SymPy 也能符号化求解微分方程，以下是一个二阶常系数非齐次方程的例子：

from sympy import Function, dsolve, Eq, Derivative

# 定义未知函数 f(x) 及微分方程: f''(x) - 2f'(x) + f(x) = sin(x)
f = Function('f')
diffeq = Eq(Derivative(f(x), x, x) - 2*Derivative(f(x), x) + f(x), sin(x))

# 求解
sol = dsolve(diffeq, f(x))

假设、精确值与结果检查

在已有 SymPy 的 Python 环境中按顺序运行这些片段。它们演示符号结果，不保证每个积分或微分方程都有闭式解。未求值的 Integral 表示积分尚未算出，不表示结果为零。官方微积分教程说明了这类对象和级数的余项记号。

若不指定假设，符号就没有相应约束。例如对负实数 $r$，不能把 $\sqrt{r^2}$ 化简为 $r$，正确结果是 $|r|$。需要精确代数时，应使用 SymPy 的精确数值：

from sympy import Rational, sqrt, simplify, Abs

r = symbols('r', real=True)
p = symbols('p', positive=True)
assert sqrt(r**2) == Abs(r)
assert sqrt(p**2) == p
assert integrate(x**2, (x, 0, 1)) == Rational(1, 3)
print(Rational(1, 3).evalf(20))

先写 Python 的 1/3 会得到浮点近似，之后提高显示精度无法恢复已丢失的精确性。integrate(x**2, x) 返回一个原函数；完整的原函数族还包含任意加法常数。可以在定义域内求导来检查原函数。

在有限点求极限时，如果左右两侧可能不同，就应明确指定方向：

assert limit(1/x, x, 0, dir='+') == oo
assert limit(1/x, x, 0, dir='-') == -oo
assert simplify(diff(integral_f, x) - x**2) == 0

这里左右极限不同，因此双侧极限不存在。以 O(x**5) 结尾的泰勒展开保留低于五次的幂；删除余项符号只得到局部多项式近似，不是全局恒等式。

前面的微分方程通解为 $f(x)=(C_1+C_2x)e^x+\tfrac12\cos x$。需要两个初始条件才能确定两个常数。应同时检查方程与初始条件，而非只看打印形式：

from sympy import cos

C1, C2 = symbols('C1 C2')
g = (C1 + C2*x)*exp(x) + cos(x)/2
assert simplify(diff(g, x, 2) - 2*diff(g, x) + g - sin(x)) == 0
sol_ivp = dsolve(diffeq, f(x), ics={f(0): 0, diff(f(x), x).subs(x, 0): 0})
assert simplify(sol_ivp.rhs - ((x - 1)*exp(x) + cos(x))/2) == 0

我的阅读 · 打开过的笔记会记在当前浏览器里。