15B — Les algorithmes du programme

1Recherche de seuil — les suites

On connaît une suite \((u_n)\) et un seuil \(S.\) La question est : à partir de quel rang a-t-on \(u_n\geqslant S\) ? On ne sait pas d'avance combien de tours il faudra : c'est exactement ce qu'on cherche, donc la boucle est non bornée.

La même chose en Python, écrite comme une fonction réutilisable, avec le premier terme, la raison et le seuil en paramètres :

def seuil(u0, q, S):
    u = u0
    n = 0
    while u < S:
        u = q * u
        n = n + 1
    return n

print(seuil(1, 2, 1000))

Contrôle à la main : \(2^9=512<1000\) et \(2^{10}=1024\geqslant1000\), donc le premier rang qui atteint \(1000\) est bien \(10.\)

2Dichotomie — continuité et valeurs intermédiaires

Le théorème des valeurs intermédiaires garantit qu'une fonction continue qui change de signe entre \(a\) et \(b\) s'annule quelque part entre les deux. Il ne dit pas où. La dichotomie le trouve : on coupe l'intervalle en deux, on garde la moitié où le changement de signe subsiste, et on recommence.

def f(x):
    return x**3 + x - 1

def dichotomie(a, b, e):
    n = 0
    while b - a > e:
        m = (a + b) / 2
        if f(a) * f(m) <= 0:
            b = m
        else:
            a = m
        n = n + 1
    return a, b, n

print(dichotomie(0, 1, 0.001))

La solution \(\alpha\) de \(x^3+x-1=0\) vérifie donc \(0{,}681640625 < \alpha < 0{,}6826171875\), et il a fallu \(10\) étapes. Les deux bornes affichées sont exactes : ce sont les fractions \(\dfrac{698}{1024}\) et \(\dfrac{699}{1024}.\)

3Méthode de Newton — dérivation et continuité

La dichotomie n'utilise que le signe de \(f.\) La méthode de Newton utilise en plus la tangente : on remplace la courbe par sa tangente, et on prend le point où cette droite coupe l'axe des abscisses.

def f(x):
    return x * x - 2

def fp(x):
    return 2 * x

def newton(x, n):
    for k in range(n):
        x = x - f(x) / fp(x)
        print(k + 1, x)
    return x

newton(2, 4)

On cherche \(\sqrt2\), c'est-à-dire la solution positive de \(x^2-2=0.\) En quatre tours, on obtient \(1{,}41421356\ldots\) : la valeur \(\sqrt2\approx1{,}4142135624\) est retrouvée avec dix décimales justes, alors que la dichotomie en aurait demandé une trentaine d'étapes.

4Méthode d'Euler — exponentielle et équations différentielles

On connaît la pente de la courbe en chaque point, par l'équation \(y'=f(x,y)\), et un point de départ. La méthode d'Euler avance à petits pas en suivant à chaque fois la tangente : c'est encore une histoire de tangente, mais pour construire une courbe entière au lieu de chercher un zéro.

a) Construire l'exponentielle. La fonction exponentielle est définie par \(y'=y\) et \(y(0)=1.\) La pente en \(y_k\) vaut donc \(y_k\) lui-même, et \(y_{k+1}=y_k+h\,y_k=(1+h)\,y_k.\)

def euler_exp(n):
    h = 1 / n
    y = 1
    for k in range(n):
        y = y + h * y
    return y

print(euler_exp(4))

Avec quatre pas de \(h=0{,}25\) sur \([0\,;1]\), on obtient \((1{,}25)^4=2{,}44140625\) : c'est une valeur exacte de la ligne brisée, et une valeur approchée de \(\mathrm e\approx2{,}718.\) La ligne brisée reste en dessous de la courbe, parce que l'exponentielle est convexe.

b) Résoudre \(y'=ay+b.\) La même boucle sert pour n'importe quelle équation de ce type : il suffit de changer la pente.

def euler(a, b, y0, h, n):
    y = y0
    for k in range(n):
        y = y + h * (a * y + b)
    return y

print(euler(-2, 6, 1, 0.25, 4))

Pour \(y'=-2y+6\) avec \(y(0)=1\), la récurrence devient \(y_{k+1}=y_k+0{,}25\,(-2y_k+6)=0{,}5\,y_k+1{,}5.\) Les valeurs successives sont \(1\), \(2\), \(2{,}5\), \(2{,}75\), \(2{,}875\) : elles montent vers la solution constante \(y=3\), qui est bien la valeur d'équilibre \(-\dfrac ba=-\dfrac{6}{-2}=3.\)

5Méthode des rectangles — calcul intégral

L'intégrale \(\displaystyle\int_a^b f(x)\,\mathrm dx\) est une aire. On la découpe en \(n\) tranches de largeur \(h=\dfrac{b-a}{n}\), et on remplace chaque tranche par un rectangle : une fois en prenant la hauteur à gauche de la tranche, une fois à droite.

def f(x):
    return x * x

def rectangles(a, b, n):
    h = (b - a) / n
    g = 0
    d = 0
    for k in range(n):
        g = g + h * f(a + k * h)
        d = d + h * f(a + (k + 1) * h)
    return g, d

for n in [4, 16, 64, 256]:
    print(n, rectangles(0, 1, n))

La valeur exacte est \(\displaystyle\int_0^1 x^2\,\mathrm dx=\dfrac13\), et elle est bien à l'intérieur des quatre encadrements. Les quatre nombres affichés sont exacts, sans arrondi : les bornes de l'encadrement sont des fractions de dénominateur une puissance de \(2.\)

6Simulation d'un échantillon, marche aléatoire — probabilités

Simuler, c'est demander à l'ordinateur de jouer l'expérience aléatoire, un grand nombre de fois, pour observer des fréquences. L'instruction de base est random.random(), qui rend un nombre au hasard entre \(0\) et \(1\), chaque intervalle ayant une probabilité égale à sa longueur. Un succès de probabilité \(p\) se simule donc par le test random.random() < p.

a) Un échantillon de taille \(n.\) Un échantillon de taille \(n\), c'est la liste de \(n\) observations indépendantes de la même variable aléatoire.

import random

def tirage(p):
    if random.random() < p:
        return 1
    return 0

def echantillon(n, p):
    return [tirage(p) for k in range(n)]

random.seed(2026)
E = echantillon(20, 0.5)
print(E)
print(sum(E) / 20)
print(sum(echantillon(1000, 0.5)) / 1000)

Sur \(20\) tirages, la fréquence observée vaut \(\dfrac{6}{20}=0{,}3\) : loin de \(0{,}5.\) Sur \(1000\) tirages elle vaut \(\dfrac{501}{1000}=0{,}501\) : la loi des grands nombres est à l'œuvre. Une fréquence observée n'est jamais la probabilité ; elle s'en rapproche quand \(n\) grandit.

b) Marche aléatoire. Un pion part de \(0\) ; à chaque étape il avance de \(1\) ou recule de \(1\), à pile ou face. Où est-il après \(n\) pas ?

import random

def marche(n):
    x = 0
    for k in range(n):
        if random.random() < 0.5:
            x = x + 1
        else:
            x = x - 1
    return x

random.seed(7)
print([marche(10) for i in range(5)])

Les cinq positions sont paires, et ce n'est pas un hasard : après \(10\) pas, la position est \(a-b\) avec \(a+b=10\), donc \(a-b=2a-10\) est toujours de la même parité que \(10.\) C'est un contrôle de vraisemblance immédiat du programme (méthode D).

7Fréquence simulée contre majoration de Bienaymé-Tchebychev

Le programme cite nommément cet algorithme : « calculer la probabilité de \(\left(|S_n-pn|>\sqrt n\right)\), où \(S_n\) est une variable aléatoire qui suit une loi binomiale \(\mathcal B(n,p)\). Comparer avec l'inégalité de Bienaymé-Tchebychev ». On simule la probabilité, on calcule la majoration, et on regarde l'écart entre les deux.

import random

def binomiale(n, p):
    s = 0
    for k in range(n):
        if random.random() < p:
            s = s + 1
    return s

def frequence(N, n, p, a):
    c = 0
    for i in range(N):
        if abs(binomiale(n, p) - n * p) >= a:
            c = c + 1
    return c / N

random.seed(2026)
print(frequence(10000, 100, 0.5, 10))
print(100 * 0.5 * 0.5 / (10 * 10))

Avec \(n=100\), \(p=0{,}5\) et \(\delta=\sqrt{100}=10\) (le texte officiel écrit l'inégalité stricte \(>\) ; la fonction frequence teste l'inégalité large >=, celle qui figure dans Bienaymé-Tchebychev), la fréquence observée de l'événement \(|S_n-50|\geqslant10\) vaut \(\dfrac{581}{10\,000}=0{,}0581\), alors que la majoration vaut \(\dfrac{100\times0{,}5\times0{,}5}{10^2}=\dfrac{25}{100}=0{,}25.\)