如何通过编程方法推导泰勒级数？以反正切函数为例的实践分析

我最近在琢磨怎么不靠现成的泰勒级数公式，而是通过反复求导来推导它，进而计算像反正切这样的超越函数，这里把我的思路和实现过程分享出来：

核心思路 #

从物理里的位置与导数的关系，我联想到这和泰勒级数本质是一回事——某一时刻的位置可以通过导数展开，这完全对应泰勒级数的形式。所以我想试试：从目标函数的导数出发，一步步计算高阶导数，进而推导出泰勒展开的系数，而不是直接用已知的泰勒公式。这里我选了反正切函数来测试，它的一阶导数是1/(x²+1)。

多项式运算基础规则 #

要实现分式多项式的求导，得先搞定基础的多项式运算：

• 多项式求导：对每个(系数, 指数)对，新系数是系数×指数，新指数是指数-1；指数为0的项直接丢弃。
• 多项式减法：相同指数的项系数相减，对方有但自己没有的项，取负系数加入结果。
• 多项式乘法：每一项两两相乘，系数相乘、指数相加，最后合并相同指数的系数。
• 分式多项式求导：严格遵循商的导数法则：(f/g)' = (f'g - g'f)/g²

代码实现 #

我用Python的collections.Counter来表示多项式（键是指数，值是对应系数），写了一套完整的多项式运算和求导工具：

import sympy
from collections import Counter

def poly_add(poly1: Counter, poly2: Counter) -> Counter:
    result = poly1.copy()
    for exp, coef in poly2.items():
        result[exp] += coef
    return result

def poly_sub(poly1: Counter, poly2: Counter) -> Counter:
    result = poly1.copy()
    for exp, coef in poly2.items():
        result[exp] -= coef
    return result

def poly_mul(poly1: Counter, poly2: Counter) -> Counter:
    result = Counter()
    for exp1, coef1 in poly1.items():
        for exp2, coef2 in poly2.items():
            result[exp1 + exp2] += coef1 * coef2
    return result

def poly_dif(poly: Counter) -> Counter:
    result = poly.copy()
    result.pop(0, None)
    return Counter({exp - 1: exp * coef for exp, coef in result.items()})

def frac_poly_dif(poly_pair: tuple[Counter, Counter]) -> tuple[Counter, Counter]:
    num_poly, den_poly = poly_pair
    return poly_sub(
        poly_mul(poly_dif(num_poly), den_poly),
        poly_mul(poly_dif(den_poly), num_poly)
    ), poly_mul(den_poly, den_poly)

x = sympy.symbols("x")

def sympy_simplify(poly_pair: tuple[Counter, Counter]):
    num_poly, den_poly = poly_pair
    terms1 = iter(list(num_poly.items()))
    a, b = next(terms1)
    expr1 = f"{b}*x**{a}"
    for a, b in terms1:
        expr1 = f"{expr1}{'+-'[b < 0]}{abs(b)}*x**{a}"
    expr1 = expr1.replace("*x**0", "")
    terms2 = iter(list(den_poly.items()))
    a, b = next(terms2)
    expr2 = f"{b}*x**{a}"
    for a, b in terms2:
        expr2 = f"{expr2}{'+-'[b < 0]}{abs(b)}*x**{a}"
    expr2 = expr2.replace("*x**0", "")
    return eval(f"sympy.simplify(({expr1})/({expr2}))")

测试结果与问题 #

以反正切函数的初始导数(Counter({0: 1}), Counter({0: 1, 2:1}))为例，计算前5阶导数：

In [2]: poly_pair = (Counter({0: 1}), Counter({0: 1, 2:1}))
In [3]: derivatives = [poly_pair, *[(poly_pair := frac_poly_dif(poly_pair)) for _ in range(5)]]
In [4]: derivatives
Out[4]: [(Counter({0: 1}), Counter({0: 1, 2: 1})), (Counter({1: -2}), Counter({2: 2, 0: 1, 4: 1})), (Counter({4: 6, 2: 4, 0: -2}), Counter({4: 6, 2: 4, 6: 4, 0: 1, 8: 1})), (Counter({3: 72, 5: 48, 1: 24, 11: -24, 7: -48, 9: -72}), Counter({8: 70, 6: 56, 10: 56, 4: 28, 12: 28, 2: 8, 14: 8, 0: 1, 16: 1})), (Counter({20: 6864, 22: 3984, 18: 1320, 24: 1080, 26: 120, 0: 24, 2: 24, 4: -1200, 6: -7920, 16: -19800, 8: -25080, 14: -47520, 10: -48312, 12: -60192}), Counter({16: 12870, 14: 11440, 18: 11440, 12: 8008, 20: 8008, 10: 4368, 22: 4368, 8: 1820, 24: 1820, 6: 560, 26: 560, 4: 120, 28: 120, 2: 16, 30: 16, 0: 1, 32: 1})), (Counter({29: 11054352000, 27: 10127212800, 31: 10127212800, 25: 7764790800, 33: 7764790800, 23: 4935652800, 35: 4935652800, 21: 2549632800, 37: 2549632800, 19: 1026168000, 39: 1026168000, 17: 288116400, 41: 288116400, 15: 31574400, 43: 31574400, 1: -720, 57: -720, 3: -16320, 55: -16320, 5: -172320, 53: -172320, 7: -1110720, 51: -1110720, 9: -4758000, 49: -4758000, 11: -13353600, 47: -13353600, 13: -18657600, 45: -18657600}), Counter({32: 601080390, 30: 565722720, 34: 565722720, 28: 471435600, 36: 471435600, 26: 347373600, 38: 347373600, 24: 225792840, 40: 225792840, 22: 129024480, 42: 129024480, 20: 64512240, 44: 64512240, 18: 28048800, 46: 28048800, 16: 10518300, 48: 10518300, 14: 3365856, 50: ...})]

结果是正确的，但问题也很明显：计算成本会随着导数阶数升高迅速爆炸，分子分母的项数、系数大小都会指数级增长，这是接下来需要优化的方向。

内容的提问来源于stack exchange，提问作者Ξένη Γήινος