如何使用 sympy 找到生成函数的第 n 项?
How to find the nth term of a generating function using sympy?
我有一个有理函数:f(x) = P(x)/Q(x)
。
例如:
f(x) = (5x + 3)/(1-x^2)
因为f(x)是生成函数所以可以写成:
f(x) = a0 + a1*x + a2*x² + ... + a_n*x^n + ... = P(x)/Q(x)
如何使用 sympy 找到生成函数 f(x)
(即 a_n
)的第 nth 项?
如果 Sympy 中没有这样的实现,我也很想知道这是否在其他包中实现,例如 Maxima。
感谢任何帮助。
采用的方法,您可以尝试以下方法:
from sympy import *
from sympy.abc import x
f = (5*x + 3) / (1-x**2)
print(f.series(n=20))
k = 50
coeff50 = Poly(series(f, x, n=k + 1).removeO(), x).coeff_monomial(x ** k)
print(f"The coeffcient of x^{k} of the generating function of {f} is {coeff50}")
# to get the first 100 coeffcients (reversing the list to get a[0] the
# coefficient of x**0 etc.):
a = Poly(series(f, x, n=100).removeO(), x).all_coeffs()[::-1]
输出:
3 + 5*x + 3*x**2 + 5*x**3 + 3*x**4 + 5*x**5 + 3*x**6 + 5*x**7 + 3*x**8 + 5*x**9 + 3*x**10 + 5*x**11 + 3*x**12 + 5*x**13 + 3*x**14 + 5*x**15 + 3*x**16 + 5*x**17 + 3*x**18 + 5*x**19 + O(x**20)
The coeffcient of x^50 of the generating function of (5*x + 3)/(1 - x**2) is 3
按照 Cut The Knot 的 this 示例,该方法可用于找出可以使用 1、5、10、25 和50 美分。
f = 1/((1 - x)*(1 - x**5)*(1 - x**10)*(1 - x**25)*(1 - x**50))
a = Poly(series(f, x, n=101).removeO(), x).all_coeffs()[::-1]
print(a[50]) # there are 50 ways to pay 50 cents
print(a[100]) # there are 292 ways to pay 100 cents
你可以求第 k 个导数,用 0 代替 x
然后除以 factorial(k)
:
>>> f = (5*x + 3) / (1-x**2)
>>> f.diff(x, 20).subs(x, 0)/factorial(20)
3
参考文献here讨论有理生成函数。寻找重复,您可以使用微分很快看到模式:
[f.diff(x,i).subs(x,0)/factorial(i) for i in range(6)]
[3, 5, 3, 5, 3, 5]
最大值:
powerseries((5*x+3)/(1-x^2),x,0);
returns
使用part
提取生成器:
part(''%,1);
(4-(-1)^i1)x^i1
和coeff
得到系数:
a(i1) := coeff(''%, x, i1);
[a(0), a(1), a(2)];
[3, 5, 3]
要获得有理式母函数a_n
的通式,可以使用SymPy的rational_algorithm
。
例如:
from sympy import simplify
from sympy.abc import x, n
from sympy.series.formal import rational_algorithm
f = (5*x + 3)/(1-x**2)
func_n, independent_term, order = rational_algorithm(f, x, n, full=True)
print(f"The general formula for a_n is {func_n}")
for k in range(10):
print(f"a_{k} = {simplify(func_n.subs(n, k))}")
输出:
The general formula for a_n is (-1)**(-n - 1) + 4
a_0 = 3
a_1 = 5
a_2 = 3
a_3 = 5
a_4 = 3
a_5 = 5
a_6 = 3
a_7 = 5
a_8 = 3
a_9 = 5
我有一个有理函数:f(x) = P(x)/Q(x)
。
例如:
f(x) = (5x + 3)/(1-x^2)
因为f(x)是生成函数所以可以写成:
f(x) = a0 + a1*x + a2*x² + ... + a_n*x^n + ... = P(x)/Q(x)
如何使用 sympy 找到生成函数 f(x)
(即 a_n
)的第 nth 项?
如果 Sympy 中没有这样的实现,我也很想知道这是否在其他包中实现,例如 Maxima。
感谢任何帮助。
采用
from sympy import *
from sympy.abc import x
f = (5*x + 3) / (1-x**2)
print(f.series(n=20))
k = 50
coeff50 = Poly(series(f, x, n=k + 1).removeO(), x).coeff_monomial(x ** k)
print(f"The coeffcient of x^{k} of the generating function of {f} is {coeff50}")
# to get the first 100 coeffcients (reversing the list to get a[0] the
# coefficient of x**0 etc.):
a = Poly(series(f, x, n=100).removeO(), x).all_coeffs()[::-1]
输出:
3 + 5*x + 3*x**2 + 5*x**3 + 3*x**4 + 5*x**5 + 3*x**6 + 5*x**7 + 3*x**8 + 5*x**9 + 3*x**10 + 5*x**11 + 3*x**12 + 5*x**13 + 3*x**14 + 5*x**15 + 3*x**16 + 5*x**17 + 3*x**18 + 5*x**19 + O(x**20)
The coeffcient of x^50 of the generating function of (5*x + 3)/(1 - x**2) is 3
按照 Cut The Knot 的 this 示例,该方法可用于找出可以使用 1、5、10、25 和50 美分。
f = 1/((1 - x)*(1 - x**5)*(1 - x**10)*(1 - x**25)*(1 - x**50))
a = Poly(series(f, x, n=101).removeO(), x).all_coeffs()[::-1]
print(a[50]) # there are 50 ways to pay 50 cents
print(a[100]) # there are 292 ways to pay 100 cents
你可以求第 k 个导数,用 0 代替 x
然后除以 factorial(k)
:
>>> f = (5*x + 3) / (1-x**2)
>>> f.diff(x, 20).subs(x, 0)/factorial(20)
3
参考文献here讨论有理生成函数。寻找重复,您可以使用微分很快看到模式:
[f.diff(x,i).subs(x,0)/factorial(i) for i in range(6)]
[3, 5, 3, 5, 3, 5]
最大值:
powerseries((5*x+3)/(1-x^2),x,0);
returns
使用part
提取生成器:
part(''%,1);
(4-(-1)^i1)x^i1
和coeff
得到系数:
a(i1) := coeff(''%, x, i1);
[a(0), a(1), a(2)];
[3, 5, 3]
要获得有理式母函数a_n
的通式,可以使用SymPy的rational_algorithm
。
例如:
from sympy import simplify
from sympy.abc import x, n
from sympy.series.formal import rational_algorithm
f = (5*x + 3)/(1-x**2)
func_n, independent_term, order = rational_algorithm(f, x, n, full=True)
print(f"The general formula for a_n is {func_n}")
for k in range(10):
print(f"a_{k} = {simplify(func_n.subs(n, k))}")
输出:
The general formula for a_n is (-1)**(-n - 1) + 4
a_0 = 3
a_1 = 5
a_2 = 3
a_3 = 5
a_4 = 3
a_5 = 5
a_6 = 3
a_7 = 5
a_8 = 3
a_9 = 5