On rappelle (cf. TD 4) que les polynômes de Bell (ou de Touchard) $b_n(t)$, définis par la série génératrice exponentielle $$ b(x,t)=e^{t(e^x−1)}=\sum_{\ge 0}b_n(t)\frac{x^n}{n!}. $$ Leurs coefficients sont les nombres de Stirling de seconde espèce : $$b_n(t)=\sum_{k=0}^n S(n,k)t^k.$$
Affichez le triangle Stirling2 pour $n\le 9$ (puissances de $t$ en ordre croissant).
Écrivez un fonction calculant $S(n,k)$ au moyen de la récurrence vue en cours $$ S(n,k) = S(n-1,k-1)+kS(n-1,k),$$ avec les conditions initiales $S(0,0)=1$ et $S(n,0)=S(0,n)=0$ (voir ici).
Les nombres de Stirling de première espèce $s(n,k)$ sont les coefficients des polynômes $$t^{\underline n}=(t)_n = t(t-1)\cdots (t-n+1)$$ dont la série génératrice exponentielle est $$ (1+x)^t = \sum_{n\ge 0}{t\choose n} = \sum_{n\ge 0}(t)_n\frac{x^n}{n!}.$$ Affichez le triangle Stirling1 pour $n\le 9$.
Écrivez une fonction calculant $s(n,k)$ au moyen de la récurrence vue en cours $$ s(n,k)= s(n-1,k-1) - (n-1)s(n-1,k),$$ avec les conditions initiales $s(0,0)=1$ et $s(n,0)=s(0,n)=0$ (voir là).
# Modèle avec le triangle de Pascal
from sympy import *
var('x t')
def C(n,k):
if k==n==0: return 1
elif k>n or k<0: return 0
else: return C(n-1,k-1)+C(n-1,k)
for n in range(10):
for k in range(n+1):
print (C(n,k), end=' ')
print ()
# Bell
B = exp(t*(exp(x)-1))
bb = series(B,x,0,10).removeO().as_poly(x)
bb
bb = bb.all_coeffs()[::-1]
bb
bb = [bb[i]*factorial(i) for i in range(len(bb))]
bb
# Stirling 2
def S(n,k):
if k==n==0: return 1
elif k>n or k<0: return 0
else: return S(n-1,k-1)+k*S(n-1,k)
for n in range(10):
for k in range(n+1):
print (S(n,k), end=' ')
print()
G = (1+x)**t
gg = series(G,x,0,10).removeO().as_poly(x)
gg
gg = gg.all_coeffs()[::-1]
gg
gg = [gg[i]*factorial(i) for i in range(len(gg))]
gg
# s(n,k)
for p in gg: print (p.as_poly(t).all_coeffs()[::-1])
# Stirling 1
def s(n,k):
if k==n==0: return 1
elif k>n or k<0: return 0
else: return s(n-1,k-1)-(n-1)*s(n-1,k)
for n in range(10):
for k in range(n+1):
print (s(n,k), end=' ')
print ()
On a vu que l'opération "intégrale discrète" $$\Sigma_a^b f(t)\delta t := \sum_{k=a}^{b-1} f(k)$$ était à l'opérateur de différence finie ce que l'intégrale ordinaire est à la dérivée : $$ \Sigma_a^b f(t)\delta t = F(b)-F(a)\ {\rm si}\ \Delta F(t):=F(t+1)-F(t)=f(t)$$ et aussi que $$\Sigma_0^n t^{\underline n} \delta t = \frac {t^{\underline{n+1}}}{n+1}$$
Sachant que $$t^n = \sum_{k=1}^n S(n,k) t^{\underline k}$$ on a donc $$S_p(n):=\sum_{k=0}^{n-1}k^p = \Sigma_0^n t^p\delta t = \Sigma_0^n \sum_{j=1}^p S(p,j)t^{\underline j}\delta t$$ $$ = \left.\sum_{j=1}^p S(p,j)\frac {t^{\underline{j+1} }}{j+1}\right|_{t=n}$$ et il ne reste plus qu'à redévelopper les $t^{\underline {j+1}}$ au moyen des nombres de Stirling de première espèce (ou plus brutalement, avec la fonction expand).
Utilisez vos nombres de Stirling pour donner l'expression des polynômes $S_p(t)$ pour $p$ de 1 à 8. Présentez les résultats sous forme factorisée.
from functools import reduce
def fallpow(n):
return reduce(lambda u,v:u*v, [t-i for i in range(n)])
fallpow(4)
def sumpow(p):
return factor(sum([S(p,j)*fallpow(j+1)*Rational(1,j+1) for j in range(1,p+1)]))
Le résultat devra ressembler à ceci :
>>> for p in range(1,5): print (sumpow(p)) t*(t - 1)/2 t*(t - 1)*(2*t - 1)/6 t**2*(t - 1)**2/4 t*(t - 1)*(2*t - 1)*(3*t**2 - 3*t - 1)/30
for p in range(1,9): print (sumpow(p))
# Un test
S4 = (t*(t - 1)*(2*t - 1)*(3*t**2 - 3*t - 1)/30).as_poly()
S4
# Un petit test pour voir
for n in range(1,10):
print (sum([k**4 for k in range(n)]), S4.eval(n))
Les polynômes de Bernoulli sont définis par la série génératrice exponentielle (TD 4)
$$B(x,t)=\sum_{n\ge 0}B_n(t)\frac{x^n}{n!} = \frac{xe^{tx}}{e^x-1}$$
D'après le calcul vu en cours, on doit avoir $$S_p(n)=\frac{B_{p+1}(n)-B_{p+1}(0)}{p+1}\ \text{et aussi}\ S(x,n)=\sum_{p\ge 0}S_p(n)\frac{x^p}{p!}=\frac{e^{nx}-1}{e^x-1}.$$ Calculez les premiers polynômes de Bernoulli, recalculez les sommes $S_p(t)$ par cette méthode, et comparez à vos résultats précédents.
Be = x*exp(t*x)/(exp(x)-1)
be = series(Be,x,0,10).removeO().as_poly(x).all_coeffs()[::-1]
be = [be[i]*factorial(i) for i in range(len(be))]
be
def sumpow2(p):
return factor((be[p+1]-be[p+1].as_poly(t).eval(0))*Rational(1,p+1))
for i in range(8): print (sumpow2(i))