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\geqslant 2

, on note

a_n

le périmètre du polygone à

2^n

côtés inscrit dans le cercle et

b_n

le périmètre du polygone à

2^n

côtés circonscrit au cercle.

Exercice 1

Justifier que l'on a 

a_2=2\sqrt{2}

et

b_2=4

.

Exercice 2

On admet que, pour tout entier naturel

n\geqslant 2

, on a

b_{n+1}=\dfrac{2a_nb_n}{a_n+b_n}

et

a_{n+1}=\sqrt{a_nb_{n+1}}

.
Calculer 

b_3

,

a_3

,

b_4

et

a_4

.

Exercice 3

Compléter le programme suivant qui renvoie les 

k

premiers termes des suites

(a_n)

et

(b_n)

.
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

(a_n)

et

(b_n)

 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

(a_n)

et

(b_n)

sont dites adjacentes si l'une d'elles est croissante, l'autre décroissante et si

\displaystyle\lim_{n \to + \infty}(a_n-b_n)=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

(a_n)

et

(b_n)

définies dans cette activité sont en effet adjacentes.

On rappelle que l'on a

b_0=4

,

a_0=2\sqrt{2}

et, pour tout entier naturel

n\geqslant 2

,

b_{n+1}=\dfrac{2a_nb_n}{a_n+b_n}

et

a_{n+1}=\sqrt{a_nb_{n+1}}

.

1.Montrer par récurrence que, pour tout entier naturel

n\geqslant 2

, on a

0 \leqslant a_n \leqslant b_n

.

2.Montrer que la suite

(a_n)

est croissante et que la suite

(b_n)

est décroissante.

3.Montrer que, pour tout entier naturel

n

,

b_{n+1}-a_{n+1} \leqslant \dfrac{b_n-a_n}{2}

.

4.En déduire

\displaystyle\lim_{n \to + \infty}(a_n-b_n)

 et conclure.

5. Écrire une fonction approx_pi en Python qui prend en argument un entier 

n

et qui donne un encadrement de

\pi

d'amplitude

10^{-n}

 en utilisant la méthode d'Archimède.
Il est bien évidemment interdit de comparer les valeurs des suites

(a_n)

et

(b_n)

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 : 

\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}

Il est aussi possible de définir une suite

(p_n)

par

p_0=4

et, pour tout entier naturel non nul

n

,

p_{n}=p_{n-1} + \dfrac{4\times (-1)^n}{2n+1}

. On a alors

\pi = \displaystyle\lim_{n\to+\infty}p_n

.

Exercice

Compléter la fonction pi_somme1 qui prend en argument un entier

n

et qui renvoie la liste des

n

premiers termes de la suite

(p_n)

.

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

p_{10000}

. 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

(p_n)

.

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

n

, on pose alors

q_n = \dfrac{p_n+p_{n+1}}{2}

,

r_n = \dfrac{q_n+q_{n+1}}{2}

,

s_n = \dfrac{r_n+r_{n+1}}{2}

,

t_n = \dfrac{s_n+s_{n+1}}{2}

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

t_{30}

.
1.Combien de termes de la suite

(p_n)

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

t_{30}

?
3.Le programme suivant permet de calculer les

n

premiers termes de la suite

(t_n)

. 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

t_{30}

. Comparer la précision de l'approximation de

\pi

 entrecettevaleur et la valeur de

p_{10000}

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

\displaystyle\sum_{k=1}^{n}\dfrac{1}{k^2}=1+\dfrac{1}{2^2}+\dfrac{1}{3^2}+\dots +\dfrac{1}{n^2}

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

\displaystyle\lim_{n \to +\infty}\sum_{k=1}^{n}\dfrac{1}{k^2}=\dfrac{\pi ^2}{6}

.

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 : 

\dfrac{1}{\pi}=12\displaystyle\sum_{k=0}^{+\infty} \dfrac{(-1)^k(6k)!(545140134k+13591409)}{(3k)!(k!)^3(640320)^{3k+3/2}}

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

Activité CAPYTALE : approximation de Pi

Quelques rappels pour commencer

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

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

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

Méthode d'Archimède

Principe de la méthode d'Archimède

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

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

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

Méthode de l'arctangente

Suite moyenne

Autres méthodes

Le problème de Bâle

Algorithme de Chudnovsky