Revenir
Revenir

Approximation de Pi

En Python, il est possible de créer une liste de taille n constituée uniquement de 0 en utilisant la...

Sommaire

Activité CAPYTALE : approximation de Pi
Quelques rappels pour commencerGénérer et modifier une liste de taille donnéeGénérer une liste contenant l'ensemble des valeurs d'une suiteGénérer la liste des termes d'une suite
Méthode d'ArchimèdePrincipe de la méthode d'ArchimèdeMéthode d'Archimède - Construction de deux suitesMéthode d'Archimède - Représentation graphique des suites a et bMéthode d'Archimède - Aspect théorique
Méthode de l'arctangenteSuite moyenne
Autres méthodesLe problème de BâleAlgorithme de Chudnovsky

Activité CAPYTALE : approximation de Pi


Quelques rappels pour commencer

Générer et modifier une liste de taille donnée

En Python, il est possible de créer une liste de taille n constituée uniquement de 0 en utilisant la syntaxe suivante
L = [0] * n
On peut alors modifier les éléments individuellement en accédant à ces éléments par leur indice. Attention, on commence à compter les indices à partir de l'indice 0 !
L = [0] * 6 # crée une liste de taille 6 constituée de zéros
L[0] = 3 # modifie le premier élément et lui affecte la valeur 3
L[4] = -24 # modifie le cinquième élément et lui affecte la valeur -24
print(L)
>>> [3, 0, 0, 0, -24, 0]
Un parcours de liste à l'aide d'une boucle for permet alors de modifier l'ensemble des éléments de la liste.
L = [0] * 6
for i in range(len(L)) :
    L[i] = i ** 2
print(L)
>>> [0, 1, 4, 9, 16, 25]

Générer une liste contenant l'ensemble des valeurs d'une suite

On considère une suite
(u_n)
et un entier naturel n.
Pour stocker l'ensemble des valeurs prises par cette suite jusqu'au rang n, on peut alors utiliser une liste de taille (n+1).
L'élément d'indice i de la liste sera le terme de rang i de la suite.
def u(n) :
    return (n+3)/(2*n + 5)
def termes_u(n):"""renvoie une liste contenant tous les termes de la suite (u_n) jusqu'à celui de rang n"""
    L = [0] * (n+1)
    for i in range(len(L)):
        L[i] = u(i)
    return L
termes_u(10)
>>> [0.6, 0.5714285714285714, 0.5555555555555556, 0.5454545454545454, 0.5384615384615384, 0.5333333333333333, 0.5294117647058824, 0.5263157894736842, 0.5238095238095238, 0.5217391304347826, 0.52]
Le même procédé peut être utilisé pour la génération d'une suite définie par récurrence. Il faut pour cela commencer par stocker la valeur de 
u_0
au début de la liste.
La fonction suivante permet de générer la liste des termes de la suite définie par
u_0=3
et, pour tout entier naturel
n
, 
u_{n+1}=2u_n+n+1
.
def termes_u(n):
    L = [0] * n
    L[0] = 3     # On place la valeur de u_0 en début de liste
    for i in range(1, len(L)) :   # on commence à remplir la liste à partir de l'indice 1
        L[i]= 2 * L[i-1] + (i-1) + 1   # Attention ! on exprime u_i en fonction de u_(i-1)
    return L
termes_u(10)

Générer la liste des termes d'une suite

Exercice 1
Créer une fonction qui renvoie la liste des termes de la suite
(u_n)
définie par
u_0=7
 et pour tout entier naturel
n
, 
u_{n+1} = 3u_n -4n+ 1
.
Exercice 2
Compléter la fonction suivante qui renvoie la liste des termes des suites
(u_n)
et 
(v_n)
définie par
u_0=5
,
v_0=3
et, pour tout entier naturel
n
,
u_{n+1}=u_n+v_n
et 
v_{n+1} = u_{n+1}v_n
.
def termes_uv(n) :
    # Liste des termes de la suite (u_n)
    Lu = [0] * (n+1)
    Lu[0] = ...
    # Liste des termes de la suite (v_n)
    Lv = [0] * (n+1)
    Lv[0] = ...
    # Calcul des termes
    for i in range(1, len(Lu)):
        Lu[i] = ...
        Lv[i] = ...
    return Lu, Lv
u, v = termes_uv(5)
print(u)
print(v)

Méthode d'Archimède

Principe de la méthode d'Archimède

La méthode d'Archimède permet d'obtenir une approximation du nombre
\pi
. Pour cela, on considère un cercle de rayon 
1/2
: ce cercle a un périmètre de
\pi
. On construit alors des polygones réguliers inscrits et circonscrits à ce cercle : plus le nombre de côtés de ces polygones sera grand, plus ces polygones seront proches du cercle et, donc, plus leur périmètre sera proche de 
\pi
.

Méthode d'Archimède - Construction de deux suites

Pour tout entier naturel
n⩾2n\geqslant 2n⩾2
, on note
ana_nan​
le périmètre du polygone à
2n2^n2n
côtés inscrit dans le cercle et
bnb_nbn​
le périmètre du polygone à
2n2^n2n
côtés circonscrit au cercle.
Exercice 1
Justifier que l'on a 
a2=22a_2=2\sqrt{2}a2​=22​
et 
b2=4b_2=4b2​=4
.
Exercice 2
On admet que, pour tout entier naturel
n⩾2n\geqslant 2n⩾2
, on a
bn+1=2anbnan+bnb_{n+1}=\dfrac{2a_nb_n}{a_n+b_n}bn+1​=an​+bn​2an​bn​​
 et
an+1=anbn+1a_{n+1}=\sqrt{a_nb_{n+1}}an+1​=an​bn+1​​
.
Calculer 
b3b_3b3​
, 
a3a_3a3​
, 
b4b_4b4​
 et 
a4a_4a4​
.
Exercice 3
Compléter le programme suivant qui renvoie les 
kkk
premiers termes des suites
(an)(a_n)(an​)
et
(bn)(b_n)(bn​)
.
On rappelle que la fonction sqrt du module math permet de calculer la racine carrée d'un nombre donné en argument.
from math import sqrt
def archimede(k):
    # La liste La contiendra tous les termes de la suite (a_n)
    La = [0] * (k)
    La[0] = ...
    # La liste Lb contiendra tous les termes de la suite (b_n)
    Lb = [0] * (k)
    Lb[0] = ...
    for i in range(1, k):
        Lb[i] = ...
        La[i] = ...
    return La, Lb
La, Lb = archimede(4)
print('Polygones inscrits : ', La)
print('Polygone circonscrits : ', Lb)

Méthode d'Archimède - Représentation graphique des suites a et b

Vous pouvez utiliser l'animation suivante pour visualiser la convergence des suites
(an)(a_n)(an​)
 et
(bn)(b_n)(bn​)
 vers
π\piπ
.
import matplotlib.pyplot as plt
from matplotlib.patches import Circle,RegularPolygon
from IPython.display import HTML
import matplotlib.animation
from math import pi, cos
## Nombre d'étapes à représenter
N = 5
abscisse = [2**(i+2) for i in range(N)]
inf,sup = archimede(N)
fig, (ax1, ax2) = plt.subplots(1, 2,figsize=(12, 6))
l= .7
ax1.set_xlim(( -l, l))
ax1.set_ylim((-l, l))
ax2.set_xlim(( 4, 2**(N+1)))
ax2.set_ylim((2.8, 4))
cercle = Circle((0, 0), .5, facecolor='none',
                edgecolor=(0, 0, 0), linewidth=2, alpha=0.75)
courbeInf, = ax2.plot([],[],'-o',color="#1e7fcb")
courbeSup, = ax2.plot([],[],'-o',color='orange')
def init():
    return []
def animate(i):
    ax1.clear()
    ax1.set_xlim(( -l, l))
    ax1.set_ylim((-l, l))
    ax1.add_patch(cercle)
    long = 0.5/cos(pi/(4*2**i))
    PI = RegularPolygon(numVertices = 4*2**i,xy=(0, 0), radius=.5, orientation=0.79,edgecolor="#1e7fcb", facecolor='none',
                linewidth=2, alpha=0.5)
    PS = RegularPolygon((0, 0), 4*2**i, radius=long, orientation=.79, facecolor='none',
                edgecolor='orange', linewidth=2, alpha=0.5)
    ax1.add_patch(PI)
    ax1.add_patch(PS)
    ax1.set_title('{} côtés'.format(4*2**i),color="#1e7fcb",fontsize=14)
    courbeInf.set_data(abscisse[:i+1], inf[:i+1])
    courbeSup.set_data(abscisse[:i+1], sup[:i+1])
    return PI,
ax2.plot([4,2**(N+1)],[pi,pi],'--',color='green')
ax2.legend(['Polygones intérieurs','Polygones extérieurs','$\pi$'])
plt.close ()
ani = matplotlib.animation.FuncAnimation(fig, animate, init_func=init,  frames=N, blit=False, interval=750)
# l'un ou l'autre
HTML(ani.to_jshtml())
#HTML(ani.to_html5_video())

Méthode d'Archimède - Aspect théorique

Cette approximation met en jeu la notion de suites adjacentes.
Deux suites
(an)(a_n)(an​)
et
(bn)(b_n)(bn​)
sont dites adjacentes si l'une d'elles est croissante, l'autre décroissante et si
lim⁡n→+∞(an−bn)=0\displaystyle\lim_{n \to + \infty}(a_n-b_n)=0n→+∞lim​(an​−bn​)=0
.
On admet alors que deux suites adjacentes sont convergentes et de même limite. Ce résultat est abordé dans l'une des perles d'exercices de la section « Dans la prochaine saison ».
Cet exercice a pour objectif de démontrer que les suites
(an)(a_n)(an​)
 et 
(bn)(b_n)(bn​)
définies dans cette activité sont en effet adjacentes.
On rappelle que l'on a
b0=4b_0=4b0​=4
, 
a0=22a_0=2\sqrt{2}a0​=22​
et, pour tout entier naturel
n⩾2n\geqslant 2n⩾2
,
bn+1=2anbnan+bnb_{n+1}=\dfrac{2a_nb_n}{a_n+b_n}bn+1​=an​+bn​2an​bn​​
 et
an+1=anbn+1a_{n+1}=\sqrt{a_nb_{n+1}}an+1​=an​bn+1​​
.
1.Montrer par récurrence que, pour tout entier naturel
n⩾2n\geqslant 2n⩾2
, on a
0⩽an⩽bn0 \leqslant a_n \leqslant b_n0⩽an​⩽bn​
.
2.Montrer que la suite
(an)(a_n)(an​)
est croissante et que la suite
(bn)(b_n)(bn​)
est décroissante.
3.Montrer que, pour tout entier naturel
nnn
, 
bn+1−an+1⩽bn−an2b_{n+1}-a_{n+1} \leqslant \dfrac{b_n-a_n}{2}bn+1​−an+1​⩽2bn​−an​​
.
4.En déduire
lim⁡n→+∞(an−bn)\displaystyle\lim_{n \to + \infty}(a_n-b_n)n→+∞lim​(an​−bn​)
 et conclure.
5. Écrire une fonction approx_pi en Python qui prend en argument un entier 
nnn
et qui donne un encadrement de
π\piπ
d'amplitude
10−n10^{-n}10−n
 en utilisant la méthode d'Archimède.
Il est bien évidemment interdit de comparer les valeurs des suites
(an)(a_n)(an​)
 et 
(bn)(b_n)(bn​)
directement à la valeur de
π\piπ
!

Méthode de l'arctangente

Une approximation de 
\pi
est donnée par la formule suivante, liée à la trigonométrie : 
π=4lim⁡n→+∞(11−13+45−17+⋯+(−1)n2n+1)=4lim⁡n→+∞∑k=1n(−1)k2k+1\pi = 4 \displaystyle\lim_{n\to+\infty} \left(\dfrac{1}{1}-\dfrac{1}{3}+\dfrac{4}{5}-\dfrac{1}{7} + \dots + \dfrac{(-1)^n}{2n+1}\right) = 4\displaystyle\lim_{n\to+\infty}\sum_{k=1}^{n}\dfrac{(-1)^k}{2k+1}π=4n→+∞lim​(11​−31​+54​−71​+⋯+2n+1(−1)n​)=4n→+∞lim​k=1∑n​2k+1(−1)k​
Il est aussi possible de définir une suite
(pn)(p_n)(pn​)
par
p0=4p_0=4p0​=4
et, pour tout entier naturel non nul
nnn
,
pn=pn−1+4×(−1)n2n+1p_{n}=p_{n-1} + \dfrac{4\times (-1)^n}{2n+1}pn​=pn−1​+2n+14×(−1)n​
. On a alors
π=lim⁡n→+∞pn\pi = \displaystyle\lim_{n\to+\infty}p_nπ=n→+∞lim​pn​
.
Exercice
Compléter la fonction pi_somme1 qui prend en argument un entier
nnn
et qui renvoie la liste des
nnn
premiers termes de la suite
(pn)(p_n)(pn​)
.
def pi_somme1(n):
    Lp = [0] * n
    Lp[0] = ...
    for i in range(1, n):
        Lp[i] = ...
    return Lp
approx = pi_somme1(5)
print(approx)
Calculer alors
p10000p_{10000}p10000​
. Le comparer avec la valeur de
π\piπ
. Que peut-on dire de la vitesse de convergence de cette suite ?
Le programme suivant permet de représenter graphiquement les 100 premiers termes de la suite
(pn)(p_n)(pn​)
.
import matplotlib.pyplot as plt
from math import pi
abscisses = list(range(100))
ordonnees = pi_somme1(100)
plt.plot(abscisses, ordonnees, '.')
plt.plot([0,100],[pi,pi],'--',color='green')
plt.legend(['$p_n$','$\pi$'])
plt.title('100 premiers termes de la suite $(p_n)$')
plt.show()
plt.close()

Suite moyenne

Bien que la convergence de cette suite soit lente, il est possible de l'utiliser pour construire de nouvelles suites dont la convergence vers
π\piπ
est plus rapide.
Pour tout entier naturel
nnn
, on pose alors
qn=pn+pn+12q_n = \dfrac{p_n+p_{n+1}}{2}qn​=2pn​+pn+1​​
,
rn=qn+qn+12r_n = \dfrac{q_n+q_{n+1}}{2}rn​=2qn​+qn+1​​
,
sn=rn+rn+12s_n = \dfrac{r_n+r_{n+1}}{2}sn​=2rn​+rn+1​​
,
tn=sn+sn+12t_n = \dfrac{s_n+s_{n+1}}{2}tn​=2sn​+sn+1​​
... ce procédé pouvant être répété autant de fois que l'on souhaite.
Exercice 1
Justifier que les suites ainsi construites convergent toutes vers
π\piπ
.
Exercice 2
On souhaite déterminer la valeur de 
t30t_{30}t30​
.
1.Combien de termes de la suite
(pn)(p_n)(pn​)
faut-il calculer ?
2.En comptant l'ensemble des suites en jeu, combien de termes faut-il calculer en tout pour déterminer la valeur de
t30t_{30}t30​
?
3.Le programme suivant permet de calculer les
nnn
premiers termes de la suite
(tn)(t_n)(tn​)
. Expliquer brièvement le rôle de la fonction moyenne_suite.
def moyenne_suite(L):
    n = len(L)
    suite = [0] * (n-1)
    for i in range(n-1):
        suite[i] = (L[i]+L[i+1])/2
    return suite
def pi_somme2(n) :
    Lp = pi_somme1(n)
    Lq = moyenne_suite(Lp)
    Lr = moyenne_suite(Lq)
    Ls = moyenne_suite(Lr)
    Lt = moyenne_suite(Ls)
    return Lt
4.À l'aide de ce programme, calculer la valeur de
t30t_{30}t30​
. Comparer la précision de l'approximation de
π\piπ
 entrecettevaleur et la valeur de
p10000p_{10000}p10000​
précédemment calculée.

Autres méthodes

Le problème de Bâle

Le problème de Bâle est un problème mathématique qui consiste à déterminer la limite de la somme
∑k=1n1k2=1+122+132+⋯+1n2\displaystyle\sum_{k=1}^{n}\dfrac{1}{k^2}=1+\dfrac{1}{2^2}+\dfrac{1}{3^2}+\dots +\dfrac{1}{n^2}k=1∑n​k21​=1+221​+321​+⋯+n21​
.
Bien qu'établir la convergence d'une telle somme soit accessible avec les notions vues en classe de terminale, déterminer la valeur de la limite de cette somme utilise d'autres notions plus avancées.
En 1735, le mathématicien Leohnard Euler est alors le premier à démontrer que :  
lim⁡n→+∞∑k=1n1k2=π26\displaystyle\lim_{n \to +\infty}\sum_{k=1}^{n}\dfrac{1}{k^2}=\dfrac{\pi ^2}{6}n→+∞lim​k=1∑n​k21​=6π2​
.
Exercice
En utilisant le résultat précédent, écrire une fonction bale en Python qui permet d'obtenir une valeur approchée de
π\piπ
. Comparer alors la vitesse de convergence de cette somme avec les méthodes précédemment utilisées.

Algorithme de Chudnovsky

Les méthodes de calcul modernes se basent sur une autre somme dont la convergence est plus rapide. Cette formule a notamment permis d'établir un nouveau record de calcul des décimales de 
π\piπ
le 14 mars 2024, puisque ce ne sont pas moins de 105 000 milliards de chiffres du développement décimal qui ont été calculés à l'aide de la formule suivante : 
1π=12∑k=0+∞(−1)k(6k)!(545140134k+13591409)(3k)!(k!)3(640320)3k+3/2\dfrac{1}{\pi}=12\displaystyle\sum_{k=0}^{+\infty} \dfrac{(-1)^k(6k)!(545140134k+13591409)}{(3k)!(k!)^3(640320)^{3k+3/2}}π1​=12k=0∑+∞​(3k)!(k!)3(640320)3k+3/2(−1)k(6k)!(545140134k+13591409)​
Ce calcul a tout de même duré 75 jours, du 14 décembre 2023 au 27 février 2024 !
Le programme ci-dessous permet de calculer les 10 premiers termes de l'inverse de cette somme.
import decimal
def binary_split(a, b):
    if b == a + 1:
        Pab = -(6*a - 5)*(2*a - 1)*(6*a - 1)
        Qab = 10939058860032000 * a**3
        Rab = Pab * (545140134*a + 13591409)
    else:
        m = (a + b) // 2
        Pam, Qam, Ram = binary_split(a, m)
        Pmb, Qmb, Rmb = binary_split(m, b)
        Pab = Pam * Pmb
        Qab = Qam * Qmb
        Rab = Qmb * Ram + Pam * Rmb
    return Pab, Qab, Rab
def chudnovsky(n):
    """Chudnovsky algorithm."""
    P1n, Q1n, R1n = binary_split(1, n)
    return (426880 * decimal.Decimal(10005).sqrt() * Q1n) / (13591409*Q1n + R1n)
print(chudnovsky(2))  # 3.141592653589793238462643384
decimal.getcontext().prec = 100
for n in range(2,10):
    print(f"{n=} {chudnovsky(n)}")  # 3.14159265358979323846264338...