如何用程序推导泰勒级数？
所以我知道，在物理学里，某一时刻的位置可以用它的导数来计算，像这样：
现在我发现这看起来和泰勒级数非常相似：
事实上它们完全相同。这也意味着我可以用这种方法来计算超越函数（三角函数、反三角函数、双曲函数、反双曲函数），其实我已经在用泰勒级数通过整数运算来计算所有这些函数了，但我想尝试自己通过导数来程序化推导泰勒级数，看看会走向何方。
我们需要求出导数，导数其实很容易求。
例如，给定一个多项式F(x) = ax6 + bx5 + cx4 + dx3 + ex2 + fx + g，它的导数是F'(x) = 6ax5 + 5bx4 + 4cx3 + 3dx2 + 2ex + f。换句话说，将多项式表示为系数和指数对的列表时，导数就是把每个指数乘以系数得到新的系数，指数减1作为新的指数；当指数为0时，该项将被忽略。
而给定一个以分式形式表示的多项式组合f(x)/g(x)，其导数是 (f'(x) * g(x) - g'(x) * f(x)) / g(x)2
要把两个多项式相减，只需将两多项式中具有相同指数的项配对，然后相减系数，对第二个在第一项中未出现的指数用负号代入。
要把两个多项式相乘，就是把一个多项式的每一项都与另一个多项式的每一项相乘，然后求和；要得到系数-次数对的乘积，只需把系数相乘、指数相加；要把项相加，只需按指数分组、对相同指数的系数求和。
所以我选择先求出反正切的高阶导数，看看能否不借助库就推导出它的泰勒级数。反正切的导数是1/(x2 + 1)。
于是我尝试自己去做，结果是正确的，但计算量很快变得非常大：
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}))")
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: 3365856,
12: 906192,
52: 906192,
10: 201376,
54: 201376,
8: 35960,
56: 35960,
6: 4960,
58: 4960,
4: 496,
60: 496,
2: 32,
62: 32,
0: 1,
64: 1}))]
正如你所见，项的数量在求导次数增加时呈组合式爆炸，这意味着随着阶数的提升，计算量越来越大，阶数越高，项数越多，因此我进行的高阶导数计算成本也会越高，形成了一个正反馈循环。
下面的计算花了很长时间，我不知道花了多久，因为我在中途手动中断了输入：
poly_pair = (Counter({0: 1}), Counter({0: 1, 2:1}))
derivatives = [poly_pair, *[(poly_pair := frac_poly_dif(poly_pair)) for _ in range(12)]]
而且我至今还没能计算出这个：
poly_pair = (Counter({0: 1}), Counter({0: 1, 2:1}))
derivatives = [poly_pair, *[(poly_pair := frac_poly_dif(poly_pair)) for _ in range(16)]]
但我已经用这个网站验证了我的结果，它们是一致的：
In [5]: for poly in derivatives:
...: print(sympy_simplify(poly))
...:
1/(x**2 + 1)
-2*x/(x**4 + 2*x**2 + 1)
2*(3*x**2 - 1)/(x**6 + 3*x**4 + 3*x**2 + 1)
24*x*(1 - x**2)/(x**8 + 4*x**6 + 6*x**4 + 4*x**2 + 1)
24*(5*x**4 - 10*x**2 + 1)/(x**10 + 5*x**8 + 10*x**6 + 10*x**4 + 5*x**2 + 1)
(-720*x**5 + 2400*x**3 - 720*x)/(x**12 + 6*x**10 + 15*x**8 + 20*x**6 + 15*x**4 + 6*x**2 + 1)
我也不知道该如何化简多项式、加速计算，怎么才能在不使用sympy等库的情况下计算更高阶（16阶及以上）的导数呢？
解决方案
是的，泰勒级数通常计算量很大。如果我们只针对特定的一组函数，有时可以找到一个不容易爆炸的表示方式。
如果我们一般性地讨论有理函数，那么P_num / P_den的导数确实是dP_num / P_den - P_num dP_den / P_den^2 = (dP_num P_den - P_num dP_den) / P_den^2，次数会指数级增长。然而，这种把新得到的导数视为完全一般的函数来处理的方式有点多此一举——我们可以利用每次导数时分母都是P_den的幂的这一事实，因此设第k 次导数为P[k] / P_den^(k+1)，其中P[k] = P_num。
d(P[k] / P_den^(k+1)) = dP[k] / P_den^(k+1) - (k+1) P[k] dP_den / P_den^(k+2)
P[k+1] = dP[k] P_den - (k+1) P[k] dP_den
这里P[k] 的次数随k 线性增加，在每一步至多为deg(P_den) - 1。
def frac_poly_dif(poly_pair: tuple[Counter, Counter]) -> tuple[Counter, Counter]:
num_poly, (den_poly, k) = poly_pair
return poly_sub(
poly_mul(poly_dif(num_poly), den_poly), poly_mul(Counter({0: k+1}), poly_mul(poly_dif(den_poly), num_poly))
), (den_poly, k+1)
poly_pair = (Counter({0: 1}), (Counter({0: 1, 2:1}), 0))
derivatives = [poly_pair, *[(poly_pair := frac_poly_dif(poly_pair)) for _ in range(256)]]
（如上面的公式所示，我把分母从一个多项式改为一个二元组（多项式、k））
相关文章
