- Le chapitre 10 de « Methodes Numeriques, Algorithmes, analyse et applications ». Le livre est disponible sur mon compte à l’enseirb ~mfaverge/public_html/IS104/Methodes-Numeriques-Algorithmes-analyse-et-applications.pdf
- Les slides de présentation de ce chapitre sont disponibles ici. Ils sont également disponibles sur mon compte à l’enseir au même endroit que le livre.
- Quelques slides additionnels pour présenter les méthodes du point milieu, de Heun, et de Runge-Kutta 4.
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 :
-
- le modèle de biologie des populations dit de Lotka-Volterra,
- le pendule à
maillons,
- le traitement d’image avec la méthode de Perona et Malik,
- les réactions chimiques oscillantes, comme celles de Belousov-Zhabotinsky,
- le problème à trois corps.
-
- La deuxième partie est ensuite constituée d’applications de ces méthodes à un certain nombre de systèmes dynamiques :
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
.png)
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.
-
- 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.
- 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.
-
- 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.
- la méthode d’Euler
-
- Programmez les méthodes de résolution suivantes :
-
- 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 nombreN
de pas de taille constanteh
, - puis un second algorithme
meth_epsilon(y0,t0,tf,eps,f,meth)
calculant une solution approchée avec un paramètre d’erreurepsilon
.
- pour chaque méthode de résolution, une fonction
Comparer l’utilité de l’implémentation générique par rapport au projet précédent.
-
- Écrire 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
)
- Écrire 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
-
- 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 :
- En dimension 1 :
Quels sont les types d’erreurs que vous rencontrez ?
-
- 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 :
Système proie-prédateur de Lotka-Volterra
.png)
}{dt}=.png)
-dN(t)=gammaN(t)quadbd>0.png)
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 :
}{dt}=displaystylegammaN(t)left(1-frac{N(t)}{kappa}right)quadgammakappa>0.png)
-
- 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.
-
- 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 : est une population de proies et
est une population de prédateurs. Le modèle de Lotka-Volterra est le suivant :
}{dt}=N(t)(a-bP(t))frac{dP(t)}{dt}=P(t)(cN(t)-d)end{cases}quadabcd>0.png)
-
- Expliquer en quoi ces deux équations peuvent-elles modéliser un écosystème proie/prédateur.
-
- 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.
-
- Tracer sur une courbe les variations du couple
au cours du temps. Quels types de solutions fait-on apparaître ? Quelles sont les solutions constantes ?
- Tracer sur une courbe les variations du couple
-
- 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.
-
- Écrire un algorithme qui permette de tracer les solutions autour d’un point de départ donné. En fonction d’une condition initiale
, il s’agit de tracer la courbe solution en partant de
, puis les courbes solutions en partant de conditions initiales qui sont proches de
(il est possible de les choisir dans un petit disque/rectangle centré sur
).
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.)
- Écrire un algorithme qui permette de tracer les solutions autour d’un point de départ donné. En fonction d’une condition initiale
-
- 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 ?
-
- 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 ?
- 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 ?
Pendule à N maillons



![{\begin{pspicture}(-2.2,-1.2)(2.2,2.8) \pspolygon[fillstyle=solid,linestyle=none,fillcolor=gray] (-2,2)(2,2)(2,2.5)(-2,2.5) \psline(-2,2)(2,2) \psline(0,2)(1,0.5)(-0.5,-0.5) \psline[linestyle=dashed](0,2)(0,0.5) \psarc{->}(0,2){0.7}{270}{303} \rput(0.3,1){\small$\theta_1$} \pscircle*(1,0.5){0.1} \psline[linestyle=dashed](1,0.5)(1,-1) \pscircle*(-0.5,-0.5){0.1} \psarcn{->}(1,0.5){0.6}{270}{214} \rput(0.6,-0.3){\small$\theta_2$} \pcline[offset=8pt,linecolor=gray]{<->}(-0.5,-0.5)(1,0.5) \lput*{:U}{\small$l$} \end{pspicture}}](wp-content/cours/IS104/images/p6_double_pendulum.png)
Ici, les angles


-
- 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
. En traçant cette fréquence en fonction de
, montrer que cette fréquence pour de petites oscillations est égale à
.
- 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
-
- 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.
-
- 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.
- 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)
-
- 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..
Filtrage de Perona et Malik



.png)
-Deltau(xyt)=0text{pour}t>0u(0xy)=u_0(xy)end{array}right..png)
… pour



Pour résoudre effectivement ce problème, on considère que l’image d’entrée est une version discrétisée de la fonction au temps
. Le laplacien
est calculé, à temps
fixé, via
des opérateurs discrets de divergence et de gradient : . À temps
fixé, le gradient est un champ de vecteur de taille
, et chaque vecteur a deux dimensions.
![]() |
|
![]() |
La divergence est alors l’opérateur conjugué. Il prend un champ de vecteur à deux dimensions en entrée (le résultat du gradient par exemple), et est défini comme suit :
_{kl}=.png)
![\left\lbrace \begin{array}{c} p^1_{k,l}- p^1_{k-1,l} \text{ si } 1<k<N \text{,}\\ \rule[-2ex]{0pt}{5.5ex} p^1_{k,l} \text{ si } k=1 \text{,} \\ \rule[-2ex]{0pt}{5.5ex} -p^1_{k-1,l} \text{ si } k=N \text{,} \end{array} \right.](wp-content/cours/IS104/images/leftlbracebegin{array}{c}p^1_{kl}-p^1_{k-1l}text{si}1<k<Ntext{}rule[-2ex]{0pt}{5.5ex}p^1_{kl}text{si}k=1text{}rule[-2ex]{0pt}{5.5ex}-p^1_{k-1l}text{si}k=Ntext{}end{array}right..png)

![\left\lbrace \begin{array}{c} p^2_{k,l}- p^2_{k,l-1} \text{ si } 1<l<M \text{,}\\ \rule[-2ex]{0pt}{5.5ex} p^2_{k,l} \text{ si } l=1 \text{,} \\ \rule[-2ex]{0pt}{5.5ex} -p^2_{k,l-1} \text{ si } l=M \text{,} \end{array} \right.](wp-content/cours/IS104/images/leftlbracebegin{array}{c}p^2_{kl}-p^2_{kl-1}text{si}1<l<Mtext{}rule[-2ex]{0pt}{5.5ex}p^2_{kl}text{si}l=1text{}rule[-2ex]{0pt}{5.5ex}-p^2_{kl-1}text{si}l=Mtext{}end{array}right..png)
avec
_{kl}.png)


.png)
![]() |
![]() |
![]() |
Image originale | Équation de la chaleur | Filtre de Perona Malik |
La méthode de Perona et Malik évite de flouter l’image en préservant les contours. Celle-ci repose sur l’EDP suivante :
=mathrm{div}(f(xy).nablau(xyt))text{pour}t>0u(0xy)=u_0(xy)end{array}right..png)


![[0,1]](wp-content/cours/IS104/images/[01].png)
La fonction se calcule à partir de la fonction initiale et permet de limiter l’action du filtrage (floutage) au niveau des zones de contours. En pratique, c’est un tableau de coefficients de la taille de l’image. L’équation de la chaleur correspond à celle de Perona-Malik où
.
Afin de ne pas flouter uniformément l’image, on réduit l’action de l’opérateur laplacien au niveau des zones de contours en multipliant le gradient par une fonction ayant une valeur proche de 1 aux endroits où il n’y a pas de contour (norme du gradient faible) et ayant une valeur proche de 0 aux endroits où il y a des contours (norme du gradient élevée). Une telle fonction peut être choisie de la forme avec
une convolution de l’image initiale avec une gaussienne. Cette convolution vise à élargir la zone de coefficients faibles afin de mieux préserver les contours.
-
- Implémenter les gradients et divergence de l’image.
-
- Implémenter le schéma d’Euler sur l’équation de la chaleur (
, prendre un pas de temps
)
- Implémenter le schéma d’Euler sur l’équation de la chaleur (
-
- Implémenter la convolution de l’image avec une gaussienne en prenant soin de considérer que les bords de l’image soient symétriques (utiliser
scipy.signal.convolve2d
). Construire la fonctionet la visualiser.
- Implémenter la convolution de l’image avec une gaussienne en prenant soin de considérer que les bords de l’image soient symétriques (utiliser
-
- Implémenter le schéma d’Euler correspondant à l’équation de Perona et Malik avec un pas de temps de 0.2 et 200 itérations.
Réactions chimiques oscillantes
Le système auquel nous nous intéressons ici est un sytème oscillant simplifié, à trois dimensions, modélisé par les équations suivantes :
}{dt}=a+bZ-X-XY^2frac{dY(t)}{dt}=c(X+XY^2-Y)frac{dZ(t)}{dt}=d(Y-Z)end{cases}quadabcd>0.png)
Pour les conditions de la réaction, les valeurs intéressantes se situent pour



![b\in [0.1,0.2]](wp-content/cours/IS104/images/bin[0.10.2].png)
![[1; 0; 0]](wp-content/cours/IS104/images/[100].png)

-
- Résoudre les équations précédentes pour différentes valeurs du paramètre
.
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.
- Résoudre les équations précédentes pour différentes valeurs du paramètre
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 : étant donné une solution , la période de cette solution est égale au nombre de maxima locaux de
restreinte à sa phase périodique.
-
- Comment ne considérer que la phase périodique d’une solution ?
- Écrire une fonction qui, étant donné une valeur du paramètre
, 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 ?
- 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
.
- 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 ?
- En quoi la modélisation numérique d’un système dynamique chaotique est elle difficile ?
Modélisation du problème à trois corps
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 dûe à 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
.
-
- 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(ce qui réduit proportionnellement la force subie par
A
) et que les masses et distances sont des valeurs simples.
- Ecrire les équations du mouvement associées au corps
-
- 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 ?
- 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.
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 :

Exemple :

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 
-
- 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 deA
, 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.
- 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
Comme le mouvement de B
est connu, il est possible d’effectuer un changement de base apres la fin du calcul (une rotation), de maniere 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
.
-
- Transformer les trajectoires obtenues en appliquant à chaque point la rotation faisant que la planète
B
reste fixe.
- Transformer les trajectoires obtenues en appliquant à chaque point la rotation faisant que la planète
-
- 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 ?
- Que se passe-t’il lorsque l’on choisit deux points de départ pour le corps
-
- Rechercher dans ce système les points particuliers évoqués au cours d’un précédet 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 equilateral dont la base est formée par le segment Soleil-Jupiter. Que constatez-vous ?
- Rechercher dans ce système les points particuliers évoqués au cours d’un précédet 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.
Noter qu’il existe des solutions périodiques au problème général à n corps, dont certaines qu’il est possible de visualiser là.