Résolution approchée d'équations différentielles / Modélisation de systèmes dynamiques
Sujets de David Renault
Pour cet ultime projet, nous allons nous intéresser aux méthodes de résolution d'équations différentielles ordinaires (EDO). L'intérêt de ces équations réside dans le fait qu'elles permettent de modéliser des systèmes complexes relativement facilement. Ce qui ne signifie pas forcément que leur résolution ne pose pas de problèmes. Ce projet comporte deux parties :
Dans un premier temps, il s'agit de programmer les méthodes de résolution numérique à un pas, adaptées à des équations différentielles ordinaires en dimension finie, et de les tester dans un cadre simple.
La deuxième partie est ensuite constituée d'applications de ces méthodes à un certain nombre de systèmes dynamiques :
Dans un premier temps, il s'agit de programmer les méthodes de résolution numérique à un pas, adaptées à des équations différentielles ordinaires en dimension finie, et de les tester dans un cadre simple.
La deuxième partie est ensuite constituée d'applications de ces méthodes à un certain nombre de systèmes dynamiques :
- le modèle de biologie des populations dit de Lotka-Volterra,
- le pendule à N maillons,
- les réactions chimiques oscillantes, comme celles de Belousov-Zhabotinsky,
- le problème à trois corps.
Ce projet laisse une part importante à l'expérimentation et à la recherche de paramètres mettant en évidence des comportements particuliers des équations différentielles rencontrées. Il est donc plus que conseillé de faire preuve d'initiative et de tester plusieurs exemples par soi-même.
Il est demandé de réaliser au moins deux applications pour mettre en valeur les problèmes numériques liés à la résolution d'équations différentielles.
Il est demandé de réaliser au moins deux applications pour mettre en valeur les problèmes numériques liés à la résolution d'équations différentielles.
Méthodes numériques de résolution d'équations différentielles
Dans cette partie, on programme les différentes méthodes numériques de résolution des équations différentielles que nous testerons par la suite. Ces algorithmes sont tous des méthodes à un pas, de la forme suivante :
Les algorithmes reprenant tous le même principe, il est possible de répondre à une grande partie des questions sans programmer toutes les méthodes demandées, pour y revenir plus tard lorsqu'elles auront été vues en cours.
1. Définissez la façon dont vous allez représenter un problème de Cauchy dans votre code. Il s'agit de fixer la façon dont vous allez représenter la condition initiale et la fonction différentielle représentant l'équation, afin d'unifier le code écrit dans chacune des parties suivantes.
1. Définissez la façon dont vous allez représenter un problème de Cauchy dans votre code. Il s'agit de fixer la façon dont vous allez représenter la condition initiale et la fonction différentielle représentant l'équation, afin d'unifier le code écrit dans chacune des parties suivantes.
Votre représentation doit pouvoir prendre en compte des équations différentielles en dimension arbitraire.
2. Programmez les méthodes de résolution suivantes :
2. Programmez les méthodes de résolution suivantes :
- la méthode d'Euler :
- la méthode du point milieu,
- la méthode de Heun,
- et la méthode de Runge-Kutta d'ordre 4.
Remarque :
Pour chacune de ces méthodes on programmera, comme pour les méthodes d'intégration, de la manière la plus générique possible :
Pour chacune de ces méthodes on programmera, comme pour les méthodes d'intégration, de la manière la plus générique possible :
- pour chaque méthode de résolution, une fonction step_<name>(y,t,h,f) qui calcule un seul pas de la méthode,
- un premier algorithme meth_n_step(y0,t0,N,h,f,meth) calculant un nombre N de pas de taille constante h,
- puis un second algorithme meth_epsilon(y0,t0,tf,eps,f,meth) calculant une solution approchée avec un paramètre d'erreur epsilon.
3. Ecrire une fonction permettant de dessiner le champ des tangentes de l'équation différentielle (on ne s'intéressera qu'au cas de la dimension 2, et on utilisera la fonction quiver).
4. Testez vos implémentation sur des équations différentielles connues, et comparez-les aux résultats exacts. En particulier, on tracera les courbes correspondant aux deux équations différentielles suivantes :
En dimension 1 :
4. Testez vos implémentation sur des équations différentielles connues, et comparez-les aux résultats exacts. En particulier, on tracera les courbes correspondant aux deux équations différentielles suivantes :
En dimension 1 :
En dimension 2 :
Quels sont les types d'erreurs que vous rencontrez ?
Système proie-prédacteur de Lotka-Volterra
Dans cette partie, on s'intéresse à la modélisation de l'évolution de la population d'une ou plusieurs espèces dans un milieu naturel. Les équations différentielles décrites ici sont des modèles de population. Considérons la fonction N(t) décrivant la variation de la population d'un ensemble d'individus au cours du temps. Les premiers modèles de population (dits malthusiens, car dus à Malthus à la fin du XVIIIème siècle) décrivaient ainsi les variations de cette fonction :
Néanmoins, cette approche n'est pas vraiment réaliste (même si la croissance de la population humaine ne contredit pas effectivement le modèle), et Verhulst a proposé un autre modèle de population auto-limitant :
1. Expliquer en quoi ces deux équations différentielles peuvent modéliser les variations d'une population, et en particulier à quoi correspondent les constantes apparaissant à l'intérieur.
2. Résoudre ces équations différentielles en utilisant une méthode de votre choix, et en choisissant différents points de départ crédibles.
Ces modèles peuvent être modifiés à loisir, mais lorsque l'on veut décrire un écosystème de manière un peu plus riche, il est possible de modéliser les interactions de deux populations proies/prédateurs : N(t) est une population de proies et P(t) est une population de prédateurs. Le modèle de Lotka-Volterra est le suivant :
2. Résoudre ces équations différentielles en utilisant une méthode de votre choix, et en choisissant différents points de départ crédibles.
Ces modèles peuvent être modifiés à loisir, mais lorsque l'on veut décrire un écosystème de manière un peu plus riche, il est possible de modéliser les interactions de deux populations proies/prédateurs : N(t) est une population de proies et P(t) est une population de prédateurs. Le modèle de Lotka-Volterra est le suivant :
3. Expliquer en quoi ces deux équations peuvent-elles modéliser un écosystème proie/prédateur.
4. Résoudre numériquement ce système d'équations différentielles en utilisant une méthode de votre choix, et en choisissant différents points de départ crédibles. Tracer les variations des deux populations au cours du temps.
5. Tracer sur une courbe les variations du couple (N(t), P(t)) au cours du temps. Quels types de solutions fait-on apparaître ? Quelles sont les solutions constantes ?
6. Mettre en place un algorithme permettant de calculer la période de ces solutions (de manière approchée).
Lors de l'étude d'un système d'équations différentielles, on s'intéresse au comportement global des solutions, mais aussi à leur comportement local. En particulier, il est intéressant de regarder les variations entre des solutions partant de conditions initiales très proches.
7. Ecrire un algorithme qui permette de tracer les solutions autour d'un point de départ donné. En fonction d'une condition initiale y0, il s'agit de tracer la courbe solution en partant de y0, puis les courbes solutions en partant de conditions initiales qui sont proches de y0 (il est possible de les choisir dans un petit disque/rectangle centré sur y0).
Le résultat pourra ressembler à :
4. Résoudre numériquement ce système d'équations différentielles en utilisant une méthode de votre choix, et en choisissant différents points de départ crédibles. Tracer les variations des deux populations au cours du temps.
5. Tracer sur une courbe les variations du couple (N(t), P(t)) au cours du temps. Quels types de solutions fait-on apparaître ? Quelles sont les solutions constantes ?
6. Mettre en place un algorithme permettant de calculer la période de ces solutions (de manière approchée).
Lors de l'étude d'un système d'équations différentielles, on s'intéresse au comportement global des solutions, mais aussi à leur comportement local. En particulier, il est intéressant de regarder les variations entre des solutions partant de conditions initiales très proches.
7. Ecrire un algorithme qui permette de tracer les solutions autour d'un point de départ donné. En fonction d'une condition initiale y0, il s'agit de tracer la courbe solution en partant de y0, puis les courbes solutions en partant de conditions initiales qui sont proches de y0 (il est possible de les choisir dans un petit disque/rectangle centré sur y0).
Le résultat pourra ressembler à :
(On rappelle qu'il s'agit bien d'expliquer le comportement local de l'équation différentielle, et donc que la fenêtre considérée doit être petite. De plus, il est important de ne considérer que les EDO en dimension 2, afin de montrer des courbes solutions qui ne se croisent pas.)
8. Un point singulier de l'équation différentielle est un point où toutes les dérivées s'annulent. Quels sont les points singuliers du système de Lotka-Volterra ?
9. Quelle est la forme générale des diagrammes obtenus à la question 7 (en partant d'une condition initiale différente d'un point singulier) ? Quelle est la forme de ces diagrammes lorsque l'on part d'un des points singuliers de l'équation ?
(Difficile) Si l'on se replace dans le contexte d'une équation différentielle générale, à quels types de diagrammes peut-on s'attendre ?
8. Un point singulier de l'équation différentielle est un point où toutes les dérivées s'annulent. Quels sont les points singuliers du système de Lotka-Volterra ?
9. Quelle est la forme générale des diagrammes obtenus à la question 7 (en partant d'une condition initiale différente d'un point singulier) ? Quelle est la forme de ces diagrammes lorsque l'on part d'un des points singuliers de l'équation ?
(Difficile) Si l'on se replace dans le contexte d'une équation différentielle générale, à quels types de diagrammes peut-on s'attendre ?
Pendule à N maillons
Un système dynamique très simple peut être construit en considérant un pendule constitué d'un ensemble de N solides de masse m reliés entre eux par des tiges sans masse de même longueur l, les tiges étant supposées en rotation les unes par rapport aux autres (le système est supposé confiné à un plan). Par exemple, pour un pendule à deux maillons, le système est modélisé par le dessin suivant :
Ici, les angles theta_1 et theta_2 permettent de spécifier entièrement le système. On se propose alors d'examiner les évolutions de ce système par rapport au temps, à partir d'une position à vitesse nulle.
1. Modéliser le pendule à un unique maillon. Implémenter une fonction permettant de déterminer la fréquence d'oscillation du système, en fonction de l'angle initial theta. En traçant cette fréquence en fonction de theta, montrer que cette fréquence pour de petites oscillations est égale à :
1. Modéliser le pendule à un unique maillon. Implémenter une fonction permettant de déterminer la fréquence d'oscillation du système, en fonction de l'angle initial theta. En traçant cette fréquence en fonction de theta, montrer que cette fréquence pour de petites oscillations est égale à :
2. Modéliser le pendule à deux maillons. On pourra s'inspirer de la page suivante pour les équations.
Il se trouve que le pendule à deux maillons est un système dynamique chaotique, ce qui peut se voir à travers la sensibilité des trajectoires aux conditions initiales.
3. Dessiner la trajectoire de l'extrémité du pendule à deux maillons en fonction du temps (il est possible d'obtenir des séquences du style de cette vidéo) En considérant plusieurs trajectoires avec des conditions initiales différentes mais très proches, mettre en exergue la dépendance aux conditions initiales de ce système.
4. Une autre manière de mettre en valeur les propriétés du système consiste à déterminer le temps de premier retournement. Ecrire une fonction permettant de calculer ce temps, afin de tracer une carte du type de celle dessinées sur cette page...
Il se trouve que le pendule à deux maillons est un système dynamique chaotique, ce qui peut se voir à travers la sensibilité des trajectoires aux conditions initiales.
3. Dessiner la trajectoire de l'extrémité du pendule à deux maillons en fonction du temps (il est possible d'obtenir des séquences du style de cette vidéo) En considérant plusieurs trajectoires avec des conditions initiales différentes mais très proches, mettre en exergue la dépendance aux conditions initiales de ce système.
4. Une autre manière de mettre en valeur les propriétés du système consiste à déterminer le temps de premier retournement. Ecrire une fonction permettant de calculer ce temps, afin de tracer une carte du type de celle dessinées sur cette page...
Réactions chimiques oscillantes
Usuellement, les réactions chimiques ont un comportement simple d'un point de vue thermodynamique : la réaction arrive à un équilibre stable à partir duquel on peut considérer qu'il n'y a plus d'interaction entre les réactifs. Néanmoins, il existe des réactions chimiques dites oscillantes, qui ont un comportement périodique. La plus connue de ces réactions est la réacition de Belousov-Zhabotinsky, ou encore celle de Briggs-Schauer.
Le système auquel nous nous intéressons ici est un système oscillant simplifié, à trois dimensions, modélisé par les équations suivantes :
Le système auquel nous nous intéressons ici est un système oscillant simplifié, à trois dimensions, modélisé par les équations suivantes :
Pour les conditions de la réaction, les valeurs intéressantes se situent pour :
Les valeurs initiales sont fixées à [1; 0; 0]. Dans cette partie, on s'intéresse au comportement de la réaction en fonction du paramètre b.
1. Résoudre les équations précédentes pour différentes valeurs du paramètre b.
1. Résoudre les équations précédentes pour différentes valeurs du paramètre b.
Remarque : Les solutions obtenues doivent être périodiques, ou ultimement périodiques. Aucune solution ne diverge vers l'infini. Veillez à vérifier que vos calculs ne divergent pas, en prenant des pas de calculs suffisamment faibles.
Toutes les solutions ont un comportement ressemblant : premièrement, une phase de démarrage, puis une phase stable dans laquelle la fonction a un comportement périodique, ou quasi-périodique. Pour mesurer la complexité de ce comportement, nous allons définir une notion de période d'une solution :
Etant donné une solution [X(t); Y(t); Z(t)], la période de cette solution est égale au nombre de maxima locaux de Y(t) restreinte à sa phase périodique.
2. Comment ne considérer que la phase périodique d'une solution ?
3. Ecrire une fonction qui, étant donné une valeur du paramètre b, calcule la période de la solution associée.
Pouvez-vous exhiber des solutions de période 2 ? de période 4 ? de période 5 ? de période 3 ? Pouvez-vous exhiber des solutions de période infinie ?
4. Tracer le diagramme de bifurcation de ce système dynamique, à savoir une courbe traçant les pics de la solution en fonction du paramètre b.
Etant donné une solution [X(t); Y(t); Z(t)], la période de cette solution est égale au nombre de maxima locaux de Y(t) restreinte à sa phase périodique.
2. Comment ne considérer que la phase périodique d'une solution ?
3. Ecrire une fonction qui, étant donné une valeur du paramètre b, calcule la période de la solution associée.
Pouvez-vous exhiber des solutions de période 2 ? de période 4 ? de période 5 ? de période 3 ? Pouvez-vous exhiber des solutions de période infinie ?
4. Tracer le diagramme de bifurcation de ce système dynamique, à savoir une courbe traçant les pics de la solution en fonction du paramètre b.
5. Dans leur article Period 3 implies Chaos, Li et Yorke montrent qu'un système dynamique exhibant des solutions de période 3 devait nécessairement contenir des solutions chaotiques. Est-ce le cas du présent système ? Et sinon, est-ce tout de même un système chaotique ?
6. En quoi la modélisation numérique d'un système dynamique chaotique est elle difficile ?
6. En quoi la modélisation numérique d'un système dynamique chaotique est elle difficile ?
Modélisation du problème à trois corps
Dans cette partie, on s'intéresse à la modélisation d'un problème apparemment simple de dynamique du solide : le problème à trois corps. Dans un premier temps, on s'attache à modéliser les interactions gravitationnelles classiques entre deux astres, puis on introduit un troisième corps massif pour mettre en évidence les modifications du système.
Considérons un système physique composé de deux corps massifs A, B ponctuels placés dans un plan. Chacun de ces deux corps est soumis uniquement à la force d'interaction gravitationnelle due à la présence de l'autre corps. Les caractéristiques de chaque corps sont donc :
Considérons un système physique composé de deux corps massifs A, B ponctuels placés dans un plan. Chacun de ces deux corps est soumis uniquement à la force d'interaction gravitationnelle due à la présence de l'autre corps. Les caractéristiques de chaque corps sont donc :
Ces corps sont supposés de masse non nulle, placés initialement à des positions distinctes. Pour rappel, la force d'interaction gravitationnelle entre deux corps (par exemple la force que A subit sous l'attraction de B) est donnée par :
Par la suite, nous supposerons que la masse du corps A est largement supérieure à la masse du corps B. Cela simplifie le problème en impliquant que le corps A peut être supposé immobile, et insensible à l'interaction gravitationnelle de B.
1. Ecrire les équations du mouvement associées au corps B en utilisant l'équation fondamentale de la dynamique. De combien d'inconnues avez-vous besoin pour modéliser entièrement votre système ?
1. Ecrire les équations du mouvement associées au corps B en utilisant l'équation fondamentale de la dynamique. De combien d'inconnues avez-vous besoin pour modéliser entièrement votre système ?
Il est plus que fortement conseillé de rendre vos équations adimensionnelles, afin de simplifier les calculs. Cela consiste à simplifier les équations en éliminant les constantes inutiles. Ici en particulier, une bonne idée consiste à considérer que la constante G = 1 (ce qui réduit proportionnellement la force subie par A) et que les masses et distances sont des valeurs simples.
2. Résoudre ce système d'équations différentielles en utilisant une méthode de votre choix, et en choisissant différents points de départ crédibles. Tracer les évolutions de vos planétoïdes au cours du temps.
Quelles sont les courbes que l'on est censé obtenir ? Quels types de comportement erronés arrivez-vous à obtenir ?
La partie suivante s'intéresse à une modification du modèle précédent, appelée le problème à trois corps restreint. Dans ce nouveau système, on suppose que le système contient un troisième corps C, lui aussi de masse négligeable devant les masses de A et de B.
Quelles sont les courbes que l'on est censé obtenir ? Quels types de comportement erronés arrivez-vous à obtenir ?
La partie suivante s'intéresse à une modification du modèle précédent, appelée le problème à trois corps restreint. Dans ce nouveau système, on suppose que le système contient un troisième corps C, lui aussi de masse négligeable devant les masses de A et de B.
On pourra par exemple s'inspirer des valeurs suivantes :
Le centre de gravité du système se situe donc à peu près au centre de A, ce qui permet de se placer dans un référentiel placé au centre du soleil, supposé fixe. Le corps B est supposé en rotation circulaire autour du soleil, avec une période de 2Pi jours terrestres.
3. Modifier votre système d'équations pour modéliser ce nouveau système en 4 dimensions d'espace et une dimension de temps (Indice : on pourra supposer le mouvement de B connu et circulaire autour de A, et ne modéliser que le mouvement de l'astéroïde).
Résoudre ce système d'équations différentielles en utilisant une méthode de votre choix.
Comme le mouvement de B est connu, il est possible d'effectuer un changement de base après la fin du calcul (une rotation), de manière a ce que les corps A et B restent fixes. Il est alors possible de ne s'intéresser qu'aux mouvements du corps C.
4. Transformer les trajectoires obtenues en appliquant à chaque point la rotation faisant que la planète B reste fixe.
5. Que se passe-t-il lorsque l'on choisit deux points de départ pour le corps C qui sont initialement très proches l'un de l'autre ? En particulier, que deviennent les trajectoires au bout d'un certain temps ?
6. Rechercher dans ce système les points particuliers évoqués au cours d'un précédent projet, appelés points de Lagrange, qui correspondent à des points d'équilibre du système physique et des points singuliers du système d'équations différentielles.
Placez le corps C initialement proche du dernier sommet du triangle équilatéral dont la base est formée par le segment Soleil-Jupiter. Que constatez-vous ?
3. Modifier votre système d'équations pour modéliser ce nouveau système en 4 dimensions d'espace et une dimension de temps (Indice : on pourra supposer le mouvement de B connu et circulaire autour de A, et ne modéliser que le mouvement de l'astéroïde).
Résoudre ce système d'équations différentielles en utilisant une méthode de votre choix.
Comme le mouvement de B est connu, il est possible d'effectuer un changement de base après la fin du calcul (une rotation), de manière a ce que les corps A et B restent fixes. Il est alors possible de ne s'intéresser qu'aux mouvements du corps C.
4. Transformer les trajectoires obtenues en appliquant à chaque point la rotation faisant que la planète B reste fixe.
5. Que se passe-t-il lorsque l'on choisit deux points de départ pour le corps C qui sont initialement très proches l'un de l'autre ? En particulier, que deviennent les trajectoires au bout d'un certain temps ?
6. Rechercher dans ce système les points particuliers évoqués au cours d'un précédent projet, appelés points de Lagrange, qui correspondent à des points d'équilibre du système physique et des points singuliers du système d'équations différentielles.
Placez le corps C initialement proche du dernier sommet du triangle équilatéral dont la base est formée par le segment Soleil-Jupiter. Que constatez-vous ?