Python · Niveau 5 — Pour aller plus loin

Mode d'emploi

Statut de ce niveau

Ce niveau est en cours d'audit par M. Van Hoorde. Il est utilisable, mais pas encore garanti à 100 %.

Ce niveau suppose acquis les niveaux 1 à 4. Ses sept étapes sont des ouvertures : elles ne seront évaluées ni à l'ECE ni au bac. Chacune porte une étiquette qui indique la formation où cette notion est au programme.

Les étapes se suivent : l'étape 3 réutilise l'étape 2, l'étape 4 l'étape 3, l'étape 6 l'étape 5.
Tous les signaux sont fabriqués par le programme, avec numpy : aucun fichier à charger.
L'étape 7 utilise la bibliothèque scipy. Elle se télécharge à son premier usage (14 Mo) : comptez quelques dizaines de secondes, une seule fois.
Comme aux niveaux précédents : 💡 Aide pour un indice, ⓘ Rappel de cours pour la notion du cours, solution après un premier essai.

Numériser une tension

BUT Mesures physiques · R2.08

Un capteur délivre une tension. Pour la traiter en Python, une carte d'acquisition la convertit en nombre entier : c'est le rôle du convertisseur analogique-numérique (CAN). Le programme du BUT Mesures physiques demande de « savoir numériser un signal analogique et connaître les caractéristiques métrologiques du signal numérisé ».

Rappel de cours — Numérisation (enseignement scientifique de 1re)

Une information codée sur n bits peut prendre 2ⁿ valeurs. La résolution de la numérisation est le nombre de bits sur lesquels la valeur est codée.

L'écart entre deux niveaux voisins s'appelle le pas p. Pour une plage de mesure de 0 à U_max :

p = U_max2ⁿ − 1

L'imprécision due à la numérisation n'excède jamais la moitié du pas.

À retenir

Le CAN remplace la tension u par l'entier le plus proche de u / p : N = np.round(u / p). La tension restituée vaut N * p.
np.round(x) arrondit à l'entier le plus proche, même un tableau entier d'un coup.
np.abs(x) donne la valeur absolue, np.max(x) la plus grande valeur d'un tableau.

1.1 La carte micro:bit de l'ECE

D'après ECE PC, banque 2026, sujet n° 62. Une carte micro:bit mesure la tension aux bornes d'un condensateur. Son CAN code la tension sur 10 bits, pour une plage de 0 à 3,3 V. La carte renvoie un entier ; le programme du sujet le convertit en volts.

Combien de valeurs entières différentes ce CAN peut-il renvoyer ?

Compléter la ligne p = …… qui calcule le pas de ce CAN.
La carte renvoie l'entier 512. Compléter la ligne qui calcule la tension correspondante.
Quelle est la plus petite variation de tension que la carte peut détecter ?

n = 10           # nombre de bits du CAN
U_max = 3.3      # plage de mesure : de 0 à 3,3 V

p = ……           # pas, en V
print('p =', p, 'V')

N = 512          # entier renvoyé par la carte
u = ……           # tension correspondante, en V
print('u =', round(u, 3), 'V')

Voir une solution

p = U_max / (2 ** n - 1)
u = N * p

p ≈ 0,003 23 V, soit 3,2 mV : c'est la plus petite variation détectable. L'entier 512 correspond à u ≈ 1,652 V. C'est exactement la ligne du sujet : Uc = Uc_NUM * 3.3 / 1023.

1.2 Un son en escalier

Un microphone délivre une tension sinusoïdale de fréquence 440 Hz, comprise entre 0,15 et 3,15 V. On la numérise d'abord avec un CAN de 3 bits seulement.

Compléter les lignes N = …… et u_num = ……, puis exécuter.
Compléter la ligne erreur = …… : le plus grand écart, en valeur absolue, entre u_num et u. Le comparer à la moitié du pas.
Passer à 10 bits. Que deviennent l'escalier et l'erreur ?

Toutes les lignes s'écrivent avec les fonctions de l'encadré À retenir.
L'écart entre les deux tableaux est u_num - u ; on en prend la valeur absolue, puis la plus grande valeur.
Les fonctions s'emboîtent : np.max(np.abs(……)).

import numpy as np
import matplotlib.pyplot as plt

n = 3                                # nombre de bits du CAN
U_max = 3.3                          # plage : de 0 à 3,3 V
p = U_max / (2 ** n - 1)             # pas, en V

t = np.linspace(0, 0.005, 500)                  # 5 ms, en s
u = 1.65 + 1.5 * np.sin(2 * np.pi * 440 * t)    # tension du micro, en V

N = ……           # entiers renvoyés par le CAN
u_num = ……       # tensions restituées, en V
erreur = ……      # plus grand écart entre u_num et u, en V

print('pas :', round(p, 4), 'V')
print('erreur maximale :', round(erreur, 4), 'V')

plt.plot(t * 1000, u, 'b-', label='signal analogique')
plt.plot(t * 1000, u_num, 'r-', label='signal numérisé')
plt.xlabel('t (ms)')
plt.ylabel('u (V)')
plt.legend()
plt.show()

Voir une solution

N = np.round(u / p)
u_num = N * p
erreur = np.max(np.abs(u_num - u))

Sur 3 bits, le pas vaut 0,47 V : l'escalier n'a que 8 marches et l'erreur atteint 0,236 V, la moitié du pas. Sur 10 bits, le pas tombe à 3,2 mV et l'erreur à 1,6 mV : la courbe rouge se confond avec la bleue. En contrepartie, chaque mesure occupe plus de place en mémoire.

Échantillonner un signal

attendu de PC · MP · PSI · BUT Mesures physiques (R2.08)

Numériser, c'est aussi choisir quand mesurer. Mesurer trop rarement ne fait pas seulement perdre des détails : cela peut faire apparaître une fréquence qui n'existe pas. Le programme de PC (2e année) demande de « mettre en évidence le phénomène de repliement du spectre provoqué par l'échantillonnage » et d'« utiliser la condition de Nyquist-Shannon ».

Rappel de cours — Échantillonnage (enseignement scientifique de 1re)

La carte d'acquisition mesure le signal à intervalles de temps réguliers : la période d'échantillonnage T_e. Le nombre de mesures par seconde est la fréquence d'échantillonnage f_e.

f_e = 1T_e

Le format CD : f_e = 44,1 kHz, sur 16 bits.

À retenir

np.arange(0, duree, Te) crée les valeurs 0, T_e, 2T_e… jusqu'à duree exclue : ce sont les instants des mesures.
Condition de Nyquist-Shannon : pour reproduire fidèlement un signal, il faut f_e > 2 × f_max, où f_max est la plus grande fréquence qu'il contient.
Sinon, il y a repliement de spectre : une sinusoïde de fréquence f comprise entre f_e / 2 et f_e est vue à la fréquence f_e − f.

2.1 Un La mal échantillonné

Le programme trace un La de 440 Hz avec 20 000 points : c'est le signal « analogique ». On l'échantillonne ensuite à la fréquence fe.

Compléter les lignes t_ech = …… et u_ech = ……, puis exécuter avec f_e = 8 000 Hz.
Remplacer f_e par 500 Hz. Les points rouges dessinent-ils encore le La ?

import numpy as np
import matplotlib.pyplot as plt

f = 440            # fréquence du La, en Hz
duree = 0.05       # durée observée, en s

t = np.linspace(0, duree, 20000)      # signal « analogique » : presque continu
u = np.sin(2 * np.pi * f * t)

fe = 8000          # fréquence d'échantillonnage, en Hz
t_ech = ……         # instants des mesures
u_ech = ……         # valeurs mesurées

plt.plot(t * 1000, u, 'b-', linewidth=0.8, label='signal')
plt.plot(t_ech * 1000, u_ech, 'ro-', label='échantillons')
plt.xlabel('t (ms)')
plt.ylabel('u')
plt.legend()
plt.show()

Avec fe = 500 Hz

À quelle fréquence la courbe rouge semble-t-elle osciller ?

Voir une solution

t_ech = np.arange(0, duree, 1 / fe)
u_ech = np.sin(2 * np.pi * f * t_ech)

À 8 000 Hz, on mesure environ 18 fois par période : les points suivent le signal. À 500 Hz, on mesure un peu plus d'une fois par période seulement : les points rouges dessinent une sinusoïde de 60 Hz qui n'existe pas. La condition de Nyquist-Shannon n'est pas respectée : 500 Hz < 2 × 440 Hz.

2.2 Le format CD et les ultrasons du TP5

Une oreille humaine jeune entend des sons de 20 Hz à 20 kHz. Utiliser la console pour les calculs si besoin.

Nyquist-Shannon et repliement

Pour enregistrer fidèlement tous les sons audibles, la fréquence d'échantillonnage doit être…

Au format CD, quelle est la plus haute fréquence enregistrée sans repliement ?

L'émetteur d'ultrasons du TP5 (40 kHz) est enregistré à 44,1 kHz, sans filtre. À quelle fréquence apparaît-il ?

Le spectre d'un son

attendu de PSI · en PC, Fourier « non exigible »

En 1re, un logiciel traçait le spectre d'un son. On le calcule maintenant soi-même avec la transformée de Fourier discrète. Seul le programme de PSI l'exige : « calculer la transformée de Fourier discrète d'un signal à valeurs réelles en utilisant la fonction rfft de la bibliothèque numpy.fft ». En PC, il précise que « la transformée de Fourier n'est pas exigible ».

Rappel de cours — Spectre d'un son (enseignement scientifique de 1re)

Un signal périodique de fréquence f₁ peut être décomposé en une somme de signaux sinusoïdaux : un de fréquence f₁, la fréquence fondamentale, et d'autres de fréquences f_n = n × f₁, les harmoniques.

L'ensemble de ces fréquences constitue le spectre du son. Un son pur, sinusoïdal, n'a qu'un pic. La hauteur d'un son est sa fréquence fondamentale ; les harmoniques font le timbre.

À retenir

X = np.fft.rfft(x) calcule la transformée de Fourier discrète du signal échantillonné x : un tableau de nombres complexes.
A = np.abs(X) * 2 / N, avec N = len(x), donne l'amplitude de chaque composante sinusoïdale.
freq = np.fft.rfftfreq(N, 1 / fe) donne les fréquences correspondantes, de 0 à f_e / 2. Le second paramètre est la période d'échantillonnage T_e.
Deux fréquences voisines du spectre sont séparées de f_e / N = 1 / durée : un enregistrement d'une seconde donne un spectre précis au hertz.
np.argmax(A) donne le rang de la plus grande valeur de A ; freq[np.argmax(A)] est donc la fréquence du pic le plus haut.

3.1 Le spectre du diapason

Un diapason émet un son pur de 440 Hz. On simule une seconde d'enregistrement, échantillonnée à 8 000 Hz.

Compléter les lignes X = …… et A = …….
Compléter la ligne freq = ……, puis exécuter. Relever la fréquence et l'amplitude du pic.

import numpy as np
import matplotlib.pyplot as plt

fe = 8000                                # fréquence d'échantillonnage, en Hz
t = np.arange(0, 1.0, 1 / fe)            # une seconde de mesures
x = 0.8 * np.sin(2 * np.pi * 440 * t)    # le son du diapason

N = len(x)
X = ……             # transformée de Fourier discrète
A = ……             # amplitudes
freq = ……          # fréquences, en Hz

print('fréquence du pic :', freq[np.argmax(A)], 'Hz')
plt.plot(freq, A, 'b-')
plt.xlim(0, 2000)                        # on n'affiche que de 0 à 2 000 Hz
plt.xlabel('f (Hz)')
plt.ylabel('amplitude')
plt.show()

Voir une solution

X = np.fft.rfft(x)
A = np.abs(X) * 2 / N
freq = np.fft.rfftfreq(N, 1 / fe)

Un seul pic, à 440 Hz, d'amplitude 0,8 : un son pur. Le spectre compte 4 001 valeurs, de 0 à 4 000 Hz par pas de 1 Hz.

3.2 Une note et ses harmoniques

On fabrique une note de fondamentale 220 Hz en additionnant ses six premières composantes. Leurs amplitudes sont rangées dans la liste a (valeurs construites pour l'exercice) : a[0] pour le rang 1, a[1] pour le rang 2…

Compléter la ligne de la boucle qui ajoute au signal l'harmonique de rang n.
Exécuter : le signal est-il sinusoïdal ? Lire le spectre.

L'harmonique de rang n est une sinusoïde de fréquence n × f₁.
Son amplitude est rangée dans a[n - 1], puisqu'une liste commence au rang 0.
La ligne a la forme x = x + amplitude * np.sin(2 * np.pi * fréquence * t).

import numpy as np
import matplotlib.pyplot as plt

fe = 8000
t = np.arange(0, 1.0, 1 / fe)
f1 = 220                                 # fréquence fondamentale, en Hz
a = [1.0, 0.5, 0.0, 0.25, 0.0, 0.12]     # amplitudes des rangs 1 à 6

x = np.zeros(len(t))                     # on part d'un signal nul
for n in range(1, 7):                    # n prend les valeurs 1, 2… 6
    x = ……

plt.plot(t[:160] * 1000, x[:160], 'b-')  # les 20 premières millisecondes
plt.xlabel('t (ms)')
plt.ylabel('signal')
plt.show()

N = len(x)
A = np.abs(np.fft.rfft(x)) * 2 / N
freq = np.fft.rfftfreq(N, 1 / fe)
plt.plot(freq, A, 'r-')
plt.xlim(0, 2000)
plt.xlabel('f (Hz)')
plt.ylabel('amplitude')
plt.show()

Lire le spectre

Quelle est la hauteur de cette note ?

Quelles harmoniques sont absentes du spectre, parmi les six premières ?

Voir une solution

for n in range(1, 7):
    x = x + a[n - 1] * np.sin(2 * np.pi * n * f1 * t)

Le signal se répète toutes les 4,5 ms (1 / 220 s), mais il n'est pas sinusoïdal : c'est un son complexe. Le spectre retrouve exactement les amplitudes de la liste : la transformée de Fourier défait l'addition faite par la boucle.

3.3 Un spectre trompeur

Une note de fondamentale 1 500 Hz possède des harmoniques jusqu'au rang 4. On l'enregistre à 8 000 Hz.

Exécuter. Quelles fréquences attendait-on ? Lesquelles s'affichent ?
Expliquer l'origine des pics intrus.
Modifier la fréquence d'échantillonnage pour obtenir le bon spectre.

import numpy as np
import matplotlib.pyplot as plt

fe = 8000
t = np.arange(0, 1.0, 1 / fe)
f1 = 1500
a = [1.0, 0.6, 0.4, 0.3]                 # amplitudes des rangs 1 à 4

x = np.zeros(len(t))
for n in range(1, 5):
    x = x + a[n - 1] * np.sin(2 * np.pi * n * f1 * t)

N = len(x)
A = np.abs(np.fft.rfft(x)) * 2 / N
freq = np.fft.rfftfreq(N, 1 / fe)
print('pics à', freq[A > 0.1], 'Hz')    # les fréquences où l'amplitude dépasse 0,1

plt.plot(freq, A, 'b-')
plt.xlabel('f (Hz)')
plt.ylabel('amplitude')
plt.show()

Le pic à 2 000 Hz

D'où vient le pic à 2 000 Hz, avec f_e = 8 000 Hz ?

Voir une solution

On attendait 1 500, 3 000, 4 500 et 6 000 Hz. On lit 1 500, 2 000, 3 000 et 3 500 Hz : les rangs 3 et 4 dépassent f_e / 2 et sont repliés. Rien ne les distingue de vraies composantes : un spectre ne se lit qu'en connaissant la fréquence d'échantillonnage.

fe = 16000

Toute valeur supérieure à 12 000 Hz convient. Le spectre s'étend alors jusqu'à 8 000 Hz et montre les quatre pics attendus.

Derrière le rideau

Calculer directement la transformée de Fourier de N valeurs demande de l'ordre de N² multiplications : 64 millions pour une seconde à 8 000 Hz. En 1965, James Cooley et John Tukey publient une méthode qui n'en demande que de l'ordre de N × log₂ N, environ 100 000 ici : la transformée de Fourier rapide, FFT en anglais. C'est le « fft » de np.fft. On a découvert plus tard que Carl Friedrich Gauss avait trouvé la même astuce vers 1805, sans la publier. La FFT tourne aujourd'hui dans chaque téléphone, dans les formats MP3 et JPEG, et dans les accordeurs de guitare.

Filtrer un signal bruité

attendu de MP · TSI2

Une mesure réelle est toujours bruitée. Un filtre numérique la nettoie par le calcul, après l'acquisition. Le programme de MP demande de « simuler un filtrage numérique et visualiser son action sur un signal périodique » ; celui de TSI2, de « réaliser […] un filtrage numérique passe-bas d'un signal issu d'une acquisition ».

À retenir

Moyenne glissante : on remplace chaque valeur par la moyenne des 2k + 1 valeurs qui l'entourent. Les k premières et les k dernières sont perdues.
T[i - k : i + k + 1] contient les valeurs de rang i - k à i + k : une tranche s'arrête juste avant son second indice.
liste.append(v) ajoute la valeur v au bout d'une liste.
Filtre passe-bas du premier ordre : il laisse passer les fréquences basses et atténue les fréquences élevées. Sa version numérique est une récurrence : chaque valeur filtrée se calcule à partir de la précédente.

4.1 Lisser une courbe de refroidissement

Rappel de cours — Refroidissement d'un système

Un système de température initiale T₀, placé dans l'air de température T_ext, se refroidit. Selon la loi de Newton, sa température évolue selon :

T(t) = T_ext + (T₀ − T_ext) × e^{−t / τ}

τ est le temps caractéristique du refroidissement.

Un thermomètre numérique relève chaque seconde, pendant 10 min, la température d'une eau chaude qui refroidit de 80 °C vers 20 °C. Le capteur est bruité. Les mesures sont simulées : le modèle de Newton, plus un bruit tiré au hasard.

Compléter la ligne de la boucle, puis exécuter.
Passer à k = 50. Qu'y gagne-t-on ? Qu'y perd-on ?

np.random.seed(2) fixe le point de départ du hasard : les « mesures » sont les mêmes à chaque exécution, ce qui permet de comparer deux réglages.

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(2)
t = np.arange(0, 601, 1.0)                    # une mesure par seconde, en s
T = 20 + 60 * np.exp(-t / 300) + np.random.normal(0, 0.4, len(t))   # en °C

k = 10                                        # demi-largeur de la fenêtre
T_liss = []
for i in range(k, len(T) - k):
    T_liss.append(……)                         # moyenne des 2k + 1 valeurs autour de i
t_liss = t[k : len(T) - k]                    # les instants correspondants

plt.plot(t, T, 'b.', markersize=2, label='mesures')
plt.plot(t_liss, T_liss, 'r-', label='moyenne glissante')
plt.xlabel('t (s)')
plt.ylabel('T (°C)')
plt.legend()
plt.show()

Voir une solution

    T_liss.append(np.mean(T[i - k : i + k + 1]))

Avec k = 10, chaque point est la moyenne de 21 mesures : le bruit est divisé par environ √21 ≈ 4,6. Avec k = 50, la courbe est plus lisse encore, mais on perd 50 s au début, là où la température varie le plus vite, et 50 s à la fin. Lisser, c'est toujours choisir entre le bruit et la perte d'information.

4.2 Un filtre passe-bas, comme un circuit RC

Rappel de cours — Circuit RC

Lors de la charge d'un condensateur de capacité C à travers un conducteur ohmique de résistance R, sous une tension E, la tension u_C vérifie :

RC × du_Cdt + u_C = E

τ = RC est le temps caractéristique du circuit.

Un circuit RC alimenté par un signal x(t) quelconque vérifie τ × dy/dt + y = x, où y est la tension aux bornes du condensateur. La méthode d'Euler, vue au niveau 3, donne la récurrence :

y_i = y_i−1 + T_eτ × (x_i−1 − y_i−1)

On l'applique à un son de 200 Hz mêlé d'un sifflement de 5 000 Hz (valeurs construites pour l'exercice). Le temps caractéristique est choisi pour couper vers 1 000 Hz.

Compléter la ligne de la boucle avec la récurrence, puis exécuter.
Comparer les deux spectres : que devient chaque composante ?

import numpy as np
import matplotlib.pyplot as plt

fe = 50000                                 # fréquence d'échantillonnage, en Hz
Te = 1 / fe
t = np.arange(0, 0.04, Te)
x = np.sin(2 * np.pi * 200 * t) + 0.5 * np.sin(2 * np.pi * 5000 * t)   # son + sifflement

tau = 1 / (2 * np.pi * 1000)               # temps caractéristique, en s
y = np.zeros(len(x))                       # signal filtré, nul au départ
for i in range(1, len(x)):
    y[i] = ……

plt.plot(t[:500] * 1000, x[:500], 'b-', label='avant')
plt.plot(t[:500] * 1000, y[:500], 'r-', label='après')
plt.xlabel('t (ms)')
plt.legend()
plt.show()

freq = np.fft.rfftfreq(len(x), Te)
Ax = np.abs(np.fft.rfft(x)) * 2 / len(x)
Ay = np.abs(np.fft.rfft(y)) * 2 / len(y)
plt.plot(freq, Ax, 'b-', label='avant')
plt.plot(freq, Ay, 'r-', label='après')
plt.xlim(0, 6000)
plt.xlabel('f (Hz)')
plt.ylabel('amplitude')
plt.legend()
plt.show()

Lire les spectres

Par combien le filtre divise-t-il, à peu près, l'amplitude du sifflement ?

Voir une solution

    y[i] = y[i - 1] + Te / tau * (x[i - 1] - y[i - 1])

Le signal filtré est presque une sinusoïde de 200 Hz, légèrement retardée. Ce calcul est exactement la méthode d'Euler appliquée au circuit RC : un filtre numérique imite un circuit électrique. Pour qu'il reste fidèle, il faut T_e petit devant τ, ici 8 fois plus petit : c'est la limite introduite par l'échantillonnage dont parle le programme de TSI2.

La diffusion thermique : un tableau à deux dimensions

attendu de PC · MP · TSI2

Comment le froid de l'hiver traverse-t-il un mur ? La température dépend de deux variables : la position dans le mur et le temps. On la range dans un tableau à deux dimensions. Le programme de PC demande de « résoudre l'équation de la diffusion thermique à une dimension par une méthode des différences finies dérivée de la méthode d'Euler explicite ».

À retenir

T = np.zeros((3, 5)) crée un tableau de 3 lignes et 5 colonnes rempli de zéros. T.shape vaut (3, 5).
T[j, i] est l'élément de la ligne j, colonne i, en comptant à partir de 0.
T[j, :] est toute la ligne j ; T[:, i] toute la colonne i.
L'indice -1 désigne le dernier élément. T[j, 1:-1] est la ligne j sans son premier ni son dernier élément.
Une ligne entière s'affecte d'un coup : T[0, :] = 20 met 20 dans toute la ligne 0.

M = np.array([[1, 2, 3],
              [4, 5, 6]])     # 2 lignes, 3 colonnes
print(M[1, 0])                # affiche 4
print(M[:, 2])                # affiche [3 6]

5.1 Se repérer dans un tableau

Dans la console, taper import numpy as np, puis :

M = np.array([[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12]])

Prédire chaque résultat, puis le vérifier dans la console.

Que vaut…

M.shape

M[1, 2]

M[:, 0]

M[0, 1:-1]

M[-1, -1]

5.2 Un mur en hiver

Rappel de cours — Transfert thermique par conduction

Un transfert thermique se fait spontanément du corps le plus chaud vers le corps le plus froid. Dans un solide, il se fait par conduction, de proche en proche, sans déplacement de matière.

À travers une paroi de résistance thermique R_th, dont les faces sont aux températures T₁ et T₂, le flux thermique vaut, en régime permanent :

Φ = T₁ − T₂R_th

Un mur de béton de 20 cm d'épaisseur est à 20 °C. Brusquement, l'air extérieur passe à 0 °C ; à l'intérieur, le chauffage maintient 20 °C. On découpe l'épaisseur en points espacés de Δx, le temps en pas Δt, et on range dans T[j, i] la température du point i à l'instant j. En 2e année, on montre que chaque point intérieur évolue selon :

T[j+1, i] = T[j, i] + r × (T[j, i+1] − 2 T[j, i] + T[j, i−1]) avec r = D × ΔtΔx²

D est la diffusivité thermique du matériau, en m²·s⁻¹ ; pour le béton, de l'ordre de 8 × 10⁻⁷ m²·s⁻¹. Chaque point se rapproche de la moyenne de ses deux voisins. Pour calculer tous les points intérieurs d'un coup, T[j, 2:] contient les voisins de droite et T[j, :-2] les voisins de gauche.

Compléter les trois lignes qui fixent la température initiale et celle des deux faces.
Compléter la ligne de la boucle, puis exécuter.
Décrire l'évolution du profil de température.

À l'instant 0, c'est la ligne 0 qui vaut 20 : T[0, :] = 20.
La face extérieure est la colonne 0, à tous les instants ; la face intérieure est la dernière colonne.
Dans la boucle, chaque T[j, i] de la formule devient T[j, 1:-1], chaque T[j, i+1] devient T[j, 2:], chaque T[j, i-1] devient T[j, :-2].

import numpy as np
import matplotlib.pyplot as plt

L = 0.20                   # épaisseur du mur, en m
nx = 41                    # nombre de points dans l'épaisseur
dx = L / (nx - 1)          # pas d'espace : 5 mm
D = 8e-7                   # diffusivité thermique du béton, en m²/s
dt = 12.5                  # pas de temps, en s
r = D * dt / dx ** 2
nt = int(12 * 3600 / dt)   # nombre de pas pour 12 h

T = np.zeros((nt + 1, nx))     # T[j, i] : point i, instant j
……                             # à l'instant 0, tout le mur est à 20 °C
……                             # face extérieure (colonne 0) : 0 °C à tout instant
……                             # face intérieure (dernière colonne) : 20 °C à tout instant

for j in range(nt):
    T[j + 1, 1:-1] = ……

x = np.linspace(0, L, nx)
for h in [0, 0.25, 1, 3, 6]:
    j = int(h * 3600 / dt)                    # rang de l'instant h heures
    plt.plot(x * 100, T[j, :], label=str(h) + ' h')
plt.xlabel('position dans le mur (cm)')
plt.ylabel('T (°C)')
plt.legend()
plt.show()
print('milieu du mur après 1 h :', round(T[288, 20], 1), '°C')

Lire le graphique

Au bout de 6 h, le profil de température dans le mur est…

Voir une solution

T[0, :] = 20
T[:, 0] = 0
T[:, -1] = 20

for j in range(nt):
    T[j + 1, 1:-1] = T[j, 1:-1] + r * (T[j, 2:] - 2 * T[j, 1:-1] + T[j, :-2])

Au début, seuls les premiers centimètres se refroidissent. Au bout d'une heure, le milieu du mur est à 16,2 °C ; au bout de quelques heures, le profil devient une droite de 0 à 20 °C : le régime permanent. Le tableau T compte 3 457 lignes et 41 colonnes : toute l'histoire du mur.

5.3 Quand le calcul explose

On reprend le mur, en cherchant à calculer plus vite avec un pas de temps plus grand.

Exécuter avec dt = 12.5, puis avec dt = 16. Observer la valeur de r et le graphique.
Trouver le plus grand pas de temps utilisable.

Le schéma explicite n'est stable que si r ≤ 0,5. Au-delà, les petites erreurs d'arrondi sont amplifiées à chaque pas.

import numpy as np
import matplotlib.pyplot as plt

L, nx, D = 0.20, 41, 8e-7
dx = L / (nx - 1)
dt = 12.5                  # pas de temps, en s
r = D * dt / dx ** 2
print('r =', round(r, 3))
nt = int(12 * 3600 / dt)

T = np.zeros((nt + 1, nx))
T[0, :] = 20
T[:, 0] = 0
T[:, -1] = 20
for j in range(nt):
    T[j + 1, 1:-1] = T[j, 1:-1] + r * (T[j, 2:] - 2 * T[j, 1:-1] + T[j, :-2])

x = np.linspace(0, L, nx)
for h in [1, 3, 6, 12]:
    plt.plot(x * 100, T[int(h * 3600 / dt), :], label=str(h) + ' h')
plt.xlabel('position dans le mur (cm)')
plt.ylabel('T (°C)')
plt.legend()
plt.show()

Stabilité

Avec dt = 16, que se passe-t-il ?

Quel est le plus grand pas de temps utilisable ?

Derrière le rideau

C'est pour résoudre cette équation de la chaleur que Joseph Fourier, préfet de l'Isère, invente en 1807 la décomposition d'une fonction en somme de sinusoïdes : celle de l'étape 3. Son mémoire est d'abord jugé peu rigoureux par l'Académie des sciences ; il ne paraît sous sa forme définitive qu'en 1822, dans la Théorie analytique de la chaleur. Deux siècles plus tard, les séries de Fourier servent partout où il y a des ondes, du son à l'imagerie médicale.

La marche aléatoire

attendu de PC

Une goutte de sirop de menthe déposée dans de l'eau immobile finit par la colorer entièrement, sans qu'on remue. Chaque molécule de colorant est bousculée au hasard par les molécules d'eau. Le programme de PC demande de « simuler la marche au hasard d'un grand nombre de particules à partir d'un centre et caractériser l'étalement spatial ».

À retenir

Modèle à une dimension : à chaque pas, une particule avance ou recule d'une longueur a, au hasard, avec la même probabilité.
np.random.choice([-1, 1], size=(n, p)) crée un tableau de n lignes et p colonnes, rempli au hasard de −1 et de +1 : une ligne par particule, une colonne par pas.
np.sum(M, axis=1) fait la somme de chaque ligne de M ; axis=0 ferait la somme de chaque colonne.
L'étalement du nuage de particules se mesure par l'écart-type de leurs positions.

6.1 Une goutte de colorant

2 000 molécules de colorant partent toutes de la position 0 et font 100 pas de longueur a = 1 (en unité arbitraire).

Compléter les lignes pas = …… et x = ……, puis exécuter.
Quelle est la forme de l'histogramme ? Combien vaut l'écart-type ?

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(3)
n_part = 2000          # nombre de particules
n_pas = 100            # nombre de pas

pas = ……               # tableau n_part × n_pas de −1 et de +1
x = ……                 # position finale de chaque particule

print('écart-type des positions :', round(np.std(x, ddof=1), 1))
plt.hist(x, bins=30, edgecolor='black')
plt.xlabel('position finale')
plt.ylabel('nombre de particules')
plt.show()

Voir une solution

pas = np.random.choice([-1, 1], size=(n_part, n_pas))
x = np.sum(pas, axis=1)

L'histogramme a la forme d'une cloche centrée sur 0, celle de la loi normale du niveau 4. L'écart-type vaut environ 10 : après 100 pas, les molécules se sont étalées de 10 pas en moyenne, et non de 100. Elles vont et viennent.

6.2 L'étalement grandit comme √t

On recommence pour 25, 100, 400 et 1 600 pas, et on range les écarts-types obtenus dans la liste ecarts.

Compléter les deux lignes de la boucle, puis exécuter.
Que montre le graphique ? Combien de fois plus de temps faut-il pour s'étaler deux fois plus ?

Pour chaque valeur de n_pas, refaire la simulation de l'exercice 6.1.
La position finale peut s'écrire en une ligne : x = np.sum(np.random.choice(……), axis=1).
Ajouter ensuite l'écart-type de x à la liste : ecarts.append(……).

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(4)
n_part = 2000
liste_pas = [25, 100, 400, 1600]
ecarts = []
for n_pas in liste_pas:
    x = ……                      # positions finales après n_pas pas
    ecarts.append(……)           # écart-type de ces positions

print(np.round(ecarts, 1))
plt.plot(np.sqrt(liste_pas), ecarts, 'bo')
plt.xlabel('racine du nombre de pas')
plt.ylabel('écart-type des positions')
plt.grid()
plt.show()

Voir une solution

    x = np.sum(np.random.choice([-1, 1], size=(n_part, n_pas)), axis=1)
    ecarts.append(np.std(x, ddof=1))

Les écarts-types valent environ 5, 10, 20 et 40 : les points sont alignés sur une droite passant par l'origine. L'étalement est proportionnel à la racine du temps : pour s'étaler deux fois plus loin, il faut quatre fois plus de temps. C'est la signature de toute diffusion, celle du mur de l'étape 5 comprise.

Derrière le rideau

En 1905, Albert Einstein explique ainsi le mouvement désordonné de minuscules particules en suspension dans l'eau, observé au microscope par le botaniste Robert Brown en 1827 sur des particules issues de grains de pollen : elles sont bousculées par des molécules, invisibles. Son calcul prévoit un étalement en √t. Jean Perrin le vérifie expérimentalement en 1908 et en déduit une valeur de la constante d'Avogadro. Ces mesures achèvent de convaincre les derniers sceptiques que les atomes existent ; Perrin reçoit le prix Nobel de physique en 1926.

Loi de Student et limite de détection

BUT Chimie · R3.08

Au laboratoire d'analyse, on répète rarement une mesure plus de cinq fois. Avec si peu de valeurs, l'intervalle de confiance de la loi normale est trop optimiste. Et quelle est la plus petite concentration qu'une méthode peut détecter ? Le programme du BUT Chimie cite la « loi de Student » et les « notions de limites de détection et de quantification ».

Rappel de cours — Incertitude de type A

Pour n mesures répétées, de moyenne x̄ et d'écart-type expérimental s(x) :

u(x̄) = s(x)√n

À retenir

Avec n mesures, l'intervalle de confiance à 95 % s'écrit x̄ ± t × u(x̄). Le coefficient de Student t dépend de n − 1, le nombre de degrés de liberté. Il vaut 4,30 pour 3 mesures, 2,78 pour 5, 2,26 pour 10, et tend vers 1,96 ≈ 2, la valeur de la loi normale, quand n devient grand.
from scipy import stats, puis stats.t.ppf(0.975, n - 1) donne t pour 95 % : on laisse 2,5 % de chaque côté de l'intervalle.
On mesure plusieurs fois un blanc, solution sans l'espèce dosée ; s_blanc est l'écart-type de ces mesures et k la pente de la droite d'étalonnage. Selon la convention la plus courante, la limite de détection vaut LD = 3 s_blanc / k et la limite de quantification LQ = 10 s_blanc / k.
Sous LD, on ne peut pas affirmer que l'espèce est présente. Entre LD et LQ, elle est détectée, mais sa concentration n'est pas fiable.

7.1 Cinq mesures seulement

On mesure cinq fois l'absorbance du même sirop de menthe dilué, comme au TP3 (valeurs construites pour l'exercice).

Compléter la ligne t = …… : le coefficient de Student à 95 %.
Compléter la ligne demi = …… : la demi-largeur de l'intervalle de confiance. La comparer à 2u.

import numpy as np
from scipy import stats

A = np.array([0.471, 0.468, 0.475, 0.466, 0.472])   # absorbances
n = len(A)
A_moy = np.mean(A)
u = np.std(A, ddof=1) / np.sqrt(n)

t = ……             # coefficient de Student à 95 %
demi = ……          # demi-largeur de l'intervalle de confiance

print('A =', round(A_moy, 4), '±', round(demi, 4), '(Student, 95 %)')
print('avec 2u :', round(2 * u, 4))

Voir une solution

t = stats.t.ppf(0.975, n - 1)
demi = t * u

A = 0,4704 ± 0,0044 avec Student, contre ± 0,0031 avec 2u. Avec cinq mesures seulement, l'écart-type est lui-même mal connu : Student élargit l'intervalle pour en tenir compte. Avec 30 mesures, t = 2,05 et la différence devient négligeable.

7.2 Jusqu'où peut-on doser ?

Rappel de cours — Loi de Beer-Lambert

L'absorbance A d'une solution est proportionnelle à la concentration c de l'espèce absorbante : A = k × c. Une concentration inconnue se lit sur la droite d'étalonnage.

On dose le colorant bleu d'un sirop de menthe. Le programme contient la gamme étalon et dix mesures du blanc (valeurs construites pour l'exercice). np.polyfit, vu au niveau 3, donne la pente k et l'ordonnée à l'origine b de la droite d'étalonnage.

Compléter les lignes s_blanc, LD et LQ, puis exécuter.
Un échantillon très dilué a une absorbance de 0,012. Conclure.

import numpy as np

c = np.array([0.0, 1.0, 2.0, 3.0, 4.0, 5.0])                # en mg/L
A = np.array([0.001, 0.126, 0.247, 0.373, 0.495, 0.621])    # absorbances
k, b = np.polyfit(c, A, 1)                                  # droite d'étalonnage

blanc = np.array([0.002, -0.001, 0.003, 0.000, 0.001,
                  -0.002, 0.002, 0.001, -0.001, 0.000])     # 10 mesures du blanc

s_blanc = ……       # écart-type expérimental du blanc
LD = ……            # limite de détection, en mg/L
LQ = ……            # limite de quantification, en mg/L

A_ech = 0.012
c_ech = (A_ech - b) / k
print('k =', round(k, 4), 'L/mg')
print('LD =', round(LD, 3), 'mg/L   LQ =', round(LQ, 3), 'mg/L')
print('échantillon :', round(c_ech, 3), 'mg/L')

Conclure

L'échantillon, à environ 0,09 mg/L, est…

Voir une solution

s_blanc = np.std(blanc, ddof=1)
LD = 3 * s_blanc / k
LQ = 10 * s_blanc / k

k ≈ 0,124 L·mg⁻¹, s_blanc ≈ 0,0016, LD ≈ 0,038 mg·L⁻¹ et LQ ≈ 0,128 mg·L⁻¹. L'échantillon, à 0,089 mg·L⁻¹, est détecté mais non quantifiable. Valider une méthode d'analyse, c'est justement établir ces limites avant de rendre le moindre résultat.

Bilan du niveau 5

Ce que je sais faire

Calculer le pas d'un CAN et numériser une tension, en connaissant l'erreur commise.
Choisir une fréquence d'échantillonnage avec la condition de Nyquist-Shannon, et reconnaître un repliement.
Calculer et lire le spectre d'un son avec np.fft.rfft et np.fft.rfftfreq.
Lisser une mesure par moyenne glissante et programmer un filtre passe-bas par récurrence.
Ranger une grandeur dans un tableau à deux dimensions et résoudre l'équation de la diffusion par différences finies, en respectant la condition de stabilité.
Simuler une marche aléatoire et montrer que l'étalement croît comme √t.
Écrire un intervalle de confiance avec le coefficient de Student, et calculer les limites de détection et de quantification d'un dosage.

Ces outils sont ceux de la 2e année de classe préparatoire et des formations de techniciens supérieurs. Vous les retrouverez avec leur documentation, comme ici : un bon programmeur ne retient pas tout, il sait où chercher.