Mathématiques pour l'informatique 4

Solution du TD 5

Nombres de Stirling et calcul de sommes

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 ).

In [1]:
from sympy import *
var('x t')
Out[1]:
(x, t)
In [2]:
B = exp(t*(exp(x)-1))
bb = series(B,x,0,10).removeO().as_poly(x)
bb
Out[2]:
Poly((t**9/362880 + t**8/10080 + 11*t**7/8640 + 7*t**6/960 + 331*t**5/17280 + 37*t**4/1728 + 605*t**3/72576 + 17*t**2/24192 + t/362880)*x**9 + (t**8/40320 + t**7/1440 + 19*t**6/2880 + 5*t**5/192 + 27*t**4/640 + 23*t**3/960 + 127*t**2/40320 + t/40320)*x**8 + (t**7/5040 + t**6/240 + t**5/36 + 5*t**4/72 + 43*t**3/720 + t**2/80 + t/5040)*x**7 + (t**6/720 + t**5/48 + 13*t**4/144 + t**3/8 + 31*t**2/720 + t/720)*x**6 + (t**5/120 + t**4/12 + 5*t**3/24 + t**2/8 + t/120)*x**5 + (t**4/24 + t**3/4 + 7*t**2/24 + t/24)*x**4 + (t**3/6 + t**2/2 + t/6)*x**3 + (t**2/2 + t/2)*x**2 + t*x + 1, x, domain='QQ[t]')
In [3]:
bb = bb.all_coeffs()[::-1]
bb
Out[3]:
[1,
 t,
 t**2/2 + t/2,
 t**3/6 + t**2/2 + t/6,
 t**4/24 + t**3/4 + 7*t**2/24 + t/24,
 t**5/120 + t**4/12 + 5*t**3/24 + t**2/8 + t/120,
 t**6/720 + t**5/48 + 13*t**4/144 + t**3/8 + 31*t**2/720 + t/720,
 t**7/5040 + t**6/240 + t**5/36 + 5*t**4/72 + 43*t**3/720 + t**2/80 + t/5040,
 t**8/40320 + t**7/1440 + 19*t**6/2880 + 5*t**5/192 + 27*t**4/640 + 23*t**3/960 + 127*t**2/40320 + t/40320,
 t**9/362880 + t**8/10080 + 11*t**7/8640 + 7*t**6/960 + 331*t**5/17280 + 37*t**4/1728 + 605*t**3/72576 + 17*t**2/24192 + t/362880]
In [4]:
bb = [bb[i]*factorial(i) for i in range(len(bb))]
bb
Out[4]:
[1,
 t,
 t**2 + t,
 t**3 + 3*t**2 + t,
 t**4 + 6*t**3 + 7*t**2 + t,
 t**5 + 10*t**4 + 25*t**3 + 15*t**2 + t,
 t**6 + 15*t**5 + 65*t**4 + 90*t**3 + 31*t**2 + t,
 t**7 + 21*t**6 + 140*t**5 + 350*t**4 + 301*t**3 + 63*t**2 + t,
 t**8 + 28*t**7 + 266*t**6 + 1050*t**5 + 1701*t**4 + 966*t**3 + 127*t**2 + t,
 t**9 + 36*t**8 + 462*t**7 + 2646*t**6 + 6951*t**5 + 7770*t**4 + 3025*t**3 + 255*t**2 + t]
In [5]:
for p in bb: print p.as_poly(t).all_coeffs()[::-1]
[1]
[0, 1]
[0, 1, 1]
[0, 1, 3, 1]
[0, 1, 7, 6, 1]
[0, 1, 15, 25, 10, 1]
[0, 1, 31, 90, 65, 15, 1]
[0, 1, 63, 301, 350, 140, 21, 1]
[0, 1, 127, 966, 1701, 1050, 266, 28, 1]
[0, 1, 255, 3025, 7770, 6951, 2646, 462, 36, 1]

In [6]:
# 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)
In [7]:
for n in  range(10):
    for k in range(n+1): 
        print S(n,k),
    print 
        
1
0 1
0 1 1
0 1 3 1
0 1 7 6 1
0 1 15 25 10 1
0 1 31 90 65 15 1
0 1 63 301 350 140 21 1
0 1 127 966 1701 1050 266 28 1
0 1 255 3025 7770 6951 2646 462 36 1

In [8]:
G = (1+x)**t
gg = series(G,x,0,10).removeO().as_poly(x)
gg
Out[8]:
Poly((t**9/362880 - t**8/10080 + 13*t**7/8640 - t**6/80 + 1069*t**5/17280 - 89*t**4/480 + 29531*t**3/90720 - 761*t**2/2520 + t/9)*x**9 + (t**8/40320 - t**7/1440 + 23*t**6/2880 - 7*t**5/144 + 967*t**4/5760 - 469*t**3/1440 + 363*t**2/1120 - t/8)*x**8 + (t**7/5040 - t**6/240 + 5*t**5/144 - 7*t**4/48 + 29*t**3/90 - 7*t**2/20 + t/7)*x**7 + (t**6/720 - t**5/48 + 17*t**4/144 - 5*t**3/16 + 137*t**2/360 - t/6)*x**6 + (t**5/120 - t**4/12 + 7*t**3/24 - 5*t**2/12 + t/5)*x**5 + (t**4/24 - t**3/4 + 11*t**2/24 - t/4)*x**4 + (t**3/6 - t**2/2 + t/3)*x**3 + (t**2/2 - t/2)*x**2 + t*x + 1, x, domain='QQ[t]')
In [9]:
gg = gg.all_coeffs()[::-1]
gg
Out[9]:
[1,
 t,
 t**2/2 - t/2,
 t**3/6 - t**2/2 + t/3,
 t**4/24 - t**3/4 + 11*t**2/24 - t/4,
 t**5/120 - t**4/12 + 7*t**3/24 - 5*t**2/12 + t/5,
 t**6/720 - t**5/48 + 17*t**4/144 - 5*t**3/16 + 137*t**2/360 - t/6,
 t**7/5040 - t**6/240 + 5*t**5/144 - 7*t**4/48 + 29*t**3/90 - 7*t**2/20 + t/7,
 t**8/40320 - t**7/1440 + 23*t**6/2880 - 7*t**5/144 + 967*t**4/5760 - 469*t**3/1440 + 363*t**2/1120 - t/8,
 t**9/362880 - t**8/10080 + 13*t**7/8640 - t**6/80 + 1069*t**5/17280 - 89*t**4/480 + 29531*t**3/90720 - 761*t**2/2520 + t/9]
In [10]:
gg = [gg[i]*factorial(i) for i in range(len(gg))]
gg
Out[10]:
[1,
 t,
 t**2 - t,
 t**3 - 3*t**2 + 2*t,
 t**4 - 6*t**3 + 11*t**2 - 6*t,
 t**5 - 10*t**4 + 35*t**3 - 50*t**2 + 24*t,
 t**6 - 15*t**5 + 85*t**4 - 225*t**3 + 274*t**2 - 120*t,
 t**7 - 21*t**6 + 175*t**5 - 735*t**4 + 1624*t**3 - 1764*t**2 + 720*t,
 t**8 - 28*t**7 + 322*t**6 - 1960*t**5 + 6769*t**4 - 13132*t**3 + 13068*t**2 - 5040*t,
 t**9 - 36*t**8 + 546*t**7 - 4536*t**6 + 22449*t**5 - 67284*t**4 + 118124*t**3 - 109584*t**2 + 40320*t]
In [11]:
for p in gg: print p.as_poly(t).all_coeffs()[::-1]
[1]
[0, 1]
[0, -1, 1]
[0, 2, -3, 1]
[0, -6, 11, -6, 1]
[0, 24, -50, 35, -10, 1]
[0, -120, 274, -225, 85, -15, 1]
[0, 720, -1764, 1624, -735, 175, -21, 1]
[0, -5040, 13068, -13132, 6769, -1960, 322, -28, 1]
[0, 40320, -109584, 118124, -67284, 22449, -4536, 546, -36, 1]

In [12]:
# 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)
In [13]:
for n in  range(10):
    for k in range(n+1): 
        print s(n,k),
    print 
        
1
0 1
0 -1 1
0 2 -3 1
0 -6 11 -6 1
0 24 -50 35 -10 1
0 -120 274 -225 85 -15 1
0 720 -1764 1624 -735 175 -21 1
0 -5040 13068 -13132 6769 -1960 322 -28 1
0 40320 -109584 118124 -67284 22449 -4536 546 -36 1

Calcul de sommes

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\] \[ = \sum_{j=1}^p S(p,j)\left[\frac {t^{\underline{j+1} }}{j+1}\right]_0^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 sommes \(S_p(n)\) pour \(p\) de 1 à 8. Présentez les résultats sous forme factorisée.

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
In [14]:
def fallpow(n):
    return reduce(lambda u,v:u*v, [t-i for i in range(n)])

fallpow(4)
Out[14]:
t*(t - 3)*(t - 2)*(t - 1)
In [15]:
def sumpow(p):
    return factor(sum([S(p,j)*fallpow(j+1)*Rational(1,j+1) for j in range(1,p+1)]))
In [16]:
for p in range(1,9): 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
t**2*(t - 1)**2*(2*t**2 - 2*t - 1)/12
t*(t - 1)*(2*t - 1)*(3*t**4 - 6*t**3 + 3*t + 1)/42
t**2*(t - 1)**2*(3*t**4 - 6*t**3 - t**2 + 4*t + 2)/24
t*(t - 1)*(2*t - 1)*(5*t**6 - 15*t**5 + 5*t**4 + 15*t**3 - t**2 - 9*t - 3)/90

Sommes de puissances et polynômes de Bernoulli

Les polynômes de Bernoulli sont définis par la série génératreice exponentielle (TP 3)

\[B(x,t)=\sum_{n\ge 0}B_n(t)\frac{x^n}{n!} = \frac{xe^{tx}}{e^x-1}\]

Calculez les premiers polynômes \(B_n(t)\).

Les nombres de Bernoulli sont les termes constants \(B_n(0)\). Ils ne sont pas entiers (ni positifs). Affichez les 20 premiers, et consultez la page qui leur est consacrée sur Wikipedia.

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.

In [17]:
Be = x*exp(t*x)/(exp(x)-1)
be = series(Be,x,0,10).removeO().as_poly(x).all_coeffs()[::-1]
In [18]:
be
Out[18]:
[1,
 t - 1/2,
 t**2/2 - t/2 + 1/12,
 t**3/6 - t**2/4 + t/12,
 t**4/24 - t**3/12 + t**2/24 - 1/720,
 t**5/120 - t**4/48 + t**3/72 - t/720,
 t**6/720 - t**5/240 + t**4/288 - t**2/1440 + 1/30240,
 t**7/5040 - t**6/1440 + t**5/1440 - t**3/4320 + t/30240,
 t**8/40320 - t**7/10080 + t**6/8640 - t**4/17280 + t**2/60480 - 1/1209600,
 t**9/362880 - t**8/80640 + t**7/60480 - t**5/86400 + t**3/181440 - t/1209600]
In [19]:
be = [be[i]*factorial(i) for i in range(len(be))]
be
                    
Out[19]:
[1,
 t - 1/2,
 t**2 - t + 1/6,
 t**3 - 3*t**2/2 + t/2,
 t**4 - 2*t**3 + t**2 - 1/30,
 t**5 - 5*t**4/2 + 5*t**3/3 - t/6,
 t**6 - 3*t**5 + 5*t**4/2 - t**2/2 + 1/42,
 t**7 - 7*t**6/2 + 7*t**5/2 - 7*t**3/6 + t/6,
 t**8 - 4*t**7 + 14*t**6/3 - 7*t**4/3 + 2*t**2/3 - 1/30,
 t**9 - 9*t**8/2 + 6*t**7 - 21*t**5/5 + 2*t**3 - 3*t/10]
In [20]:
def sumpow2(p):
    return factor((be[p+1]-be[p+1].as_poly(t).eval(0))*Rational(1,p+1))
In [21]:
for i in range(8): print sumpow2(i)
t
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
t**2*(t - 1)**2*(2*t**2 - 2*t - 1)/12
t*(t - 1)*(2*t - 1)*(3*t**4 - 6*t**3 + 3*t + 1)/42
t**2*(t - 1)**2*(3*t**4 - 6*t**3 - t**2 + 4*t + 2)/24

In []: