Introduction à la simulation

Etude de la chute d'un corps - Une introduction à la simulation

Introduction

La simulation est un outil puissant pour étudier et comprendre les phénomènes physiques. Cependant, la plupart des programmes de simulation, y compris ceux présentés sur ce site, font appel à des outils mathématiques qui sont à la limite du cursus du secondaire. Par exemple, les équations différentielles sont abordées en terminale S. Et pourtant, la physique enseignée dès la classe de seconde décrit des phénomènes intéressants, en mécanique en particulier...
Comment faire! C'est le "défi" que m'a récemment lancé, sans le savoir, un élève de seconde , passionné par la physique. Il voulait simuler le mouvement d'une balle lancée avec une vitesse initiale non nulle, en négligeant les frottements et autres désagréments, bien sur. Mon premier réflexe fut de lui indiquer l'EDO qui va bien... Las, les EDO (équations différentielles ordinaires) ne sont pas au programme de seconde!
Et puis, comment aborder le concept de simulation avec un élève de seconde! C'est à cette question que tente de répondre cette page.
Je suis preneur de tous les commentaires, en particulier ceux des élèves de seconde et de premières, et leurs professeurs. Il m'intéresse de savoir si la méthode que j'utilise est compréhensible au niveau de la seconde. Car hélas, je suis comme Groucho Marx, je n'ai pas d'élève de seconde sous la main...

Autre chose, les habitués de TangenteX savent que je programme habituellement en FORTRAN ou en C. Ici, j'ai changé. J'ai décidé d'utiliser Scilab. Pour une raison essentielle : sa riche librairie graphique, qui permet de faire des tracés de courbes sans se poser de question! Sa syntaxe est simple et cerise sur le gâteau, c'est un logiciel libre français (INRIA/ENPC). J'aurais pu utiliser Matlab, mais lorsque vous verrez la facture, vous comprendrez...

Comment construire une simulation

Simuler un phénomène physique, c'est construire un modèle dont le comportement sera le plus proche possible de celui du phénomène. Ce modèle produit des résultats numériques qui doivent être comparables avec les données réelles d'évolution du système.
Par exemple, simuler la chute d'un corps c'est écrire un modèle qui calcule la vitesse et la position par rapport au temps. Ces informations doivent être les plus proches possibles de celles que l'on pourrait mesurer en faisant l'expérience réellement.
Pour construire une simulation, nous devons suivre quatre étapes:

Modéliser:

C'est la partie physicienne du boulot. Il s'agit de déterminer les lois de la physique qui sont applicables et de définir leur domaine d'application. Cela peut être très simple :définir un référentiel galiléen et établir l'équation de mouvement d'un corps dans ce référentiel. Cela peut être très compliqué, comme par exemple étudier la turbulence dans différentes parties du domaine de vol d'un avion.
A l'issue de cette étape, j'ai obtenu un système d'équations et de conditions initiales et aux limites qui constituent le modèle physique de mon phénomène.

Elaborer l'algorithme:

C'est la partie mathématicienne! Il s'agit de trouver le ou les bons schémas numériques pour résoudre nos équations dans leur domaine d'application. Et là, ça peut être du sport. Souvent, on peut utiliser des schémas standards: tout va bien! Mais dans certains cas, ces schémas divergent ou ne répondent pas au besoin: il faut se débrouiller ou bien faire appel au mathématicien de service! Après tout il est là pour ça...
Après avoir trouvé le bon schéma, il faut s'assurer qu'il est optimisé et performant. Si votre algorithme met 100 ans à produire un résultat (cela pourrait arriver...), cela ne vous avancera à rien! Là encore, il faudra ressortir le cours d'analyse numérique ou d'analyse tout court!
Bref, une fois cette cuisine accomplie, je suis en possession d'un splendide algorithme, c'est à dire une liste d'opérations mathématiques à exécuter, qui en principe me donneront au final le résultat attendu...

Ecrire le programme:

Et enfin la partie informaticienne, car vous l'avez compris, la physique numérique est à la croisée de trois disciplines : la physique, les mathématiques et l'informatique. C'est ce qui en fait tout son intérêt, à mes yeux du moins!
Il s'agit d'une autre cuisine: muni de mon algorithme et d'un langage informatique, j'écris un programme qui implémente mon algorithme. Cela signifie que je transforme ma suite d'opérations mathématiques en une suite d'instructions informatiques.
Savez-vous que FORTRAN (FORmula TRANslator) a été créé pour cela spécialement! Mais on peut aussi utiliser C, C++ ou Python (laissez tomber java, php ou autres langages du genre! Ils ne sont pas conçus pour faire du calcul, n'en déplaise à certains...)
On peut aussi se tourner vers des langages de calcul comme Matlab ou, en logiciels libres, Octave ou Scilab. C'est ce que j'ai fait ici, en faisant des infidélités à mon FORTRAN favori.
Comme ici je ne fais pas cours d'informatique, je vous laisse, si besoin, vous former ailleurs à cette science admirable, bien que vous trouverez ailleurs sur mon site une initiation au FORTRAN et au C.
Me voici donc en possession d'un splendide programme, conforme à mon algorithme et débuggé....

Exploiter:

Revenons à la physique! Il s'agit maintenant, en utilisant notre programme, d'observer et d'étudier le comportement de notre système physique à travers le programme de simulation.
Il est en effet possible de faire varier toutes sortes de paramètres (si vous l'avez prévu...) pour observer les réactions du système. On peut faire avec un programme des choses qui ne sont pas faisables dans un labo (faire décoller une fusée, former une galaxie, mettre le feu dans un tunnel, créer un cyclone, etc.). C'est toute la puissance de la simulation!
Un ordinateur, un compilateur, des algos et un manuel d'analyse numérique : voilà les seuls outils nécessaires pour faire de la physique numérique: Ah si, j'oubliais, il faut un cerveau aussi...
Essayons maintenant d'appliquer ces principes très théoriques...

Analyser le cas de la chute d'un corps

Quelques généralités d'abord....

Nous voulons simuler le mouvement d'un corps solide, une balle de tennis par exemple, dans le champ de gravité terrestre. En d'autres termes, voir comment se comporte une balle qui tombe! Pour corser le problème et le rendre plus réaliste, nous allons supposer que la balle est lancée horizontalement, avec une certaine vitesse. On ne s'occupera pas de détails comme la résistance de l'air, la poussée d'Archimède ou autres joyeusetés du genre...
Un corps en mouvement dans l'espace change de position au cours du temps. On appelle trajectoire le lieu de l'ensemble des positions dans l'espace en fonction du temps.
Pour simuler le mouvement de la balle de tennis, nous allons donc construire sa trajectoire. Et pour ce faire, nous devons déterminer deux objets:

un repère, qui nous permettra de noter la position de la balle par ses coordonnées. En principe, la balle bouge dans un espace à 3 dimensions! Mais pour simplifier, on se contentera des 2 dimensions du plan. Notre repère sera donc le repère Oxy que vous connaissez.
une horloge. Pour calculer chaque point de la trajectoire, il faut mesurer le temps qui passe! Et donc posséder une horloge que l'on déclenchera au début du mouvement.

Facile, jusque là, n'est-ce pas!
Mais on a un problème! Dans une simulation, le temps et l'espace ne sont pas continus. Cela veut dire qu'entre l'instant t1 et l'instant t2, il n'y a RIEN. De même, entre l'abscisse x1 et l'abscisse x2, il n'y a rien non plus!
Pour faire savant, on dit que l'on "discrétise" le temps et l'espace: ils ne sont plus continus, ils sont discrets! Vous verrez plus tard que c'est un aspect très important de la simulation numérique...
En ayant présent à l'esprit ces quelques points, voyons comment construire cette simulation.

Le modèle

Définissons d'abord un repère orthonormé, muni de deux vecteurs unitaires u et v (selon la convention en vigueur chez les physiciens, j'écris les vecteurs en gras).
A l'origine du mouvement, les coordonnées de la balle sont (x0,y0). x0 vaut 0 et y0 est un paramètre du modèle : je pourrais le faire varier comme je veux.
Je lance la balle à une vitesse v₀,horizontalement.

Je vais maintenant observer ce qui se passe peu de temps après l'instant où j'ai lancé la balle (le t₀). C'est là qu'intervient notre horloge. Je vais fixer un intervalle de temps petit, que j'appelle deltat, et je vais observer ce qui se passe à t₀+deltat.

Les nouvelles coordonnées de la balle sont (x,y). Et notre problème est de déterminer, à l'instant t₀+dt, les valeurs de x et y.

Faisons le bilan des forces qui s'exercent sur la balle:

1 - son poids.

2 - la vitesse horizontale. Elle ne varie pas.

La balle subit une accélération qu'il faut décomposer en une:

accélération verticale, qui est égale à g. Vous savez en effet que le poids est proportionnel à g, selon la célèbre formule p = mg
accélération horizontale, qui est nulle. En effet, la vitesse horizontale est constante, ce qui implique d'après la première loi de Newton, que l'accélération est nulle.

Voilà le bilan. Mais quoi faire de ces données pour évaluer le mouvement de la balle dans un court intervalle de temps?
Il faut se souvenir de plusieurs choses, qui découlent des définitions de l'accélération et de la vitesse:

sur un court intervalle de temps, l'augmentation (ou la diminution) de la vitesse est égale au produit de l'accélération par l'intervalle de temps. Par exemple, une voiture qui accélère à 5 m.s-1 par seconde voit sa vitesse augmenter de 5 m.s-1 au bout d'une seconde.
sur un court intervalle de temps, la distance parcourue est égale au produit de la vitesse par le temps. Par exemple, la même voiture qui roule à 30 m.s-1 parcoure 30 m en 1 seconde.

Essayons de traduire ces informations sous forme d'équations.

Je sais que la vitesse horizontale est constante et égale à v₀.
Si j'appelle v_i la vitesse verticale de la balle à l'instant t_i, et v_i+1 la vitesse verticale de la balle à l'instant t_i+1, je peux écrire v_i+1 - v_i = -g*deltat, deltat étant l'intervalle de temps (0,1 seconde par exemple). Vous noterez le signe moins qui affecte g. Vous avez remarqué que le poids est orienté vers le bas, et donc que g est aussi orienté vers le bas : la balle tombe, n'est-ce pas! Plus mathématiquement, le vecteur p s'écrit -mgv, puisque les vecteurs p et v sont de sens opposé. Et donc la norme du vecteur p est égale à -mg.
Si j'appelle x_i l'abscisse de la balle à l'instant t_i, et x_i+1 l'abscisse de la balle à l'instant t_i+1, je peux écrire x_i+1 - x_i = v₀*deltat,
Si j'appelle y_i l'ordonnée de la balle à l'instant t_i, et y_i+1 l'ordonnée de la balle à l'instant t_i+1, je peux écrire y_i+1 - y_i = v*deltat, où v est la vitesse verticale de la balle.

Voilà tous les éléments pour établir l'algorithme de simulation.

L'algorithme

Je veux calculer la trajectoire de la balle depuis sa position d'origine jusqu'à ce qu'elle atteigne le sol.
Pour cela, je vais calculer les coordonnées (x,y) de chaque point de la trajectoire de proche en proche en partant de l'origine (x₀,y₀) à t₀, par petits intervalles de temps, jusqu'à ce que l'abscisse de la balle soit nulle.
Voyons le structure globale de l'algorithme, qui reprend toutes les données exposées ci-dessus:

x(1) = 0

y(1) = y₀

VitesseVerticale(1) = 0

TantQue y(i-1) > 0 Faire

VitesseVerticale(i) = VitesseVerticale(i-1) - g*deltat

x(i) = x(i-1) + v₀*deltat

y(i) = y(i-1) + VitesseVerticale(i)*deltat

Plotter (x,y)

FinFaire

Cet algorithme, même simple, mérite quelques explications:

Chaque fois qu'une donnée est calculée, elle est conservée dans un tableau. C'est le cas de la vitesse verticale, qui évolue dans le temps ainsi que les coordonnées x et y. La raison de ce stockage est simple. Au temps i, la vitesse et la position sont calculées à partir des données calculées au temps i-1. Les conditions initiales sont stockées dans la premier case du tableau correspondant.
La boucle TantQue s'achève lorsque y devient nul ou négatif.
Plotter(x,y) signifie tracer le point (x,y) dans le repère.
Pour la petite histoire, les habitués auront reconnu l'utilisation de la méthode des différences finies, au premier ordre (ou méthode d'Euler). Cette méthode donne des résultats acceptables si le pas de temps est petit. On le verra en faisant les manip. plus loin.

le programme

Voici le petit programme Scilab qui implémente l'algorithme ci-dessus.

// *********************************

// Etude de la chute d'un corps

// Programme de mécanique de seconde

// Dominique LEFEBVRE

// 26 février 2007

// *********************************

clear; // Initialisation de toutes les variables

// Initialisation des paramètres de simulation

deltaT = 0.01; // base de temps = 0,01 seconde

// Initialsation des paramètres du mouvement

g = 9.81; // accélération de la pesanteur 9,81 N.m-2 (ou N.kg)

v0 = 2 ; // vitesse initiale du corps (en m.s-1)

y0 = 10; // ordonnée initiale du corps (en m)

// Initialisation des variables de calcul

x(1) = 0; // abscisse initiale du corps = 0;

y(1) = y0; // ordonnée initiale = y0, paramètre de la simulation

vX = v0; // vitesse horizontale initiale du corps - cette vitesse est constante

vY(1) = 0; // vitesse verticale initiale du corps

// Initialisation graphique

xbasc();

// Calcul du mouvement

i = 2; // indice de boucle

while (y(i-1) > 0) // la boucle est exécutée tant que le corps n'est pas sur le sol

vY(i) = vY(i-1) - g*deltaT;

x(i) = x(i-1) + vX*deltaT;

y(i) = y(i-1) + vY(i)*deltaT;

i = i + 1; // on passe à l'indice de tableau suivant

end

// Affichage du mouvement dans le repère choisi

plot2d(x,y,rect=[0,0,3,y0],style = 5);

xtitle('Chute d''un corps','X','Y');

Pour charger le programme ChuteCorpsXY.

Quelques expériences

Nous pouvons maintenant passer à la dernière étape, la plus intéressante : l'exploitation du programme. J'ai donc fait tourner le programme avec différentes valeurs de v₀ et de pas de temps pour apprécier la précision de mon modèle.
Nous avons ici la chance de connaitre la solution exacte du calcul. La portée est donnée par x_max = v₀*sqrt(2*y₀/g) où sqrt signifie "racine carrée" ou "square root".

Tracer la courbe y =f(x)

Dans le programme, je fixe la valeur de deltat = 0.1s et v₀= 2 m.s-1. Voilà la courbe calculée par le programme:

J'obtiens une portée de 2,8m. Pour information, dans les conditions initiales du calcul, la valeur exacte de la portée est 2,856 m.
Pour deltat = 0.01s v₀= 2 m.s-1

Comme vous le constatez, avec un pas de temps de 0,01s, on obtient 2,86 m: pratiquement le résultat exact!
Pour deltat = 0.1s v₀= 10 m.s-1

Pour information, dans les conditions initiales du calcul, la porté exacte est 14,28 m.
Deltat = 0.01s v₀= 10 m.s-1

On est très près là aussi du résultat exact.

Tracer la courbe z=f(t)

Autre facette du problème: on peut chercher à visualiser l'altitude de la balle en fonction du temps! Il suffit d'adapter un peu le programme, qui devient:

// *********************************

// Etude de la chute d'un corps

// Tracé de la courbe z = f(t)

// Programme de mécanique de seconde

// Dominique LEFEBVRE

// 26 février 2007

// *********************************

clear; // Initialisation de toutes les variables

// Initialisation des paramètres de simulation

deltaT = 0.01; // base de temps = 0,1 seconde

// Initialsation des paramètres du mouvement

g = 9.81; // accélération de la pesanteur 9,81 N.m-2 (ou N.kg)

v0 = 2 ; // vitesse initiale du corps (en m.s-1)

z0 = 10; // ordonnée initiale du corps (en m)

// Initialisation des variables de calcul

z(1) = z0; // ordonnée initiale = y0, paramètre de la simulation

vZ(1) = 0; // vitesse verticale initiale du corps

t(1) = 0;

// Initialisation graphique

xbasc();

// Calcul du mouvement

i = 2; // indice de boucle

while (z(i-1) > 0) // la boucle est exécutée tant que le corps n'est pas sur le sol

vZ(i) = vZ(i-1) - g*deltaT;

z(i) = z(i-1) + vZ(i)*deltaT;

t(i) = t(i-1) + deltaT;

i = i + 1; // on passe à l'indice de tableau suivant

end

// Affichage du mouvement dans le repère choisi

plot2d(t,z,rect=[0,0,3,z0],style = 5);

xtitle('Chute d''un corps - Pas de temp: 0.01 ','Temps','Altitude');

Pour charger le programme ChuteCorpsZ.
Si je le fait tourner avec les conditions initiales deltat = 0,01 s et v₀= 2m.s-1, j'obtient la courbe:

Le temps exact est donné par t = sqrt(2*y₀/g) soit 1,428 s. Le programme nous le donne à 1,42 s. Pas mal, vu les moyens mathématiques mis en jeux...

Quelques remarques physiciennes:

1 - le mouvement est uniformément accéléré: on constate que la vitesse de la balle augmente à chaque intervalle de temps.

2 - le temps de chute ne dépend pas de la vitesse initiale horizontale. Pour vous en convaincre, regardez le calcul de z en fonction de t: la vitesse horizontale n'intervient pas.

3 - le temps de chute ne dépend pas de la masse de la balle! Elle n'intervient pas dans le calcul: regardez le programme... Galilée avait donc raison : si l'on néglige les frottements et autres (dans le vide, par exemple) la bille de plomb et la plume atteignent le sol simultanément.

4 - la portée dépend uniquement de la vitesse initiale. Faites varier v₀ dans le programme ChuteCorpsXY, vous verrez varier la portée du tir. Là encore, la masse de la balle n'intervient pas!

5 - la forme de la courbe dépend de la vitesse initiale. Si elle est nulle, la courbe est une droite (qui se confond d'ailleurs avec l'axe Oy car x₀ = 0). Si elle est non nulle, c'est une parabole.