Physique numérique en Python

Journée de formation

Lycée Camille Jullian

Mars 2018
Table des matières

1 Loi horaire 1
1.1 Position du problème . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Mise en équation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.3 Méthodologie de tracé . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.4 Modules et packages . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.5 Codage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.6 Compléments graphiques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.7 À vous de jouer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

2 Frottements fluides 5
2.1 Position du problème . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.2 Étude de la vitesse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.3 Schémas d’Euler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.4 Méthodologie de tracé . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.5 Codage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.6 Solution exacte . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.7 Solution Python . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.8 À vous de jouer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

3 Trajectoire 11
3.1 Position du problème . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
3.2 Mise en équation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
3.3 Méthodologie de tracé . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
3.4 Codage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
3.5 Système différentiel et Python . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
3.6 À vous de jouer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

4 Oscillateurs 15
4.1 Position du problème . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
4.2 Notations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
4.3 À vous de jouer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

5 Canon de newton 17
5.1 Position du problème . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
5.2 Mise en équation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
5.3 À vous de jouer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
5.4 Données numériques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

6 Petit vade-mecum 19
6.1 Tracé d’une courbe à partir d’une fonction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
6.2 Tracé d’une courbe à partir de données tabulaires . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
6.3 Ajustements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
Introduction

L’objet de ce document est d’illustrer l’utilisation de Python dans le cadre de problèmes de physique. Les thèmes
proposés s’appuient largement sur les contenus du programme de physique de TS. Partant de situations simples
rencontrées en TS, on présente quelques pistes d’utilisation de résolution numérique avec Python. Certains exemples
admettant une solution analytique, une confrontation avec les résultats numériques permet de mesurer la pertinence
d’une approche informatique. D’autres exemples plus complexes, sans solution analytique simple pour des élèves de
TS, voire sans solution analytique du tout, sont traités par approche numérique. Ces exemples peuvent être envisagés
comme des prolongements du programme.

Une première série de trois activités propose de se familiariser avec le graphisme élémentaire de Python et avec
les méthodes numériques de résolution numérique d’équations différentielles. Le problème de la chute sert de fil
conducteur à ces quatre activités.

Une quatrième activité vise à mettre en application les compétences acquises avec des activités précédentes. Le fil
conducteur est celui de l’oscillateur : oscillateur linéaire d’abord, oscillateur forcé ensuite, oscillateur non-linéaire
enfin.

La cinquième activité illustre le problème du canon de Newton. À quelle condition une masse lancée depuis une
certaine altitude, dans une certaine direction, peut-elle se satelliser ou échapper à l’attraction terrestre ? Cette activité
est l’occasion d’illustrer la mécanique gravitationnelle dans la situation simple à deux corps en interaction. Elle requiert
toute les compétences et la technicité des activités précédentes.

Un très grand merci à mes collègues Christophe Casseau et Marc Eldin pour leur aide à la construction et à la relecture
de ce document.
Un grand merci également à David Boyer, IA-IPR de sciences-physiques, à l’initiative de ces journées, pour nous avoir
réuni et donné l’occasion de monter cette première formation qui, nous l’espérons, sera suivie d’autres rencontres
autour de la physique et de l’informatique. À suivre. . .

Laurent Sartre
[Link]@[Link]

iii
Physique numérique et Python iv

Journée de formation - mars 2018 cbea

 1.4 Modules et packages
Python est un langage de programmation dont les fonc-
tionnalités peuvent facilement étendues à l’aide de nou-
velles instructions. Par défaut, il met à la disposition de
l’utilisateur un jeu significativement important d’instruc-
tions mais pas toujours suffisant pour des besoins plus
Activité 1 spécifiques. Tracer des courbes avec Python entre dans la
catégorie des besoins particuliers.
Pour étendre les fonctionnalités du langage, de nouvelles
commandes peuvent être créées et conservées dans des
Loi horaire fichiers particuliers appelés modules. L’intérêt de ces mo-
dules est de pouvoir être réutilisés en fonction des besoins.
Parfois, des ensembles de modules sont regroupés pour
former des packages. Il s’agit simplement d’une hiérarchi-
sation des modules.
1.1 Position du problème En pratique, dès qu’un module connu est nécessaire, il
doit être importé au début du script grâce à l’instruction
L’objectif de cette activité est de se familiariser avec
import.
quelques méthodes graphiques Python en vue de tracer
Dans ce qui suit, le module numpy, dédié au calcul numé-
une courbe représentative de la loi horaire d’un phéno-
rique, est importé. Il permet en particulier la manipulation
mène physique.
de tableaux multidimensionnels pour stocker des données.
Pour illustrer notre propos, on s’intéresse à la chute libre
Ce module fournit également un grand nombre de fonc-
d’un corps lâché sans vitesse initiale et pour lequel l’in-
tions dédiées aux calculs mathématiques 1 et numériques.
fluence des frottements peut être négligée. Cette situation,
à seul degré de liberté, peut être modélisée par une équa- # importation de numpy
tion différentielle. Sa solution analytique peut être aisé- # alias np
ment trouvée en TS. Python est utilisé pour tracer la loi import numpy as np
horaire de l’altitude du corps.
Ainsi, dans le script suivant, le module pyplot du package
1.2 Mise en équation matplotlib est importé. Ce module enrichit Python d’un
Un point matériel M de masse m est soumis à la seule certain nombre de fonctionnalités graphiques dont cer-
force de pesanteur. On note g la norme de l’accélération taines d’entre elles sont présentées dans la suite de l’acti-
de la pesanteur. Le mouvement de M est étudié dans un vité.
référentiel local supposé galiléen. La verticale du lieu défi- # importation de [Link]
nit l’axe (Oz ), vertical ascendant. # alias plt
Le point M est initialement lâché d’une hauteur h > 0, sans import [Link] as plt
vitesse initiale. Le schéma ci-dessous illustre la situation.
À un instant t ⩾ 0, on note z (t ) la cote du point M . Le L’ajout de l’instruction as suivie d’un nom particulier auto-
principe fondamental de la dynamique permet d’établir la rise la création d’un alias qui remplace avantageusement
loi d’évolution de z sous la forme d’une équation différen- le nom du module. Son intérêt réside dans le fait que les
tielle. instructions des modules sont accessibles par une nota-
∀ t ∈ R+
tion pointée : chaque instruction d’un module donné ne
z̈ (t ) = −g (1.1)
peut être appelée qu’en tapant son nom précédé du nom
La solution de (1.1) est : du module suivi d’un point. Par exemple, pour tracer la
courbe d’un tableau de valeurs y en fonction d’un tableau
de valeurs x, on code comme suit.
∀ t ∈ R+
1
z (t ) = h − g t 2 (1.2)
2 [Link](x,y)

1.3 Méthodologie de tracé L’affichage de la courbe se fait en codant comme suit.

Pour tracer la courbe représentative de la loi horaire de z
[Link]()
avec Python, on doit :

◾ définir les paramètres physiques g et h en précisant Il convient de savoir qu’il est possible d’éviter la définition
leurs valeurs numériques ; et l’utilisation d’alias en important globalement toutes les
instructions d’un module. Il suffit de réaliser l’import de la
◾ définir la fonction z ∶ t ↦ h − g t 2 /2 ;
manière suivante.
◾ préciser l’intervalle de temps [t min , t max ] sur lequel
from [Link] import *
on désire tracer la courbe ;
from numpy import *
◾ tracer la courbe et l’afficher.
1. comme les fonctions trigonométriques, la racine carrée, etc.

1
Physique numérique et Python 2

Cette importation est à utiliser avec prudence. Certains # tableau des instants
modules définissent de nouvelles commandes qui portent tab_t = [Link](t_min,t_max,n_t)
le même nom que d’autres commandes qui préexistent
dans Python. Dès lors, la dernière définition écrase toute
L’une des forces du module numpy est de permettre la
autre définition préalable. Ce qui peut présenter des incon-
construction d’un tableau en appliquant une fonction sur
vénients. Dans la suite, nous utilisons systématiquement
un autre tableau. Le script suivant affecte à la variable
la notation pointée associée à la définition d’alias.
tab_z le tableau obtenu en appliquant la fonction z au
1.5 Codage tableau tab_t et aux variables g et h. Cette facilité est l’un
des points forts de Python pour le calcul scientifique.
Pour tracer la courbe représentative de la loi horaire de
z, les modules nécessaires sont d’abord chargés. Le script # tableau des z
suivant permet le tracé de la loi horaire avec h = 5 m sur un tab_z = z(tab_t,g,h)
intervalle de temps [0, 1 s] 2 . La courbe obtenue est celle
de la figure 1.1.
Il reste à tracer la courbe en portant, dans un système
# importation de numpy d’axes orthogonaux, les valeurs présentes dans la tableau
# alias np tab_z en fonction des valeurs présentes dans le tableau
import numpy as np tab_t. Noter que le nombre de valeurs dans chaque ta-
# importation de [Link] bleau doit être le même pour que cette opération ait du
# alias plt sens.
import [Link] as plt
# tracé de tab_z en fonction de tab_t
On définit ensuite une fonction qui rend compte de la loi # option 'r-' : tracé en rouge (r)
horaire. Dans le script suivant, cette fonction est notée z # en reliant les points (-)
et trois arguments doivent lui être transmis pour qu’elle [Link](tab_t,tab_z,'r-')
renvoie un résultat. Cela permettra de faire tracer des tra- # affichage de la courbe
jectoires en faisant varier les paramètres physiques h et [Link]()
3
g .
Le résultat de l’exécution de la séquence des instructions
# fonction associée à la loi horaire
précédentes est le graphique suivant.
def z(t,g,h):
return h - g * t**2 / 2

Vient à présent la définition des paramètres physiques,

dont les valeurs sont affectées à des variables informa-
tiques. Le nom de ces variables est arbitraire mais il ne
peut pas débuter par un chiffre. En pratique, on attribue
des noms pertinents au regard du rôle que les variables
son censées jouer.

# paramètres physiques
g = 9.81 # accélération de la pesanteur
h = 1.0 # hauteur initiale en m

Afin de tracer la loi horaire, on doit se fixer un intervalle de

temps. C’est généralement le sens physique qui prévaut
pour déterminer cet intervalle. Dans le script qui suit, on
choisit des instants compris entre 0 s et 5 s.
F IGURE 1.1 – Tracé d’une chute libre.
# tracé des trajectoires
t_min = 0.0 # instant initial en s
t_max = 1.0 # instant final en s
n_t = 100 # nombre de points de calculs 1.6 Compléments graphiques
Sur le graphique précédent, aucune information ne per-
Pour tracer la loi horaire, on choisit un nombre de points met de connaître la nature des informations portées sur les
n_t de calculs dans l’intervalle de temps précédent. Il ne axes. Il n’y a également pas de titre, ni d’information sur les
reste plus qu’à faire calculer par Python les différents ins- valeurs des paramètres h et t max . Ces information peuvent
tants équirépartis dans cet intervalle et les valeurs de z être ajoutées en faisant précéder l’instruction [Link]
associées. L’instruction linspace construit un tableau de d’un certain nombre d’instructions dédiées à l’affiche de
n_t instants allant de t_min à t_max. Ce tableau est af- commentaires.
fecté à la variable tab_t.
2. Pourquoi ces choix de valeurs ?
3. Il est possible de travailler avec des variables informatiques globales, mais là encore, pour des raisons de prudence, il convient re réduire leur
utilisation au maximum.

Journée de formation - mars 2018 cbea

Physique numérique et Python 3

[Link]('Temps (en seconde)') [Link](tab_t,tab_z,'r-',label="h = " +

[Link]('Position de $M$ (en mètre)') ↪ str(h) + " m")
[Link]("Chute libre sans frottement et
↪ sans vitesse initiale") Cette option reçoit trois chaînes chaîne de caratères mises
bout à bout 4 : la chaîne "h = ", la chaîne str(h) où l’ins-
[Link](0.7,4.8,r'$h = 5$ mètres') truction str convertit la valeur contenue dans la variable
[Link](0.7,4.5,r'$t_{max} = 1$ seconde')
h en une chaîne de caractères, la chaîne " m". Noter la
présence des espaces.
Le résultat graphique est alors le suivant. L’instruction [Link](loc=4) saisie juste avant la
commande d’affichage [Link]() place les commen-
taires précédents dans le coin inférieur droit du graphique.
▶ Question 2. Modifier le nombre de points du tracé
pour afficher seulement quelques points de la courbe,
comme sur la figure 1.4.

F IGURE 1.2

De nombreuses autres options sont disponibles avec

le pyplot du package matplotlib. Le lecteur inté-
ressé peut se reporter à la documentation et aux nom-
breux exemples disponibles sur le site de matplotlib F IGURE 1.4
([Link]
▶ Question 3. Comment obtenir la figure 1.5 ?
1.7 À vous de jouer
▶ Question 1. Modifier le script précédent pour afficher
les deux courbes associées aux valeurs h = 5 m et h = 10 m.

F IGURE 1.5

▶ Question 4. Sur un même graphique, tracer l’évolu-

F IGURE 1.3 tion temporelle de l’énergie cinétique par unité de masse
et de l’énergie potentielle par unité de masse. Tracer égale-
Pour afficher chaque valeur de h et le code couleur de la ment la somme de ces deux quantités.
courbe associée, ajouter l’option label dans l’instruction
plot.
4. En informatique, la mise bout à bout de symboles est appelée la concaténation.

Journée de formation - mars 2018 cbea

Physique numérique et Python 4

Journée de formation - mars 2018 cbea

 De fait, l’équation différentielle précédente peut se mettre
sous la forme suivante.

∀ t ∈ R+
1
v̇ (t ) = g − v (t ) (2.2)
τ
On s’intéresse à présent à la recherche de solutions numé-
riques approchées à cette équation, v (0) = 0 (corps lâché
Activité 2 sans vitesse initiale). Pour mettre en œuvre les schémas
numériques d’Euler, la relation précédente est mise sous
forme intégrée.

Frottements fluides ∀ t ∈ R+ ∫
0
t
v̇ (u ) du = ∫
0
t 1
(g − v (u )) du
τ
Le membre de gauche s’intègre immédiatement sous la
forme :
t
∫ v̇ (u ) du = v (t ) − v (0)
2.1 Position du problème 0
L’objectif de cette activité est de se familiariser avec des Le membre de droite ne s’intègre que partiellement.
méthodes numériques de résolution d’équations différen- L’équation finalement obtenue se met sous la forme in-
tielles. La méthode d’Euler est introduite de trois façons tégrale suivante.
pour en montrer les avantages et les inconvénients. Puis t
∀ t ∈ R+
1
une méthode interne à Python, qui s’appuie sur des mé- v (t ) − v (0) = g t − ∫ v (u ) du (2.3)
τ 0
thodes numériques classiques, est présentée comme un
prolongement à même de répondre aux besoins numé- La condition initiale v (0) = 0 permet de simplifier le
riques du physicien. membre de gauche. Mais pour des raisons pédagogiques, il
L’étude de la chute libre est encore le fil conducteur de peut être intéressant de conserver la présence de ce terme
cette activité. Des frottements fluides sont ajoutés à la si- en vue de préparer la mise en œuvre des techniques nu-
tuation précédente. A priori, cette situation, à un degré de mériques.
liberté encore, peut être modélisée par une équation diffé-
rentielle du second ordre dont l’inconnue est l’altitude du 2.3 Schémas d’Euler
corps. Bien qu’une solution analytique existe, sa détermi- L’idée des schémas numériques d’Euler est de remplacer
nation ne peut cependant pas être trouvée par des élèves l’intégrale précédente par une valeur approchée. Cepen-
de TS. En introduisant la vitesse du corps, on se ramène dant, une telle approximation ne peut être faite sans pré-
à une équation différentielle du premier ordre. Des mé- caution. Elle sera d’autant plus raisonnable que le terme
thodes numériques fournissent des solutions approchées sous le signe intégral varie peu sur l’intervalle d’intégra-
à cette équation. L’enseignant peut alors fournir l’expres- tion. En pratique, il y a peu de chance que cette propriété
sion de la solution analytique afin d’observer l’efficacité soit vérifiée. La vitesse peut varier de manière importante
des méthodes numériques. Une discussion peut s’engager sur l’intervalle d’observation du phénomène. En revanche,
sur la confiance à accorder aux résultats numériques en si cet intervalle est divisé en de nombreux sous-intervalles
vue de sensibiliser les élèves aux problèmes numériques suffisamment petits pour que la vitesse y varie peu, l’ap-
d’une part, de montrer la nécessité d’une attitude critique proximation devient plus raisonnable.
à l’égard des résultats fournis par l’outil numérique d’autre Notons t max l’instant final d’observation de la chute. Dé-
part. Fort de ces compétences, le travail peut être poursuivi coupons la durée ∆t = t max − 0 en n t intervalles de mêmes
pour étudier l’évolution de l’altitude du corps et pour faire durées δt = ∆t /n t . Définissons les instants :
un bilan énergétique numérique.
∀i ∈ {0, 1, . . . , n }, t i = i × δt
2.2 Étude de la vitesse Plus n t est grand, plus δt est petit et plus il y a de valeurs
La force de frottements fluides est modélisée par le vecteur d’instants t i . En augmentant la valeur de n t , on peut rai-
force F⃗f = −λż ê z , λ étant un coefficient strictement positif
sonnablement penser que les intervalles temporels de-
lié aux frottements, ê z désignant un vecteur unitaire ver- viennent suffisamment petits pour que la vitesse y varie
tical ascendant. L’équation donnant l’évolution de z est à peu. L’idée des schémas d’Euler est d’exploiter cette ana-
présent : lyse pour donner une forme pratique à l’équation (2.3)
λ
∀ t ∈ R+ z̈ (t ) = −g − ż (t ) en intégrant sur un intervalle de la forme [t i , t i +1 ], i =
(2.1)
m {0, 1, . . . , n − 1}. Ainsi, il est possible d’écrire, pour tout en-
En raison du choix d’orientation de l’axe vertical, ż est né- tier i prenant ses valeurs dans l’ensemble {0, 1, . . . , n − 1} :
gative : la vitesse du corps augmente en descendant. Il peut 1 t i +1
être commode de poser v = −ż de manière à raisonner sur v (t i +1 ) − v (t i ) = g × δt − ∫ v (u ) du
τ ti
une vitesse positive dont la valeur croît au fur et à mesure
que le corps perd de l’altitude. Par ailleurs, une analyse Si sur chaque intervalle [t i , t i +1 ], la vitesse varie peu, elle
dimensionnelle de l’équation précédente montre que le varie quand même. La question est alors de savoir quelle
rapport m /λ est homogène à un temps caractéristique τ. valeur choisir pour l’approcher.

5
Physique numérique et Python 6

◾ Le schéma d’Euler explicite propose de remplacer illustrant la construction du tableau à l’aide du schéma
v (u ) par v (t i ) sur cet intervalle. L’équation inté- d’Euler explicite, la boucle for parcourt une liste d’entiers
grale prend alors la forme approchée suivante. de 0 à n_t-1. C’est la fonction range(n_t) qui définit
cette liste. Noter que la dernière valeur définie par cette
v (t i ) fonction est n_t-1 et non n_t. D’autres utilisations de la
v (t i +1 ) − v (t i ) ≈ (g − ) × δt
τ fonction range reçoivent plusieurs arguments. Si a et b
sont deux nombres entiers :
◾ Le schéma d’Euler implicite propose de remplacer
v (u ) par v (t i +1 ) sur cet intervalle. L’équation inté- ◾ range(a,b) définit un objet itérable
grale prend alors la forme approchée suivante. [a,a+1,a+2,...,b-1] ;
◾ range(b) définit un objet itérable [0,1,2,...,b-1] ;
v ( t i +1 )
v (t i +1 ) − v (t i ) ≈ (g − ) × δt ◾ range(a,b,c) définit un objet itérable
τ [a,a+c,a+2c,...] ;
Ces relations permettent de calculer, de proche en proche, La création d’un tableau se fait à l’aide de l’instruction
les valeurs approchées de la vitesse aux différents instants [Link], du module numpy. En tapant
t i . Ces valeurs ne sont pas les valeurs exactes et seront no-
tab_v = [Link](n_t)
tées v i ≈ v (t i ) de sorte que les relations établies ci-dessus
se traduisent par les égalités suivantes.
on crée un tableau vide à n_t+1 éléments. Les éléments
v i +1 − v i = (g − τ ) × δt
vi
(Euler explicite) d’un tableau sont accessibles par l’intermédiaire de leur
{ rang dans le tableau, en commençant par le rang 0. Ainsi, à
v i +1 − v i = (g − iτ+1 ) × δt (Euler implicite)
v
chaque rang dans le tableau tab_v est associée une valeur
En isolant v , on obtient des relations de récurrence qui approchée de vitesse selon les correspondances suivantes.
i +1
permettent le calcul de proche en proche des valeurs ap- tab_v[0] ←→ v 0
prochées.
tab_v[1] ←→ v 1
⎧ δt
⎪v i +1 = (1 − τ ) v i + g × δt (Euler explicite)
⎪ ...
⎨
⎪
⎪ v i +1 = δt (v i + g × δt ) (Euler implicite)
1
tab_v[n_t] ←→ v n t
⎩ 1+ τ
Le script suivant présente la fonction euler_explicite
2.4 Méthodologie de tracé qui calcule et retourne le tableau des vitesses.
Pour tracer la courbe représentative de la loi horaire de v
def eulerExp(t_min,t_max,n_t,g,h,v0,tau):
avec Python, on doit : # création d'un tableau vide
◾ définir les paramètres physiques g et h en précisant # des vitesses
leurs valeurs numériques ; # à (n_t + 1) éléments
tab_v = [Link](n_t)
◾ préciser l’intervalle de temps [t min , t max ] sur lequel
# affectation de la vitesse initiale
on désire tracer la courbe ;
# au premier élement du tableau
◾ discrétiser cet intervalle de temps ; tab_v[0] = v0
◾ évaluer de proche en proche les valeurs approchées # calcul de delta_t
de vitesses ; delta_t = (t_max - t_min) / n_t
# boucle de calcul
◾ tracer la courbe et l’afficher. # des vitesses approchées
for i in range(n_t-1):
2.5 Codage tab_v[i+1] = (1 - delta_t / tau) *
Le début du script comporte les mêmes importations que ↪ tab_v[i] + g * delta_t
celles rencontrées dans l’activité 1. return tab_v

# importation de numpy D’autres choix de codage sont possibles. Le choix fait ici
# alias np
est celui de la simplicité et d’un lien étroit avec les expres-
import numpy as np
sions mathématiques établies plus haut.
# importation de [Link]
# alias plt Le script suivant définit les paramètres physiques du pro-
import [Link] as plt blème.
# paramètres physiques
À présent, il nous faut construire l’équivalent de la loi ho- g = 9.81 # accélération de la pesanteur
raire définie par une fonction dans l’activité 1. Une fonc- m = 1.0 # masse
tion Python remplit ce rôle. Elle reçoit un ensemble d’argu- lamb = 0.1 # coef. de frottement
ments nécessaires au calcul des différentes valeurs appro- tau = m / lamb
chées de la vitesse et retourne un tableau de ces valeurs. v0 = 0.0 # vitesse initiale
D’un point de vue informatique, le calcul de proche en h = 5.0 # hauteur initiale
proche associé aux relations de récurrence fait appel à
une structure répétitive (ou boucle). Dans le script suivant Puis les tableaux des instants et des vitesses sont calculés.

Journée de formation - mars 2018 cbea

Physique numérique et Python 7

t_min = 0.0 # instant initial # fonction associée à la loi horaire

t_max = 5.0 * tau # instant final # de la vitesse
n_t = 100 def v(tab_t,g,h,tau):
return g * tau * (1.0 - [Link](-tab_t /
tab_t = [Link](t_min,t_max,n_t) ↪ tau))
tab_v = eulerExp(t_min,t_max,n_t,g,h,v0,tau)

# paramètres physiques
Enfin, la courbe et quelques informations sont affichées. g = 9.81 # accélération de la pesanteur
m = 1.0 # masse
[Link](tab_t,tab_v,'r-',label="h = " +
lamb = 0.1 # coef. de frottement
str(h) + " m")
tau = m / lamb
↪
[Link]('Temps (en s)')
v0 = 0.0 # vitesse initiale
[Link]('Vitesse de $M$ (en m/s)')
h = 5.0 # hauteur initiale
[Link]("Chute libre avec frottements
↪ fluides et sans vitesse initiale")
[Link](0.7 * t_max,10.0,r'$h = ' + str(h) + # tracé des trajectoires
↪ '$ m') t_min = 0.0 # instant initial
[Link](0.7 * t_max,5.0,r'$t_{max} = ' + t_max = 5.0 * tau # instant final
↪ str(t_max) + '$ s') n_t = 100
[Link](0.7 * t_max,15.0,r'$m = ' + str(m) + # tableaux
↪ '$ kg') tab_t = [Link](t_min,t_max,n_t)
[Link](0.7 * t_max,20.0,r'$\lambda = ' + tab_v = v(tab_t,g,h,tau)
↪ str(lamb) + '$ kg/s')# h = 10.0 # hauteur # tracé
↪ initiale [Link](tab_t,tab_v,'g--',label="Solution
[Link]() ↪ exacte")
[Link]('Temps (en s)')
Le résultat graphique est présenté figure 2.1. [Link]('Vitesse de $M$ (en m/s)')
[Link]("Chute libre avec frottements
↪ fluides et sans vitesse initiale")
[Link](0.7 * t_max,40.0,r'$h = ' + str(h) +
↪ '$ m')
[Link](0.7 * t_max,35.0,r'$t_{max} = ' +
↪ str(t_max) + '$ s')
[Link](0.7 * t_max,45.0,r'$m = ' + str(m) +
↪ '$ kg')
[Link](0.7 * t_max,50.0,r'$\lambda = ' +
↪ str(lamb) + '$ kg/s')
[Link](loc=4)
[Link]()

Le résultat graphique est donné figure 2.2.

F IGURE 2.1

La mise en œuvre du schéma d’Euler implicite suit exacte-

ment la même procédure de codage. Ce travail est proposé
en fin de cette activité.

2.6 Solution exacte

L’équation (2.2) admet une solution analytique de la
forme :
v (t ) = g τ (1 − e −t /τ )

Son évolution temporelle peut être tracée en Python en

adoptant la même méthodologie que dans l’activité 1. Le
F IGURE 2.2
script suivant présente une solution dans laquelle la fonc-
tion v est à compléter.

import numpy as np
2.7 Solution Python
import [Link] as plt D’autres méthodes numériques, plus efficaces que les
schémas d’Euler, permettent d’approcher la solution d’une

Journée de formation - mars 2018 cbea

Physique numérique et Python 8

équation différentielle. C’est le cas de la commande # résolution numérique

odeint du module integrate du package scipy 1 . t_min = 0.0
odeint met en œuvre un mélange de méthodes nu- t_max = 5 * tau
mériques parmi lesquelles la méthode de Runge Kutta n_t = 100
d’ordre 4 et des méthodes à pas adaptatifs pour résoudre tab_t = [Link](t_min,t_max,n_t)
des équations différentielles du premier ordre. Dans la tab_v = odeint(eqDif,v0,tab_t, args=(g,h))
suite, nous n’utiliserons que cette commande du module
[Link]. C’est pourquoi nous pouvons charger La fin du script permet l’affichage de la courbe solution et
directement l’instruction pour éviter le recours à la nota- de quelques commentaires. Noter l’utilisation de la fonc-
tion pointée. tion max pour trouver la valeur maximale présente dans
le tableau tab_v. Cette valeur est utilisée pour placer les
from [Link] import odeint commentaires.
v_max = max(tab_v)
Cette commande reçoit au minimum trois arguments.
Le premier argument fait référence à l’équation différen- [Link](tab_t,tab_v,label="Méthode odeint")
[Link]('Temps (en s)')
tielle 2 . Le deuxième paramètre est un tableau de la ou des
[Link]('Vitesse de $M$ (en m/s)')
conditions initiales. Le troisième argument est le tableau
[Link]("Évolution de la vitesse")
des instants où la solution approchée est déterminée. Le ré- [Link](0.7 * t_max,0.8 * v_max,r'$h = ' +
sultat de l’appel à la fonction odeint est un tableau dont la ↪ str(h) + '$ m')
dimension dépend du nombre d’équations différentielles [Link](0.7 * t_max,0.7 * v_max,r'$t_{max} =
passées en argument. Un quatrième argument peut être ↪ ' + str(t_max) + '$ s')
ajouté qui indique à odeint les valeurs de paramètres pré- [Link](0.7 * t_max,0.6 * v_max,r'$m = ' +
sents dans les équations différentielles. ↪ str(m) + '$ kg')
Illustrons notre propos avec l’équation 2.2 rappelée ci- [Link](0.7 * t_max,0.5 * v_max,r'$\lambda =
dessous. ↪ ' + str(lamb) + '$ kg/s')
[Link](loc=4)
∀t ∈ R+
1 [Link]()
v̇ (t ) = g − v (t )
τ
Le résultat graphique est donné figure 2.3
Le second membre de cette équation permet de définir
une fonction qui est le premier argument de odeint. En
Python, cette fonction est définie avec au minimum deux
arguments : une variable 3 associée à la grandeur recher-
chée et le tableau des instants tab_t. Les éventuels autres
arguments sont les paramètres physiques de l’équation.

def eqDif(v,tab_t,g,tau):
return g - v / tau

La résolution numérique de l’équation différentielle pro-

cède de la même logique que celle présentée jusqu’ici avec
les autres méthodes numériques.

# paramètres physiques
g = 9.81 # accélération de la pesanteur
m = 1.0 # masse
F IGURE 2.3
lamb = 0.5 # coefficient de frottement
tau = m / lamb # temps caractéristique
h = 5.0 # hauteur initiale La méthode odeint fournit des résultats numériques tout
# vitesse initiale à fait pertinents. Elle est plus précise que les schémas d’Eu-
v0 = 0.0 ler dont l’objet était de présenter une technique numé-
rique simple, à la base de nombreuses autres méthodes.
Le script suivant précise les différentes étapes menant à la
résolution. Le résultat de l’appel à odeint est un tableau
2.8 À vous de jouer
qui contient autant de valeurs que le tableau tab_t. Ce ▶ Question 1. Mettre en œuvre le schéma d’Euler impli-
tableau est affecté à la variable tab_v. cite pour tracer l’évolution temporelle de la vitesse.
Sur un même graphique, tracer les courbes obtenues par
1. Le package scipy
2. Celle-ci peut d’ailleurs être vectorielle, ce qui permet la résolution de certaines équations différentielles d’ordre supérieur à 1. C’est souvent le
cas en physique.
3. Nous verrons plus loin qu’il peut s’agir d’un ensemble de variables, sous la forme d’un tableau.

Journée de formation - mars 2018 cbea

Physique numérique et Python 9

les deux schémas d’Euler et la courbe théorique. Observer

les différences en prenant un nombre de points de cal-
culs peu élevé. La figure 2.4 est une illustration du résultat
attendu.

F IGURE 2.4

▶ Question 2. La solution de l’équation 2.2 :

v (t ) = g τ (1 − e −t /τ )

permet de trouver l’évolution de z (t ). Attention à ne pas

oublier que v = −ż de sorte qu’après intégration :

∀t ∈ R+ , z (t ) = h + g τ2 (1 − − e −t /τ )
t
τ

On peut également voir z comme la solution de l’équation

différentielle suivante :

∀t ∈ R+ , ż (t ) = −g t − (z (t ) − h )
1
(2.4)
τ

avec la condition initiale z (0) = h.

▷ 2.1. Mettre en œuvre les schémas d’Euler, la méthode
odeint pour trouver des solutions approchées de cette
équation différentielle.
▷ 2.2. Tracer les courbes associées aux trois solution nu-
mériques et à la solution exacte.
▷ 2.3. Comparer les méthodes numériques pour de pe-
tites valeurs de n_t.
▷ 2.4. Observer le comportement des solutions sur des
intervalles de temps différents.
▶ Question 3. Les résultats précédents peuvent être ex-
ploités pour faire une étude énergétique numérique.
▷ 3.1. Comment définir deux tableaux tab_Ec et
tab_Ep contenant respectivement les valeurs numériques
des énergies cinétiques et potentielles ?
▷ 3.2. Tracer les évolutions temporelles des énergies ci-
nétiques et temporelles ainsi que celle de l’énergie méca-
nique. Commenter.

Journée de formation - mars 2018 cbea

Physique numérique et Python 10

Journée de formation - mars 2018 cbea

 ◾ Préciser l’intervalle de temps [t min , t max ] sur lequel
on désire tracer la courbe.

◾ Tracer la courbe et l’afficher.

3.4 Codage
Activité 3 Le début du script reprend en grande partie celui écrit dans
l’activité 1. Il est complété par la définition de la fonction
x et par la construction du tableau tab_x.
Trajectoire
import numpy as np
import [Link] as plt
# fonctions associée aux lois horaires
def x(t,g,v0,alpha):
3.1 Position du problème return v0 * [Link](alpha) * t
Les activités précédentes ont permis d’acquérir les mé- def z(t,g,h,v0,alpha):
thodes élémentaires de tracé de courbes et de résolution return h + v0 * [Link](alpha) * t - g *
numérique d’équations différentielles simples. Ces com- ↪ t**2 / 2
pétences peuvent être mises à profit pour étudier un mou- # paramètres physiques
g = 9.81 # accélération de la pesanteur
vement plan et tracer la trajectoire d’un corps à partir de
h = 5.0 # altitude initiale
la solution approchée d’un système différentiel.
v0 = 10.0 # vitesse initiale
La trajectoire plane du point matériel peut aisément être # angle en degrés
numériquement tracée dans le cas de la chute libre avec alpha_deg = 45.0
vitesse initiale, sans frottement. La prise en compte de # Python calcule en radians
frottements fluides est abordé sous forme d’exercice pour alpha = alpha_deg * [Link] / 180.0
compléter la prise en main de Python. # intervalle temporel d'étude
t_min = 0.0 # instant initial
3.2 Mise en équation t_max = 1.0 # instant final
n_t = 100 # nombre de points
Le point matériel M de masse m, soumis à la seule force
# tableaux
de pesanteur, est initialement lancé avec une vitesse v 0
tab_t = [Link](t_min,t_max,n_t)
faisant un angle α avec l’axe horizontal (Ox ), depuis le tab_x = x(tab_t,g,v0,alpha)
point de coordonnées (0, h ). Son vecteur vitesse initiale se tab_z = z(tab_t,g,h,v0,alpha)
décompose dans la base (ê x , ê z ) sous la forme :

v⃗(0) = v 0 cos αê x + v 0 sin αê z

Tracer la trajectoire se fait en utilisant la fonction plot
À un instant t ⩾ 0, les coordonnées du point M , notées avec comme arguments les deux tableaux tab_x et tab_z.
x (t ) et (z (t ), vérifient le système différentiel suivant.

ẍ (t ) = 0 [Link](tab_x,tab_z,'r-')
{ (3.1)
z̈ (t ) = −g

La solution de (3.1) est : La fin du script ajoute les commentaires.

x (t ) = v 0 t cos α
∀ t ∈ R+ { (3.2)
z (t ) = h + v 0 t sin α − 21 g t 2 [Link]('Abscisse (en s)')
[Link]('Cote (en m)')
[Link]("Chute libre avec vitesse
3.3 Méthodologie de tracé ↪ initiale")
Tracer la trajectoire revient à tracer z (t ) en fonction [Link](0.7 * t_max,35.0,r'$t_{max} = ' +
de x (t ) pour t ⩾ 0. En Python, cela revient d’abord à ↪ str(t_max) + '$ s')
construire les tableaux de valeurs des coordonnées x et z [Link](0.7 * t_max,40.0,r'$h = ' + str(h) +
pour t variant dans un intervalle fixé. La méthodologie est ↪ '$ m')
très proche de celle présentée dans l’activité. Il convient [Link](0.7 * t_max,45.0,r'$v_{0} = ' +
↪ str(v0) + '$ m/s')
simplement d’ajouter la construction du tableau associé
[Link](0.7 * t_max,50.0,r'$\alpha = ' +
aux valeurs de x.
↪ str(alpha_deg) + '$°')
◾ Définir les paramètres physiques du problème en [Link]()
précisant leurs valeurs numériques.
◾ Définir deux fonctionsx ∶ x ↦ v 0 t cos α et z ∶ t ↦ h +
v 0 t sin α − g t 2 /2. Le résultat graphique est donné figure 3.1

11
Physique numérique et Python 12

suivant :
⎧
⎪ ẋ (t ) = p (t )
⎪
⎪
⎪
⎪
⎪
⎪ (t ) = f x (x (t ), y (t ), z (t ), p (t ), q (t ), r (t ), t )
⎪
⎪
⎪
ṗ
⎪
⎪
⎪ ẏ (t ) = q (t )
∀ t ∈ R+ , ⎨
⎪
⎪
⎪ q̇ (t ) = f y (x (t ), y (t ), z (t ), p (t ), q (t ), r (t ), t )
⎪
⎪
⎪
⎪
⎪
⎪ ż (t ) = r (t )
⎪
⎪
⎪
⎪
⎩r˙(t ) = f z (x (t ), y (t ), z (t ), p (t ), q (t ), r (t ), t )
muni des conditions initiales :

x (0 ) = x 0 y (0 ) = y 0 z (0 ) = z 0
p (0) = v 0x q (0) = v 0y r (0) = v 0z

Sous cette forme, la commande odeint est à même de

F IGURE 3.1 donner une solution numérique du problème. La princi-
pale difficulté est de définir le système différentiel en vue
de son utilisation avec odeint. En pratique, ce système
3.5 Système différentiel et Python n’est rien d’autre qu’une relation vectorielle de la forme :
L’exemple de la chute libre est suffisamment simple pour
permettre l’écriture d’une solution analytique du système ⎛ x (t ) ⎞
différentiel. Mais de nombreuses situations physiques sont ⎜ p ( t )⎟
⎜ ⎟
modélisées par plusieurs équations différentielles cou- ⎜ y (t ) ⎟
Ẋ (t ) = F ( X (t ), t ) avec X (t ) = ⎜ ⎟
⎜ q ( t )⎟
plées qui peuvent souvent se mettre sous la forme sui- ⎜ ⎟
⎜ z (t ) ⎟
vante. ⎜ ⎟
⎝ r (t ) ⎠
⎧
⎪ ẍ (t ) = f x (x (t ), y (t ), z (t ), ẋ (t ), ẏ (t ), ż (t ), t )
⎪
⎪
⎪
+
∀t ∈ R , ⎨ ÿ (t ) = f y (x (t ), y (t ), z (t ), ẋ (t ), ẏ (t ), ż (t ), t ) Pour utiliser la fonction odeint, il suffit donc de définir
⎪
⎪
⎪
⎪
⎩z̈ (t ) = f z (x (t ), y (t ), z (t ), ẋ (t ), ẏ (t ), ż (t ), t )
la fonction vectorielle F . Le système suivant définit cette
fonction dans le cas de la chute libre sans frottement.
En physique, un tel système est complété par les condi-
tions initiales. ⎛ p ( t )⎞ ⎛ x (t ) ⎞
+ ⎜ 0 ⎟ ⎜ p ( t )⎟
∀ t ∈ R , F ( X ( t ), t ) = ⎜ ⎟ avec X (t ) = ⎜ ⎟
x (0 ) = x 0 y (0 ) = y 0 z (0) = z 0 ⎜ r (t ) ⎟ ⎜ z (t ) ⎟
⎝ −g ⎠ ⎝ r (t ) ⎠
ẋ (0) = v 0x ẏ (0) = v 0y ż (0) = v 0z
Le script suivant traduit cette formulation mathématique
Dans l’exemple de la chute libre, deux équations rendent
en code Python, en définissant la fonction F.
compte du phénomène physique. Les fonctions f x et f z
sont : def F(X,t,g):
x,p,z,r = X
f x (x (t ), y (t ), z (t ), ẋ (t ), ẏ (t ), ż (t ), t ) = 0 return [p,0,r,-g]
{
f z (x (t ), y (t ), z (t ), ẋ (t ), ẏ (t ), ż (t ), t ) = −g
Dans ce script, la fonction renvoie un tableau à deux di-
Leurs expressions sont très simples. Nous rencontrerons mensions. Chaque colonne du tableau contient les valeurs
dans une activité ultérieure une situation de couplage des numériques associées respectivement à x, p, z, r et com-
coordonnées en les deux équations. porte autant de valeurs que d’instants définis dans le ta-
En Python, un tel système différentiel n’est pas résolu nu- bleau tab_t. Python offre un moyen simple de récupérer
mériquement directement sous cette forme. Le système chacune des colonnes pour les affecter à des variables
faisant intervenir des dérivées d’ordre supérieur à 1 doit idoines. Chaque colonne porte un numéro qui commence
d’abord être transformé en un système du premier ordre 1 . à 0. Ainsi, à la colonne 0 est associé le tableau des valeurs
Cette transformation est simple à réaliser en introduisant de x, à la colonne 1 celui des valeurs de p, à la colonne 2
de nouvelles fonctions 2 qui « absorbent » une partie des celui des valeurs de z et à la colonne 3 celui des valeurs
dérivées. Posons : de r . L’opération de sélection d’une colonne en Python
est appelée slicing, littéralement tranchage. Le script sui-
p = ẋ q = ẏ r = ż
vant illustre cette opération de récupération des tableaux
de valeurs de x et de z, après avoir utilisé la commande
Ainsi, le problème physique peut être formellement décrit
par le système différentiel équivalent du premier ordre
odeint.

1. Cette analyse, a priori mathématique, est à rapprocher des modalités d’étude des mouvements en mécanique analytique et de la notion
d’espace de phases.
2. auxquelles il est souvent possible de donner un sens physique

Journée de formation - mars 2018 cbea

Physique numérique et Python 13

# paramètres physiques
g = 9.81 # accélération de la pesanteur
x0 = 0.0
z0 = 5.0 # altitude initiale
v0 = 10.0 # vitesse initiale
# angle en degrés

alpha_deg = 45.0
# Python calcule en radians
alpha = alpha_deg * [Link] / 180.0
cond_init = [x0,v0 * [Link](alpha),z0,v0 *
↪ [Link](alpha)]
# intervalle temporel d'étude
t_min = 0.0 # instant initial
t_max = 1.0 # instant final
n_t = 3 # nombre de points
# tableaux
tab_t = [Link](t_min,t_max,n_t)
sol_num = odeint(F,cond_init,tab_t,args=(g,))
tab_x = sol_num[:,0] # tableaux des x
tab_z = sol_num[:,2] # tableaux des z

Le tracé avec les options est immédiat. La figure 3.2 illustre

la mise en œuvre du script précédent.

F IGURE 3.2

3.6 À vous de jouer

▶ Question 1. En vous aidant du dernier script, tracer
l’hodographe du corps. Le résultat est a priori évident.
▶ Question 2. Comment modifier le script précédent en
présence de frottements fluides ?

Journée de formation - mars 2018 cbea

Physique numérique et Python 14

Journée de formation - mars 2018 cbea

Activité 4

Oscillateurs

4.1 Position du problème F IGURE 4.1

Nous sommes à présent en mesure d’étudier des situa-
tions physiques plus complexes que la chute libre. L’étude
des oscillations occupe une place importante en physique.
Cette activité vise à réinvestir les compétences acquises ▶ Question 2. Tracer le portrait de phase de l’oscillateur
dans les précédentes activités en vue de tracer des lois en complétant le code précédent.
horaires et des portraits de phase. Une ouverture vers les
oscillations non-linéaires permet de dépasser le cadre des
situations pour lesquelles une solution analytique existe

4.2 Notations
On s’intéresse tout d’abord au mouvement à un seul degré
de liberté d’un point matériel soumis à une force de rappel
élastique et à une force de frottements fluides. Le corps
n’est soumis à aucune autre force. Cette situation peut être
modélisée par l’équation différentielle du second ordre
suivante.

∀t ∈ R+ , ẍ (t ) + 2ζω0 ẋ (t ) + ω20 x (t ) = 0 (4.1)

Dans cette relation, ω0 désigne une pulsation caractéris-

tique et ζ est une quantité positive sans dimension, appe-
lée taux d’amortissement.
Les conditions initiales sont notées :

x (0 ) = x 0 ẋ (0) = v 0
F IGURE 4.2
En posant p (t ) = ẋ (t ), l’équation (4.1) est équivalente au
système différentiel du premier ordre suivant :

ẋ (t ) = p (t )
∀t ∈ R+ , { ▶ Question 3. Pour disposer de plusieurs tracés dans
ṗ (t ) = −2ζω0 p (t ) + ω20 x (t )
une même fenêtre graphique, Python dispose de la com-
muni des conditions initiales : mande subplot. Le script suivant illustre son utilisation.

x (0 ) = x 0 p (0 ) = v 0
import numpy as np
L’objet de cette activité est de résoudre numériquement ce import [Link] as plt
système puis d’étudier des systèmes oscillants plus com- # fonction à tracer
plexes. L’outil numérique se révèlera ici particulièrement def f(x,k):
efficace pour atteindre ces objectifs. return [Link](k * x)
def g(x,k):
4.3 À vous de jouer return [Link](k * x)
# données de tracé
▶ Question 1. Proposer un script complet qui donne
x_min = 0.0
l’évolution temporelle de x pour différentes valeurs de x_max = 2 * [Link]
ζ. Ces valeurs seront choisies de manière à mettre en évi- n_x = 100
dence différents comportements oscillatoires. Les courbes tab_x = [Link](x_min,x_max,n_x)
seront tracées sur un même graphique.

15
Physique numérique et Python 16

# liste de valeurs de k sinusoïdale. Illustrer le phénomène de résonance.

lst_k = [1,2,3] ▶ Question 6. L’oscillateur de Van der Pol peut être mo-
# boucle de tracés délisé par une équation différentielle du second ordre de
# pour différentes valeurs de k la forme :
for k in lst_k:
[Link](2,1,1) ∀t ∈ R+ , ẍ (t ) − εω0 [1 − x 2 (t )] ẋ (t ) + ω20 x (t ) = 0
y = f(tab_x,k)
[Link](tab_x,y,label="$k = " + str(k) Illustrer son comportement en fonction des paramètres du
↪ +"$") problème. En particulier, mettre en évidence l’existence
[Link](2,1,2) d’un cycle limite indépendant des conditions initiales.
y = g(tab_x,k)
[Link](tab_x,y,label="$k = " + str(k)
↪ +"$")
[Link]()

La figure 4.3 illustre l’exécution de ce script.

F IGURE 4.3

Proposer un script qui affiche dans une même fenêtre,

l’évolution temporelle de x, le portrait de phase et, sur
un même graphique, l’évolution temporelle des énergies
cinétique, potentielle et mécanique.
▶ Question 4. Proposer une modélisation des oscilla-
tions d’un pendule. Exploiter toutes vos compétences nu-
mériques pour présenter une solution numérique du pro-
blème.

F IGURE 4.4

▶ Question 5. Le système est forcé par une excitation

Journée de formation - mars 2018 cbea

 y

M v⃗0
r α
Activité 5 θ
x
O M0

Canon de newton

5.1 Position du problème F IGURE 5.2 – Notations

On attribue à Newton une expérience de pensée dont l’un
des objectifs premiers était de montrer que la gravitation ◾ Le mouvement de M est plan. Le plan du mouve-
constituait la composante majeure de la force de pesan- ment est orthogonal au vecteur moment cinétique
teur. Dans cette expérience, Newton place un canon au constant L⃗ et passe par O. Par choix, ce plan consti-
sommet d’une « très haute » montagne. Ce canon tire des tue le plan (xO y ), l’axe (Oz ) étant défini comme
boulets qui, suivant leurs vitesses initiales, soit tombent ⃗
l’axe porté par L.
sur la Terre, soit se satellisent, soit partent dans l’espace. ◾ Dans le plan du mouvement, le point M peut être
repéré par ses coordonnées polaires (r, θ ). Celles-ci
vérifient les équations dynamiques suivantes.

L2 GM L
r¨ = − θ̇ =
m2r 3 r2 mr 2

où G désigne la constante de gravitation universelle.

◾ Ce sont ces équations dont il convient de donner
une solution numérique afin de construire la trajec-
toire du point M .

5.3 À vous de jouer

La Terre est assimilée à une sphère de rayon R, de masse
M . Son centre définit l’origine d’un référentiel supposé
galiléen. En outre, on suppose la Terre immobile.
Un canon est placé à une altitude h. Des boulets sont ti-
F IGURE 5.1 – Retombée sur Terre (A, B) - Satellisation (C,
rés avec un vecteur vitesse initial v⃗0 dont la norme et la
D) - Vers l’infini et au-delà (E)
direction peuvent être modifiées à souhait. Chaque boulet
a une masse m et est soumis à la seule force de gravitation
exercée par la Terre.
L’objet de cette activité est de simuler numériquement
Le travail demandé pourra se décomposer de la manière
cette expérience de pensée. Il sera l’occasion de confron-
suivante.
ter les vitesses de mise en orbite d’un corps. De nom-
breux sites web proposent une illustration de cette ▶ Question 1. Adimensionner les équations du mouve-
expérience de pensée. Vous pouvez par exemple tes- ment.
ter le site suivant : [Link]
▶ Question 2. Définir le système différentiel et les
physicsflash/home/gravity.
conditions initiales associées au problème adimensionné.

5.2 Mise en équation ▶ Question 3. Résoudre numériquement le système.

Le problème précédent est une application directe des ▶ Question 4. Tracer la planète et les trajectoires des
équations régissant le mouvement d’un corps dans le boulets.
champ gravitationnel d’un corps massif immobile, en l’oc- ▶ Question 5. Déterminer les conditions sur v⃗0 (direc-
curence la Terre dans le cas présent. En désignant par M tion et intensité) pour observer la retombée du boulet sur
la position du centre de masse du boulet et par O le centre la Terre ou sa mise en orbite. Comparer les valeurs aux
de la Terre, les lois de la mécanique permettent d’établir valeurs connues sur Terre (vitesse de satellisation, vitesse
les résultats suivants. de libération).

17
Physique numérique et Python 18

5.4 Données numériques

Masse de la Terre
M = 5,98 × 1024 kg
Rayon terrestre
R = 6370 km
Constante de gravitation universelle
G = 6,67 × 10−11 m3 ⋅ kg−1 ⋅ s−2
Vitesse de satellisation minimale — La vitesse de satellisa-
tion minimale est la vitesse théoriquement communiquée à
un corps au départ d’un astre pour le satelliser au plus près
de ce dernier sur une orbite circulaire.
v s = 7,9 km ⋅ s−1
Vitesse de libération — La vitesse de libération est la vitesse
minimale que doit atteindre un corps pour échapper défini-
tivement à l’attraction gravitationnelle d’un astre (planète,
étoile, etc.) et s’en éloigner indéfiniment.
v s = 11,2 km ⋅ s−1
Vitesse de la lumière dans le vide
c = 2,998 × 108 m ⋅ s−1

Journée de formation - mars 2018 cbea

 import numpy as np
import [Link] as plt

def f(x):
return [Link](-x * x)

Activité 6 x_min, x_max = -4.0, 5.0

n_x = 100
x = [Link](x_min,x_max,n_x)
y = f(x)

Petit vade-mecum [Link](x,y)

[Link]()

6.1 Tracé d’une courbe à partir 1.0

d’une fonction
Tracé de la courbe représentative de la fonction f ∶ x ↦ 0.8
−x 2
e sur un intervalle [x mi n , x max ]. Par défaut, Python
trace une courbe en reliant les points calculés. Noter l’ap- 0.6
pel aux deux modules numpy et [Link].

import numpy as np 0.4

import [Link] as plt
0.2
def f(x):
return [Link](-x * x)
0.0
4 3 2 1 0 1 2 3 4 5
x_min, x_max = -4.0, 5.0
n_x = 100
x = [Link](x_min,x_max,n_x) F IGURE 6.2
y = f(x)

[Link](x,y)
[Link]()
6.2 Tracé d’une courbe à partir de
données tabulaires
En physique, les données numériques sont d’origine ex-
1.0
périmentales. Obtenues à partir de mesures (acquisitions
de données via un CAN), elles peuvent être enregistrées
dans un fichier. Python peut lire les données pour ensuite
0.8 afficher une courbe.
Le format .csv est fréquemment adopté pour enregistrer
0.6 des données ou les exporter depuis un logiciel. En pratique,
les données sont stockées dans un fichier sous forme de
flottants, séparés par un point virgule. Le module csv de
0.4
Python fournit quelques fonctions pratiques pour mani-
puler les fichiers .csv.
0.2 Dans l’exemple suivant, le fichier de données [Link]
contient des lignes de deux flottants séparés par un point-
0.0 virgule.
4 3 2 1 0 1 2 3 4 5

import [Link] as plt

F IGURE 6.1 import csv

source = open('[Link]', 'r')

Pour faire apparaître les points, on peut ajouter une option
à [Link]. x, y = [], []
◾ 'b-' trace la courbe en bleu (option b) et relie les for row in [Link](source,delimiter=';'):
x1, y1 = map(float,row)
points (option -).
[Link](x1)
◾ 'ro' affiche des points (option o) en rouge (option [Link](y1)
r).

19
Physique numérique et Python 20

[Link](xFit,yFit,label='fit',color='r')
[Link](x,y,'b-') [Link](x,y,'bo')
[Link](x,y 'ro') [Link]()
[Link]()

3.0
3.0
2.5

2.5
2.0

2.0 1.5

1.5 1.0

1.0 0.5

0.5 0.0
0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4

0.0
0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 F IGURE 6.4

F IGURE 6.3

6.3 Ajustements
Il est parfois nécessaire de déterminer un modèle mathé-
matique rend compte de l’évolution d’un phénomène phy-
sique. Par exemple, une série de mesures a fourni un en-
semble de couples de points (x i , y i ) qui semblent s’aligner.
Comment déterminer l’équation de la droite de régres-
sion associée ? Python dispose d’une fonction curve_fit
du module from [Link] qui permet tout type
d’ajustement : ajustement linéaire, ajustement gaussien,
etc. L’exemple suivant présente un ajustement gaussien de
données. Il reprend une partie du code vu ci-dessus. Noter
la conversion des données en tableau (fonction ar).

import [Link] as plt

import csv
import numpy as np
from [Link] import curve_fit
from scipy import asarray as ar

def gauss(x,a,x0,sigma):
return a*[Link](- (x - x0)**2 / sigma**2
↪ / 2)

source = open('[Link]','r')
x, y = [], []
for row in [Link](source,delimiter=';'):
x1, y1 = map(float,row)
[Link](x1)
[Link](y1)

x, y = ar(x), ar(y)
xmin, xmax, n_x = min(x), max(x), 100

popt, pcov = curve_fit(gauss,x,y)

xFit = [Link](xmin,xmax,n_x)
yFit = gauss(xFit,*popt)

Journée de formation - mars 2018 cbea