SciPy 与 SymPy

发表于 更新于 分类于 Waline： 阅读次数： 本文字数： 2.6k 阅读时长 ≈ 2 分钟

SciPy 与 SymPy：数值算结果，符号推公式，各管各的

做科学计算的时候，你经常会遇到两种需求：

  • 算一个具体的数：比如 ∫₀¹ sin(x) dx = 0.4596...
  • 推一个数学公式：比如 ∫ x·sin(x) dx = -x·cos(x) + sin(x)

前者是数值计算（SciPy 的活），后者是符号运算（SymPy 的活）。数值给结果，符号给表达式。 两套工具，解决两类问题。

数值积分：函数没原函数，或者数据来自实验

很多工程问题需要算积分，但被积函数可能没有初等原函数，或者数据是测量来的（离散点，没有表达式）。这时候用数值积分。

一维定积分：quad（最常用）

from scipy.integrate import quad
import numpy as np

# 算 ∫₀¹ sin(x) dx
result, error = quad(np.sin, 0, 1)
print(f"{result:.10f}")   # 0.4596976941

quad 返回两个值：积分结果和误差估计。误差估计别看小，它告诉你结果可不可信。 如果误差比结果还大，说明积分没收敛。

带参数的积分：

def integrand(x, a, b):
    return np.exp(-a * x) * np.sin(b * x)

result, error = quad(integrand, 0, np.inf, args=(2.0, 3.0))

二重积分：dblquad（注意积分顺序）

from scipy.integrate import dblquad

# 算半球体积：∫∫ √(1-x²-y²) dy dx
def f(x, y):
    return np.sqrt(max(0, 1 - x**2 - y**2))

volume, error = dblquad(
    f,
    -1, 1,
    lambda x: -np.sqrt(1-x**2),
    lambda x: np.sqrt(1-x**2)
)
print(f"{volume:.6f}")   # 2.094395（= 2π/3）

注意： dblquad 的积分顺序跟数学表达式是反的。数学上写 ∫∫ f(x,y) dx dy，dblquad 里是先积 y 再积 x。搞反了结果完全不对。

三重及以上：nquad

from scipy.integrate import nquad

def integrand(x1, x2, x3, x4):
    return 1 if x1**2 + x2**2 + x3**2 + x4**2 <= 1 else 0

ranges = [(-1, 1)] * 4
volume, error = nquad(integrand, ranges)

离散数据积分：数据是点，不是函数

实验数据往往是一堆 (x, y) 点，没有表达式。

from scipy.integrate import simpson, trapezoid

x = np.linspace(0, np.pi, 100)
y = np.sin(x)

area_simpson = simpson(y, x)    # 更精确，但要求均匀采样
area_trapz = trapezoid(y, x)    # 更通用，不均匀也能算

选哪个？

数据特点 用哪个
有解析函数 quad（精度高，自动控制误差）
均匀采样的离散点 simpson（精度更高）
不均匀采样的离散点 trapezoid（只能用这个）

常微分方程：odeint 模拟动态系统

物理系统常用微分方程描述，但解析解往往不存在。odeint 用数值方法（Runge-Kutta）近似求解。

洛伦茨吸引子（混沌的经典例子）：

from scipy.integrate import odeint

def lorenz(state, t, sigma, rho, beta):
    x, y, z = state
    return [
        sigma * (y - x),
        x * (rho - z) - y,
        x * y - beta * z
    ]

t = np.linspace(0, 30, 3000)
track = odeint(lorenz, [0.0, 1.0, 0.0], t, args=(10.0, 28.0, 8.0/3.0))

几个实用建议：

  • 时间步长要够密。 变化快的系统步长太大会不稳定。
  • 刚性系统（快慢时间尺度共存）用 odeint 的 rtol 和 atol 调精度，或者换 solve_ivp 的 BDF 方法。
  • odeint 不支持事件检测（比如碰撞、触顶），这种情况用 solve_ivp。

符号运算：让计算机帮你推公式（SymPy）

数值计算算的是”数”，符号运算算的是”式子”。

from sympy import symbols, diff, integrate

x = symbols('x')
expr = x**3 + 2*x + 1

print(diff(expr, x))          # 3*x**2 + 2
print(integrate(x*sin(x), x)) # -x*cos(x) + sin(x)

符号的假设条件会影响化简结果：

from sympy import sqrt, simplify

x = symbols('x')
simplify(sqrt(x**2))        # sqrt(x**2) —— 没化简

x = symbols('x', positive=True)
simplify(sqrt(x**2))        # x —— 因为知道 x 是正数

解微分方程（dsolve）：

from sympy import Function, dsolve, Eq

f = Function('f')
x = symbols('x')

eq = Eq(f(x).diff(x) - f(x), 0)
dsolve(eq, f(x))   # f(x) = C1*exp(x)

实际场景：推导球体体积公式

from sympy import symbols, integrate, sqrt

r, x = symbols('r', positive=True)

# 圆的面积：2∫₀ʳ √(r²-x²) dx
circle_area = 2 * integrate(sqrt(r**2 - x**2), (x, 0, r))

# 球体体积：∫_{-r}^{r} π*(r²-x²) dx
volume = integrate(circle_area.subs(r, sqrt(r**2 - x**2)), (x, -r, r))
volume.simplify()   # 4*pi*r**3/3