PHYS710 : Simulation et modélisation en physique

Le but de ce cours est d'initier les étudiants au développement d'un logiciel de simulation ou de modélisation en physique. Ce logiciel sera modulaire et composé de sous-ensembles interagissant entre eux. Chaque module ou sous-ensemble sera développé par un binôme qui en aura la responsabilité. Le but sous-jacent est l'initiation à un travail collaboratif.

Dans cet esprit, la première séance consistera en une introduction succinte à l'environnement de développement et au langage choisis, en l'occurence ROOT et le C/C++. Les séances suivantes, espacées de 2 à 3 semaines seront des occasions de faire le point, chaque binôme expliquant aux autres l'évolution de leur travail et les points d'accès pour que tous puissent faire des tests entre les modules. Entre temps, les étudiants communiqueront avec l'enseignant et leurs camarades à travers un ensemble d'outils tels cette page wiki et une liste de diffusion. L'enseignant collectera les contributions (et suggérera/demandera éventuellement des modifications) de telle sorte qu'elles soient disponibles pour tous.

Autant que faire se peut, on envisagera une utilisation ultérieure du logiciel, soit dans le cadre des exercices "WIMS" de première année, soit pour servir de point de départ aux projets des années suivantes.

Introduction, outils de développement

Cadre général

• Simulation, modélisation en physique
• Développement d'un logiciel de modélisation
• Contraintes, travail en groupe
• Exemple: un système mécanique complexe, le moteur à explosion

Environnement

On utilisera ROOT.

• Rappel de quelques commandes Unix (si nécessaire)
• Utilisation de ROOT
• Compilation d'un petit code C
• Compilation d'un code intégrant des fonctionnalités ROOT

Langage

Nous utiliserons le C/C++. Les notions générales de C seront supposées connues (petit rappel rapide). On introduira les notions suivantes avec à chaque fois l'analyse de petits exemples de macros ROOT.

• Rappels de C : variables, fonctions, structures de contrôle, exemple de programme complet
• Pointeurs
• Le préprocesseur
• Allocation dynamique en C
• Le langage C++
• Classes, objets
  • Héritage
  • Constructeurs, destructeurs
  • Encapsulation
• Allocation dynamique en C++

Petit programme complet à la fin de cette page. Essayez-le dans ROOT.

En règle générale, on donnera les informations nécessaires pour comprendre les macros ROOT présentées. Il reste beaucoup d'auto-formation à faire de la part des étudiants !

Quelques références

• cours de C++ à partir du C.
• "mégacours" de C/C++.
• Cours sur "le site du zéro"
• Ce qu'il ne faut pas faire : les anti-patrons

Bibliothèques

gsl, modules ROOT, bibliothèques trouvées ou produites par les étudiants

Résumé du 1^er cours : rappel sur le langage C

Variables

int (entier)

float (réel 4 octets)

double (réel 8 octets)

char (caractère 1 octet)

...

Déclaration des variables

=> int i,j,k;
Pour déclarer un tableau qui va de 0 à 99 : float a[100]; <br\>Puis initialisation : a[i]=0;

Boucles

=> for, while, do...
Exemple avec la boucle for : for (i=1 ; i<36 ; i++) {......}

Tests

if (i<35 && i>27 && i!=30) {......} else if (...) {...}
On peut avoir :

&& "et"

== "égal"

|| "ou"

!= "différent"

Commentaires

en C : /*....*/

en C++ : //

Pointeurs

Déclaration : float* a;
Utilisation : (*a)=18;
Exemple :

z=23;

a=&z;

*a=18;

On obtient donc au final que z=18

Place réservée

à venir

Outils de travail collaboratif

Ce wiki : Pour réaliser des pages sur ce wiki, il faut utiliser la syntaxe habituelle des wikis ce qui inclut des formules LaTeX encadrées par <math>...</math>

Liste de diffusion : en cours de mise au point

Outil de gestion de versions : cvs ou svn. Accessible seulement à l'enseignant.

Sujets proposés

Evolution d'une galaxie

Système mécanique composé

...

Programme en C++ montrant les diverses notions vues en cours

En cours d'écriture

Vous n'êtes pas forcés de tout réécrire, vous pouvez télécharger les fichiers ici : [Disque.h] et [Disque.C]. Une fois que vous les avez téléchargés, suivez les instructions à la fin de ce paragraphe pour les utiliser dans ROOT. Avant, il faut comprendre ce qu'ils contiennent.

Les commentaires (après les //) devraient donner pas mal d'explications.

Voyons d'abord le fichier en-tête, déclarant les classes, fonctions, etc... des différents programmes. On met toujours ces déclarations dans un fichier séparé pour pouvoir les réutiliser. On nomme ce fichier "Disque.h" :

// La classe suivante definit un objet "disque" que l'on manipulera plus tard

// Declaration de la classe

class Disque
{
   private:
      // Definition des variables privees
      double mXC;        // coordonnees du centre
      double mYC;
      double mEpaisseur;  // Epaisseur si en 3D
      double mR;            // rayon du disque

   public:
      
      // Methodes "getter" pour acceder aux variables privees 
      double GetXC();
      double GetYC();
      double GetEpaisseur();
      double GetRayon();
      // Methodes "setter" pour initialiser les variables privees
      void SetXC(double x);
      void SetYC(double x);
      void SetEpaisseur(double x);
      void SetRayon(double x);
     
      // Methodes autres
      void Trace();                            // trace le disque
      void Deplace(double vX, double vY);  // Deplace le disque de (vX, vY)
      double Surface();                       // Renvoie la surface du disque
};

Voici le contenu du fichier "Disque.C". On commence par utiliser la directive "#include" pour inclure les fichiers d'en-tête.

#include "Disque.h"
#include "TCanvas.h"
#include "TEllipse.h"

Ensuite, on écrit le contenu des différentes méthodes de la classe en commençant par le constructeur

// ======================
//    Constructeur
// ======================
Disque::Disque(double xC, double yC, double epaisseur, double R)
{
   mXC = xC;
   mYC = yC;
   mEpaisseur = epaisseur;
   mR = R;
} 

Les fonctions qui permettent de récupérer les informations "privées" de la classe. Ces fonctions ne sont pas absolument indispensables pour faire tourner le programme, sauf si on a besoin d'accéder à des données privées. On est en plein dans les conséquences de la notion "d'encapsulation".

// ===========================
//    "Getter" et "Setter"
// ===========================

double Disque::GetXC()
{
   return mXC;
}

double Disque::GetYC()
{
   return mYC;
}

double Disque::GetEpaisseur()
{
   return mEpaisseur;
}

double Disque::GetRayon()
{
   return mR;
}

void Disque::SetXC(double x)
{
   mXC = x;
}

void Disque::SetYC(double x)
{
   mYC = x;
}

void Disque::SetEpaisseur(double x)
{
   mEpaisseur = x;
}

void Disque::SetRayon(double x)
{
   mR = x;
}

On passe maintenant aux méthodes spécifiques à la classe. Ici, on définit une méthode de tracé, une de déplacement et une de calcul simple.

// =========================================================
//     Methodes de visualisation, deplacement
// =========================================================

void Disque::Trace()
{
// trace le disque dans un canvas
   TEllipse* ell;
   
// **** ATTENTION, fuite de memoire !!!! ****
// on ne detruit jamais l'ellipse "ell" !
   ell = new TEllipse(mXC,mYC,mR);
   ell->Draw();
}

void Disque::Deplace(double vX, double vY)
{
// Deplace le disque de (vX, vY)
   mXC += vX;
   mYC += vY;
// Probleme : on n'est pas lie a l'objet ellipse cree dans la methode "Trace"
}

// ============================
//     Methodes de calcul
// ============================
double Disque::Surface()
{
// Renvoie la surface du disque
   double surf;
   double PI = 4*atan(1);
   
   surf = PI*mR*mR;
   return surf;
}

Utilisation de la classe dans ROOT

Avant de l'inclure dans un programme, il est bon de voir ce que l'on peut faire avec cette classe dans ROOT et comment on appelle chacune des méthodes. Comme ROOT contient un interpréteur C/C++, chaque ligne que nous écrirons sur la ligne de commande sera identique à ce que nous ferons plus tard dans un programme complet.

Première chose, il faut charger la classe dans ROOT (une fois celui-ci lancé dans le répertoire où se trouvent Disque.h et Disque.C, bien sûr). Pour celà, on utilise la commande ROOT ".L" :

root[0] .L Disque.C

Construisons ensuite un objet de la classe "Disque". On va utiliser la commande "new" pour créer l'objet et récupérer un pointeur sur cet objet :

root[1] d = new Disque(0.5, 0.5, 0.1, 0.1)

Ce faisant, on a appelé le constructeur de la classe disque. Ce constructeur, nous l'avons déclaré et défini pour qu'il prenne 4 paramètres en entrée et qu'il mette ces paramètres dans les variables privées de la classe. Si vous regardez le code, vous verrez que nous venons de définir un objet "Disque" de coordonnées du centre (0.5, 0.5), d'épaisseur 0.1 et de rayon 0.1.

Si l'on appelle maintenant la méthode de tracé, ce disque va se tracer tout seul : root[2] d->Trace()

Et voilà ! On voit apparaître un disque dans une fenêtre, appelée "Canvas" dans ROOT. Si on demande la surface du disque en appelant la méthode "Surface()", on obtient :

root [3] d->Surface()                   
(double)3.14159265358979339e-02

Là encore, on a appelé la méthode de l'objet d qui renvoie la valeur de la surface. Remarquez l'emploi de l'opérateur "->" (signe "-" suivi du signe supérieur ">" ). On appelle cet opérateur l'opérateur de déréférencement et il est employé pour indiquer "la méthode de l'objet pointé par ...". "d" est un pointeur, donc la méthode "Surface()" de l'objet pointé par "d" est appelée en utilisant "->".

Analyse d'une macro ROOT

On appelle "macro" ou "script" un programme ou morceau de programme que ROOT va interpréter. Analysons la macro suivante "graph.C" :

void graph() {
   //Draw a simple graph
   //Author: Rene Brun
  
   TCanvas *c1 = new TCanvas("c1","A Simple Graph Example",200,10,700,500);

   c1->SetFillColor(42);
   c1->SetGrid();

   const Int_t n = 20;
   double x[n], y[n];

   for (Int_t i=0;i<n;i++) {
     x[i] = i*0.1;
     y[i] = 10*sin(x[i]+0.2);
     printf(" i %i %f %f \n",i,x[i],y[i]);
   }

   gr = new TGraph(n,x,y);

   gr->SetLineColor(2);
   gr->SetLineWidth(4);
   gr->SetMarkerColor(4);
   gr->SetMarkerStyle(21);
   gr->SetTitle("a simple graph");
   gr->GetXaxis()->SetTitle("X title");
   gr->GetYaxis()->SetTitle("Y title");
   gr->Draw("ACP");
 
   // TCanvas::Update() draws the frame, after which one can change it
   c1->Update();
   c1->GetFrame()->SetFillColor(21);
   c1->GetFrame()->SetBorderSize(12);
   c1->Modified();
}

Vous pouvez l'exécuter dans ROOT, c'est une macro "standard" présente dans le répertoire ROOT. Faites simplement :

root[4] .x $ROOTSYS/tutorials/graphs/graph.C

Reprenons les lignes principales. Vous pouvez tester l'effet de certaines en recopiant le fichier chez vous et en enlevant la ou les lignes correspondantes avant de lancer la macro.

d'abord, la déclaration de la fonction. Elle ne prend pas d'argument :

void graph()
{

Ensuite, on construit un objet "TCanvas" qui est une fenêtre graphique dans laquelle on va dessiner :

   TCanvas *c1 = new TCanvas("c1","A Simple Graph Example",200,10,700,500);

On en profite pour appeler des méthodes de cet objet "TCanvas" qui remplissent le fond d'une certaine couleur

   c1->SetFillColor(42);

et qui activent une option traçant une grille pointillée

   c1->SetGrid();

viennent ensuite les déclarations de variables que l'on va utiliser. Ici, on a besoin de déclarer deux tableaux qui vont représenter des coordonnées de points à tracer. n est le nombre de points.

   const Int_t n = 20;
   double x[n], y[n];

ces points, il faut bien leur donner des valeurs, d'où la boucle ci-dessous :

   for (Int_t i=0;i<n;i++) {
     x[i] = i*0.1;
     y[i] = 10*sin(x[i]+0.2);
     printf(" i %i %f %f \n",i,x[i],y[i]);
   }

ensuite, on construit un objet "graphique" défini dans ROOT qui s'appelle "TGraph". On passe les pointeurs "x" et "y" sur les tableaux qu'on vient de remplir à ce constructeur, ainsi que le nombre de points :

   gr = new TGraph(n,x,y);

et on rend le graphe "joli". Les noms des méthodes parlent d'eux-même :

   gr->SetLineColor(2);
   gr->SetLineWidth(4);
   gr->SetMarkerColor(4);
   gr->SetMarkerStyle(21);
   gr->SetTitle("a simple graph");
   gr->GetXaxis()->SetTitle("X title");
   gr->GetYaxis()->SetTitle("Y title");
   gr->Draw("ACP");

enfin, l'auteur a décidé de changer encore un peu le fond, libre à lui...

   // TCanvas::Update() draws the frame, after which one can change it
   c1->Update();
   c1->GetFrame()->SetFillColor(21);
   c1->GetFrame()->SetBorderSize(12);
   c1->Modified();
}

Amusez-vous à modifier le code ci-dessus et à tester vos modifications dans ROOT. Par exemple, mettre plus de points ou changer certaines caractéristiques de couleur. C'est la meilleure manière d'apprendre ou de vérifier que vous avez bien compris. La classe TGraph est documentée sur le site de ROOT à la page : [http://root.cern.ch/root/html520/TGraph.html]

Il y a beaucoup de méthodes. Regardez certaines d'entre elles pour voir si vous comprenez leur utilité.

Il y a même un exemple encore plus simple de macro au début de cette page.

Il y a dans ROOT un très grand nombre de classes différentes, nous verrons celles dont nous aurons besoin au fur et à mesure.

Méthodes numériques de résolutions d'équations différentielles

I. Méthode d’Euler.

Cette méthode est la plus simple mais, en même temps, la moins précise.

Il peut néanmoins être utile de la détailler car certaines méthodes de résolution numériques ont pour base le calcul d’Euler, en y ajoutant certaines modifications (nous verrons cela plus loin). 

A. Principe.

Une fonction y quelconque, par exemple la position d’un objet en fonction du temps, admet le développement en série de Taylor suivant : y(t+∆t)=y(t)+∆t dy/dt+〖∆t〗^2/2! dy/dx+〖∆t〗^3/3! (d^2 y)/〖dx〗^2 +⋯ (I.A_1) en supposant y(0)=y0 connu. A l'ordre 1, on obtient : y(t+∆t)=y(t)+∆t dy/dt=y(t)+∆t*f(y,t) (I.A_2) Pour Δt constant, on a : 〖y_(i+1)=y〗_i+∆t*f(y_i,t_i) (I.A_3) La fonction y est connue à t=0, donc, la fonction y peut être déterminée à tout autre instant ultérieur.

B. Avantages et inconvénients.

Cette méthode est la plus simple à mettre en œuvre mais elle conduit très rapidement à des erreurs notamment parce qu’elle effectue l’intégration de y’(xi) à y’(xi+1) à partir de la valeur calculée en xi. Nous ne la retiendrons donc pas.

II. Méthodes de Runge-Kutta.

Les méthodes de Runge-Kutta sont basées sur la méthode d’Euler. Elles visent à approximer le résultat du calcul entre xi et xi+1 en estimant des valeurs intermédiaires sur cet intervalle.

1. Méthode d’ordre 2.

Sur l’intervalle [xi ; xi+1], on calcule la valeur que prend la fonction au milieu de cet intervalle et on calcul la dérivée en ce point. Supposons yi connu. On centre le calcul sur le milieu de l’intervalle. On résoudra donc :

y_(i+1)=y_i+∆t*f(y_(i+1/2)^a,t_(i+1/2)) (II.A_1) On rappelle la formule de Taylor : y(t+∆t)=y(t)+∆t dy/dt+〖∆t〗^2/2! dy/dx+〖∆t〗^3/3! (d^2 y)/〖dx〗^2 +⋯

On estime la valeur y ai+1/2 par : y_(i+1/2)^a=y_i+∆t/2 f(y_i,t_i) (II.A_2) D’autre part, on a évidemment : t_(i+1/2)=t_i+∆t/2 (II.A_3) Alors : y_(i+1)=y_i+∆t*f(y_(i+1/2)^a,t_(i+1/2)) (II.A_4)

A. Méthode d’ordre 4.

Le principe est le même qu’à l’ordre deux (‘ordre deux’ parce que deux évaluations successives de f sont effectuées : celle pour calculer yi+1/2 et celle permettant de trouver yi+1). Cette méthode consiste en fait à évaluer plusieurs fois la valeur de f(yi+1/2) entre yi et yi+1 : trois fois pour être précis. On part de la relation suivante :

y_(i+1)=y_i+∆t/6 [f(y_i,t_i )+2f(y_(i+1/2)^a,t_(i+1/2) )+2f(y_(i+1/2)^b,t_(i+1/2) )+f(y_(i+1)^c,t_(i+1) ) ] (II.B_1) Dans laquelle : y_(i+1/2)^a=y_i+∆t/2 f(y_i,t_i ) y_(i+1/2)^b=y_i+∆t/2 f(y_(i+1/2)^a,t_(i+1/2)) y_(i+1/2)^c=y_i+∆t/2 f(y_(i+1/2)^b,t_(i+1/2)) (II.B_2)

B. Avantages et inconvénients.

Les méthodes ci-dessus présentent plusieurs avantages : facilité de programmation, stabilité de la solution, modification simple du pas et la connaissance de y(0) suffit pour intégrer l'équation différentielle. Plus l’ordre de la méthode est élevé, meilleure est la précision. La conséquence immédiate est l’augmentation des calculs à effectuer. Les inconvénients de cette méthode se résument donc au temps de calcul lent et à la difficulté de l'estimation de l'erreur locale. A noter néanmoins que la méthode de Runge-Kutta d’ordre 4 semble être assez largement utilisée en simulation de phénomènes physiques (gage d’efficacité ??).

III. Méthodes d’Adams.

Les méthodes suivantes ont également pour base la méthode d’Euler, en particulier elles s’appuient sur le développement de Taylor y(t+∆t)=y(t)+∆t dy/dt+〖∆t〗^2/2! dy/dx+〖∆t〗^3/3! (d^2 y)/〖dx〗^2 +⋯ (III_1) en supposant y(0)=y0 connu. Il existe deux méthodes, ou formules d’Adams : les formules dites « ouvertes » ; les formules dites « fermées ».

A. Formules d’Adams ouvertes.

En posant f’(y,t)=dy/dt on peut réécrire immédiatement le développement de Taylor y_(i+1)=y_i+∆t*f_i+〖∆t〗^2/2! f_i^((1))+〖∆t〗^3/3! f_i^((2))+〖∆t〗^4/4! f_i^((3))+⋯ (III.A_1)

1. Formules d’Adams ouvertes d’ordre 2.

La résolution s’effectue en tenant compte du DL de Taylor jusqu’à l’ordre 2 (d’où le nom de la méthode). y_(i+1)=y_i+∆t*f_i+〖∆t〗^2/2! f_i^((1))+〖Ο(∆t)〗^3 (III.A.1_1) Or, on obtient la dérivée première d’une fonction par f^' (t)=lim┬(∆t→0)⁡〖(f(t+∆t)-f(t))/Δt〗 (III.A.1_2) On peut donc écrire f_i^((1) )=(f_i-f_(i-1))/Δt (III.A.1_3) où fi = f(yi, ti) et fi-1 = f(yi-1, ti-1). On peut alors calculer yi+1 par y_(i+1)=y_i+∆t*f_i+〖∆t〗^2/2! ((f_i-f_(i-1))/Δt) ⟺y_(i+1)=y_i+2∆t f_i/2+〖∆t〗^2/2 ((f_i-f_(i-1))/Δt) ⟺y_(i+1)=y_i+∆t/2 (〖3f〗_i-f_(i-1) ) (III.A.1_4) Remarquons au passage que le terme y1 ne peut pas être calculé par cette méthode car il suppose la connaissance du terme y-1, ce qui n’est pas le cas. Ce terme doit donc être évalué par une autre méthode.

2. Formules d’Adams ouvertes d’ordre n.

Généralisation de la méthode précédente, elle consiste à faire le calcul en prenant l’ordre du DL de Taylor à l’ordre n et à exprimer yi+1 sous forme d’une combinaison linéaire de fi-k, k variant de 0 à n. Pour cela, on utilise la même méthode que ci-dessus : on réinjecte dans le DL l’expression des n-1 dérivées à calculer, évaluées comme la relation (III.A.1_3) le fait. La relation utilisée est la suivante : y_(i+1)=y_i+∆t∑_(k=0)^n▒〖β_nk f_(i-k) 〗 (III.A.2_1) Les βnk sont des coefficients, calculables par une formule qui nous est inconnue, mais déterminables par le calcul explicite. Voici un tableau récapitulatif de ces coefficients pour une formule d’ordre 6.

(Source : http://www.sciences.univ-nantes.fr/physique/perso/aloui/m_numeri/61eqdiff/61eqdiff.htm#1)

Avantages et inconvénients. A notre connaissance, il n’y a pas d’avantage particulier déclaré hormis la précision accrue devant la méthode d’Euler. Par contre, comme vu pour la méthode d’ordre 1, il faut évaluer le terme y1 par une méthode autre. De même, pour une méthode d’ordre n, il faut évaluer les yk, k variant de 1 à n-1, par une méthode alternative, à déterminer par ailleurs (Runge-Kutta, ou pourquoi pas une utilisation successive des formules d’Adams pour les évaluer et partir construire ainsi progressivement la formule au rang n que l’on souhaite…). La plus utilisée semble être la méthode d’ordre 4. Formules d’Adams fermées. Cette méthode est similaire à la précédente, sauf que l’on utilise le DL de Taylor à ‘l’envers’.