如何不用 SymPy 高效计算 arctan 高阶导数并生成 Taylor 级数？

2026-11-04数据小许

可以不使用 SymPy。对 arctan 求导时，保留分母 (1+x**2)**n，只递推分子，就能避免当前实现中分母次数的指数级增长。如果最终只需要 Taylor 系数，也可以直接进行幂级数运算，无需构造完整的高阶导数表达式。

为什么当前算法越来越慢？

frac_poly_dif 每次都会将整个分母平方，却没有消去分子和分母的公共因子。从 arctan 的一阶导数开始，未经约分的分母依次是：

(1+x**2)**1
(1+x**2)**2
(1+x**2)**4
(1+x**2)**8
...

因此，计算 arctan 的第 n 阶导数时，当前表示中的分母为 (1+x**2)**(2**(n-1))，展开后的次数是 2**n。

这些导数可以统一写成：

D**n arctan(x) = P_n(x)/(1+x**2)**n,    n >= 1

其中 P_n 是次数为 n-1 的多项式。你给出的第六阶导数，分母只需保留为 (1+x**2)**6。

问题出在未经约分的表达式不断膨胀。仅删除零系数或优化 Counter 的遍历，无法解决这个增长问题。

需要完整导数时：递推分子，保留分母的幂

设：

q(x) = 1+x**2
D**n arctan(x) = P_n(x)/q(x)**n

对这个形式求导并消去公共因子：

D(P_n/q**n)
    = (q*P_n' - n*q'*P_n)/q**(n+1)

由于 q'(x)=2*x，递推式为：

P_1(x) = 1
P_(n+1)(x) = (1+x**2)*P_n'(x) - 2*n*x*P_n(x)

分子中的一项 c*x**k 在下一步的贡献为：

k*c*x**(k-1) + (k-2*n)*c*x**(k+1)

因此不需要进行一般的多项式乘法，直接移动指数、累加系数即可：

from collections import Counter


def arctan_derivative_numerator(n):
    """返回 arctan 的第 n 阶导数的分子，分母为 (1+x**2)**n。"""
    if n < 1:
        raise ValueError("n must be >= 1")

    p = Counter({0: 1})

    for order in range(1, n):
        result = Counter()

        for degree, coefficient in p.items():
            if degree > 0:
                result[degree - 1] += degree * coefficient

            result[degree + 1] += (degree - 2 * order) * coefficient

        p = Counter({
            degree: coefficient
            for degree, coefficient in result.items()
            if coefficient != 0
        })

    return p


numerator = arctan_derivative_numerator(16)
print(sorted(numerator.items()))
# 完整导数为该分子除以 (1+x**2)**16。

前三个分子分别为：

P_1(x) = 1
P_2(x) = -2*x
P_3(x) = 6*x**2 - 2

这个递推只使用整数运算，不需要计算多项式最大公因式。计算到第 N 阶时，累计需要 O(N**2) 次量级的系数运算。这个估算假定每次整数运算的成本固定，实际成本还受大整数位数的影响。

只需要 Taylor 系数时：直接递推级数

在零点展开 arctan，可以先展开它的一阶导数。设：

1/(1+x**2) = b_0 + b_1*x + b_2*x**2 + ...

两边乘以 1+x**2，比较相同次数的系数，得到：

b_0 = 1
b_1 = 0
b_k = -b_(k-2),    k >= 2

再逐项积分，利用 arctan(0)=0，得到：

arctan(x) = x - x**3/3 + x**5/5 - x**7/7 + ...

这个过程也可以写成一般的有理函数级数除法。设分子、分母和商的系数分别为 p_n、q_n 和 b_n，且 q_0 != 0，则：

b_n = (p_n - sum(q_k*b_(n-k), k=1..n))/q_0

超出多项式次数的系数按零处理。下面的实现只依赖 Python 标准库，用 Fraction 保留精确的有理数结果：

from fractions import Fraction


def rational_series(p, q, order):
    """返回 p(x)/q(x) 在 x=0 处到 order 次的幂级数系数。

    p、q 均按常数项、一次项、二次项……的顺序传入。
    """
    if order < 0:
        raise ValueError("order must be >= 0")

    p = [Fraction(value) for value in p]
    q = [Fraction(value) for value in q]

    if not q or q[0] == 0:
        raise ValueError("the denominator must have a nonzero constant term")

    coefficients = []

    for n in range(order + 1):
        value = p[n] if n < len(p) else Fraction(0)

        for k in range(1, min(n, len(q) - 1) + 1):
            value -= q[k] * coefficients[n - k]

        coefficients.append(value / q[0])

    return coefficients


# arctan'(x) = 1/(1+x**2)，先求到 x**15 的系数。
derivative_coefficients = rational_series([1], [1, 0, 1], 15)

# 逐项积分，得到 arctan 到 x**16 的 Taylor 系数。
taylor_coefficients = [Fraction(0)] + [
    coefficient / (degree + 1)
    for degree, coefficient in enumerate(derivative_coefficients)
]

print([
    (degree, str(coefficient))
    for degree, coefficient in enumerate(taylor_coefficients)
    if coefficient != 0
])

输出为：

[(1, '1'), (3, '-1/3'), (5, '1/5'), (7, '-1/7'),
 (9, '1/9'), (11, '-1/11'), (13, '1/13'), (15, '-1/15')]

对于这里固定的二次分母，每个系数最多依赖前两个系数。生成 N 个系数需要 O(N) 次量级的有理数运算。同样，这个估算没有将大整数运算的位数成本视为变量。

如果需要在其他点 a 展开，令 t=x-a，则：

arctan'(a+t) = 1/((1+a**2) + 2*a*t + t**2)

把分母系数换成 [1+a**2, 2*a, 1]，即可生成关于 t 的导数级数，再逐项积分。积分常数 arctan(a) 需要另外提供，上述递推不会自动求出这个常数。若要继续使用 Fraction 进行精确计算，a 应用整数或 Fraction 表示。

实现时需要注意什么？

注意区分导数阶数。初始的 1/(1+x**2) 已经是 arctan 的一阶导数，再调用 16 次求导，得到的是 arctan 的第 17 阶导数。
第 n 阶 Taylor 系数是 f^(n)(a)/n!。上面的级数除法和积分直接生成系数，不需要额外再除一次阶乘。
保留负系数。使用 Counter 时，应继续逐项累加；它的二元 +、- 运算会丢弃非正计数，不适合直接用于有符号多项式的加减法。
arctan 的 Taylor 系数包含 1/3 等有理数。Fraction 使用整数分子、分母表示它们，可以避免浮点舍入。
计算系数后，还要判断级数的收敛范围。零点处的 arctan 级数在 |x|<1 内收敛到 arctan，在实数端点 x=±1 也收敛，但较慢；不能直接用于 |x|>1 的求值。对于其他函数，也需要单独判断解析性和收敛范围。

如果要通用地生成更多函数的 Taylor 级数，可以进一步实现截断幂级数的加、减、乘、除、求导和积分，并为支持的函数增加相应的系数递推规则。当前的多项式商法则实现只覆盖有理函数，不能直接处理所有超越函数。