MPA Mathématiques en Python : CC du 5 mai 2021

Exercice 1 : Matplotlib (tracé de polygones hyperboliques)

Soit $\alpha > 0$ un paramètre. On considère la courbe paramétrée définie par :

$$ x(t) = (1 - \alpha) \cos(t) + \alpha \cos\left(\frac{(1-\alpha)}{\alpha} t\right) $$

$$ y(t) = (1 - \alpha) \sin(t) - \alpha \sin\left(\frac{(1-\alpha)}{\alpha} t\right) $$

pour $t\in\mathbb{R}$.

Question 1.

Pour $\alpha= \frac 1 5$, tracer la courbe, en prenant $0 \le t \le 2\pi$. Puis, remplacer $\frac 1 n$ par d'autres valeurs de la forme $\frac 1 n$ avec $n$ entier.

Correction.

In [11]:
from matplotlib import pyplot as plt
import numpy as np
import math

Faisons une fonction (ce n'est pas obligatoire!).

In [34]:
def courbe(alpha, nb_points, t_min, t_max):
    """
    Trace la courbe avec matplotlib, pour t entre t_min et t_max,
    en prenant 'nb_points' sur la courbe.    
    """
    
    t= np.linspace(t_min, t_max, nb_points) 
    x= (1-alpha)*np.cos(t) + alpha*np.cos( ((1-alpha)/alpha) *t)
    y= (1-alpha)*np.sin(t) - alpha*np.sin( ((1-alpha)/alpha) *t)
    plt.plot(x,y)
In [41]:
courbe(1/5, 100, 0, 2*math.pi)
In [42]:
courbe(1/7, 100, 0, 2*math.pi)

Question 2.

Pour $\alpha = \frac{1}{\sqrt{2}}$, tracer la courbe en prenant $0 \le t \le R$, où on augmentera $R$ progressivement. Commenter.

In [43]:
courbe(1/math.sqrt(2), 100, 0, 2*math.pi)
In [46]:
courbe(1/math.sqrt(2), 100, 0, 20*math.pi)
In [49]:
courbe(1/math.sqrt(2), 100, 0, 100*math.pi)

A ce stade, il est bon de remarquer que sur un intervalle aussi grand que $[0, 100\pi]$, il est probablement insuffisant de ne prendre que 100 points! essayons :

In [50]:
courbe(1/math.sqrt(2), 1000, 0, 100*math.pi)

C'est effectivement très différent! On se dit qu'on va augmenter le nombre de points avec la longueur de l'intervalle. Par exemple :

In [51]:
test= lambda N : courbe(1/math.sqrt(2), 10*N, 0, N*math.pi)
In [52]:
test(100)
In [59]:
test(250)
In [60]:
test(1000)

Exercice 2 : une méthode de Monte-Carlo

Dans cet exercice, nous aurons besoin de la fonction de Numpy pour produire des tableaux de nombres tirés au sort entre 0 et 1 :

In [68]:
import numpy as np
In [73]:
A= np.random.rand(2,6)
print(A)

[[0.0866418  0.43705133 0.39226443 0.46559073 0.49904447 0.43928325]
 [0.88586268 0.91314836 0.25214048 0.85867064 0.95048529 0.82791558]]

Ici la matrice tirée au sort comporte 2 lignes et 6 colonnes. On peut y penser comme à une suite de 6 vecteurs tirés au hasard, avec A[0,i] et A[1,i] qui sont les coordonnées du $i$-ième vecteur.

Question 1.

Ecrire une fonction compte(N) qui fait les choses suivantes :

  • on tire une matrice au hasard, comme ci-dessus, de 2 lignes et $N$ colonnes;
  • on compte combien d'indices $i$ vérifients A[0,i]**2 + A[1,i]**2 < 1, et on renvoie ce nombre.

Bonus pour ceux qui arrivent à écrire une fonction rapide.

Correction.

In [79]:
def compte(N):
    A= np.random.rand(2,N)
    # pour le plaisir, voici une version entièrement "en numpy", qui est très rapide
    # (ci-dessous une version plus simple)
    x= A[0,:] # première ligne
    y= A[1,:] # deuxième ligne
    test= x**2 + y**2 < 1
    return np.sum(test)
In [97]:
# ou alors, beaucoup plus simplement (et c'est ce que j'attends)
def compte2(N):
    A= np.random.rand(2,N)
    total= 0
    for i in range(N):
        if A[0,i]**2 + A[1,i]**2 < 1 :
            total= total + 1
    return total
In [95]:
%time compte2(1000000)

CPU times: user 1.15 s, sys: 6.55 ms, total: 1.15 s
Wall time: 1.17 s
Out[95]:
785715
In [96]:
%time compte(1000000)

CPU times: user 43.3 ms, sys: 4.68 ms, total: 48 ms
Wall time: 45.3 ms
Out[96]:
785576

Ah! la première version est beaucoup plus rapide.

Question 2.

Observez la valeur de 4*compte(N)/N, pour $N$ de plus en plus grand. Puis (et c'est beaucoup plus difficile), expliquez, même sans être trop rigoureux.

In [98]:
essai= lambda N : 4*compte(N)/N
In [99]:
essai(10)
Out[99]:
2.8
In [100]:
essai(100)
Out[100]:
3.4
In [101]:
essai(1000)
Out[101]:
3.24
In [102]:
essai(1000000)
Out[102]:
3.141608

Ça m'a tout l'air de converger vers $\pi$. Pourquoi?

Chaque nombre tiré au hasard est entre 0 et 1, donc chaque vecteur tiré au hasard est un point dans le carré $[0,1] \times [0,1]$. Sans entrer dans le détail, on imagine facilement que la probabilité que le point tiré au hasard tombe dans une zone donnée $Z$ est donné par l'aire de $Z$ (notez que l'aire du carré est 1).

En comptant le nombre de points moyen dans une zone $Z$, on devrait donc trouver une approximation de l'aire (c'est la loi des grands nombres, en fait, que vous ne connaissez pas bien).

Dans l'exercice on a compté les points dans le quart de disque de rayon 1 et de centre $(0,0)$. Son aire est $\frac \pi 4$. D'où les valeurs observées.

Et oui, cette méthode pour calculer les aires s'appelle méthode de Monte-Carlo, parce qu'elle s'appuie sur le hasard (qui est également présent au casino).

Exercice 3 (polynômes de Hermite)

Prenons une variable formelle $X$ grâce à sympy (vous devez faire la même chose dans votre notebook):

In [1]:
import sympy as sp
In [2]:
sp.init_printing()
In [3]:
X= sp.Symbol("X")

Les polynômes de Hermite $(H_n)_{n\ge 0}$ sont définis par récurrence de la manière suivante :

  • $H_0= 1$, $H_1 = X$,
  • $H_{n} = X \cdot H_{n-1} -(n-1) H_{n-2}$, pour $n\ge 2$.

Question 1.

Ecrire une fonction H(n) qui renvoie le polynôme $H_n$, comme une expression formelle sympy développée.

Note: il est normal que ça soit lent pour des valeurs de $n$ supérieures à 25. On ne vous demande pas d'optimiser tout ça.

Correction.

In [4]:
def H(n):
    if n == 0:
        return 1
    elif n == 1:
        return X
    else:
        return (X * H(n-1) - (n-1)*H(n-2)).expand()
In [5]:
for i in range(10):
    display(H(i))
$$1$$
$$X$$
$$X^{2} - 1$$
$$X^{3} - 3 X$$
$$X^{4} - 6 X^{2} + 3$$
$$X^{5} - 10 X^{3} + 15 X$$
$$X^{6} - 15 X^{4} + 45 X^{2} - 15$$
$$X^{7} - 21 X^{5} + 105 X^{3} - 105 X$$
$$X^{8} - 28 X^{6} + 210 X^{4} - 420 X^{2} + 105$$
$$X^{9} - 36 X^{7} + 378 X^{5} - 1260 X^{3} + 945 X$$
In [113]:
H(2)
Out[113]:
$$X^{2} - 1$$
In [114]:
H(3)
Out[114]:
$$X^{3} - 3 X$$
In [6]:
H(25)
Out[6]:
$$X^{25} - 300 X^{23} + 37950 X^{21} - 2656500 X^{19} + 113565375 X^{17} - 3088978200 X^{15} + 54057118500 X^{13} - 602350749000 X^{11} + 4141161399375 X^{9} - 16564645597500 X^{7} + 34785755754750 X^{5} - 31623414322500 X^{3} + 7905853580625 X$$

Question 2.

En regardant des exemples, devinez une relation entre $H_n'$ (la dérivée de $H_n$) et $H_{n-1}$. On ne vous demande pas de la démontrer.

Correction.

In [119]:
N= 10
for n in range(1,N):
    print("dérivée de H_", n,"=", sp.diff(H(n), X), "et H_", n-1, "=", H(n-1))

dérivée de H_ 1 = 1 et H_ 0 = 1
dérivée de H_ 2 = 2*X et H_ 1 = X
dérivée de H_ 3 = 3*X**2 - 3 et H_ 2 = X**2 - 1
dérivée de H_ 4 = 4*X**3 - 12*X et H_ 3 = X**3 - 3*X
dérivée de H_ 5 = 5*X**4 - 30*X**2 + 15 et H_ 4 = X**4 - 6*X**2 + 3
dérivée de H_ 6 = 6*X**5 - 60*X**3 + 90*X et H_ 5 = X**5 - 10*X**3 + 15*X
dérivée de H_ 7 = 7*X**6 - 105*X**4 + 315*X**2 - 105 et H_ 6 = X**6 - 15*X**4 + 45*X**2 - 15
dérivée de H_ 8 = 8*X**7 - 168*X**5 + 840*X**3 - 840*X et H_ 7 = X**7 - 21*X**5 + 105*X**3 - 105*X
dérivée de H_ 9 = 9*X**8 - 252*X**6 + 1890*X**4 - 3780*X**2 + 945 et H_ 8 = X**8 - 28*X**6 + 210*X**4 - 420*X**2 + 105

Ma foi, il semble que $H'_n = n H_{n-1}$. Essayons :

In [127]:
N= 30
for n in range(1,N):
    if sp.diff(H(n), X) != n*H(n-1):
        print("ça ne marche pas pour n=", n)
    else:
        print(".", end= "")
print("Toutes les valeurs jusqu'à n=", N, "ont été testées")

.............................Toutes les valeurs jusqu'à n= 30 ont été testées