Download Résolution numérique des équations différentielles ordinaires and more Exercises Design in PDF only on Docsity!
Résolution numérique des équations différentielles ordinaires
April 29, 2024
Le but de ce TP est de résoudre, numériquement, à l’aide du logiciel Python, une équation ou systèmes d’équations différentielles ordinaires du premier ordre. On commence par importer les différents modules dont nous auront besoin:
[34]: import numpy as np import matplotlib.pyplot as plt
1 Equations du premier ordre
Considérons l’équation différentielle suivante:
{ x′(t) = 5 x^2 (t) − x(t)(1 + t^3 ), t > 0 , x(0) = 1.
C’est une équation différentielle d’ordre 1, mais elle n’est pas linéaire. Nous ne savons pas la résoudre de manière exacte. Nous allons toutefois pouvoir la résoudre numériquement… Remarquons déjà que cette équation peut s’écrire sous la forme
{ x′(t) = F (t, x(t)), t > 0 , x(0) = 1 ,
avec F : R^2 → R la fonction définie par F (t, x) = 5x^2 − x(1 + t^3 ). Plus généralement, étant donné F : R^2 → R , le théorème de Cauchy-Lipschitz assure que sous des conditions raisonnables, il existe une unique application x de classe C^1 sur [t 0 , tf ] (t 0 comme temps initial et tf comme temps final) vérifiant le problème de Cauchy :
(P )
x′(t) = F (t, x(t)), t > t 0 , x(t 0 ) = x 0 ,
L’objet des schémas numériques est d’obtenir des approximations de ces solutions dont la théorie donne l’existence mais ne dit pas comment les obtenir. Elles consistent en général à approximer la solution x en un certain nombre de points répartis sur [t 0 , tf ]
1.1 Méthode d’Euler explicite
L’idée principale est que «localement la courbe de la fonction x ressemble à sa tangente». Ainsi si h est proche de 0, on a
x′(t 0 ) ≈ x(t 0 + h) − x(t 0 ) h
ainsi x(t 0 + h) ≈ x(t 0 ) + hx′(t 0 ) = x(t 0 ) + hF (t 0 , x(t 0 )).
On peut donc approcher x(t 0 + h) par la quantité x(t 0 ) + hF (t 0 , x(t 0 ))
On découpe ainsi l’intervalle de temps [t 0 , tf ] en n segments de même longueur h = (tf^ − n t^0 ): le pas. On dispose ainsi de n + 1 temps ti = t 0 + ih pour i ∈ { 0 , ..., n}. On va alors approximer la solution x à l’instant ti par le nombre xi défini par la relation de récurrence :
xi+1 = xi + hF (ti, xi), i = 0, ..., n − 1 , t > t 0 , On initialise enfin avec la condition initiale x 0 = x(t 0 ). La méthode d’Euler explicite pour la résolution numérique de (P ) sur l’intervalle [t 0 , tf ] est définie par le shéma itératif suivant:
{ x 0 donné, t = t 0 xi+1 = xi + hF (ti, xi), i = 0, ..., n − 1 , t > t 0
Question: Écrire un code de résolution numérique de ce problème de Cauchy sur l’intervalle [t 0 , tf ] utilisant la méthode d’Euler explicite. Pour cela, écrire une fonction Euler_explicite dont les arguments d’entrée sont:
- la fonction F définissant le membre de droite de l’équation différentielle,
- la valeur initiale t 0 de la variable de temps t,
- la valeur finale tf de la variable de temps t,
- la valeur de la donnée initiale du problème x 0 ,
- le nombre n de morceaux dans le découpage de l’intervalle de temps considéré, et dont l’argument de sortie est le tableau des valeurs de la solution approchée. Programme Python pour la méthode d’Euler explicite:
[35]: def Euler_explicite(F,t0,tf,x0,n): h=(tf-t0)/n x=[x0] for i in range(n): x.append(x[i]+hF(t0+ih,x[i])) return x
Pour n = 10, on obtient ici une courbe approchée vraiment très mauvaise. Bien entendu, en prenant n = 1000 sur le même intervalle [0, 5], on a quelque chose de tout à fait satisfaisant.
[40]: # pour n= t=np.linspace(t0,tf,1001) #temps xapp=Euler_explicite(F,t0,tf,x0,1000) #solution approchée xex=np.exp(-3*t) #solution exacte plt.plot(t,xex,'r',t,xapp,'b--') plt.xlabel('temps (t)') plt.ylabel('x(t)') plt.legend(('solution exacte','méthode d ' Euler explicite')) plt.title("solution de $x'=-3x$ pour n=1000")
[40]: Text(0.5, 1.0, "solution de $x'=-3x$ pour n=1000")
Remarque: Plus le pas h = tf^ − n t^0 est petit, l’approximation sera meilleure (proche de la solution exacte).
2 Systèmes d’équations du premier ordre.
Nous souhaitons étudier l’évolution simultanée de deux populations x(t) et y(t) (disons pour fixer les idées que x représente des lapins et y des renards) soumises aux équations suivantes:
(S)
x′(t) = x(t)(3 − 2 y(t)), t > 0 , y′(t) = y(t)(−4 + x(t)), t > 0 ,
associé à la donnée initiale x(0) = x 0 et y(0) = y 0.
Ce genre de modélisation est connu sous le nom de système de Lotka-Volterra.
Questions:
- On suppose que initialement x 0 = 5 et y 0 = 3, calculer une solution approchée du système (S) sur l’intervalle de temps [0, 10] en utilisant la méthode d’Euler explicite (on se servira de la fonction Euler_explicite écrite dans la section précédente) et en prenant n = 1000.
- Représenter l’évolution des valeurs des approximations obtenues des x(t) et y(t) au cours du temps.
Réponse: Il n’est pas nécessaire de changer de méthode d’Euler sous prétexte qu’on a désormais un système d’équations à résoudre, le principe de la méthode d’Euler reste tout à fait valable.
3 Déjà disponible en Python
La fonction odeint(F,x0,t) nous permet d’obtenir une résolution numérique de référence pour l’équation différentielle, et fonctionne sur le meme mode que notre Euler, on lui donne:
- une fonction F (éventuellement vectorielle) représentant l’équation (ou le système),
- une condition initiale x 0 qui sera un vecteur dans le cas d’un système,
- un tableau de nombres correspondant aux valeurs ti pour lesquelles on souhaite calculer une valeur approchée de x(ti). Exemple : On considère l’équation différentielle
x′(t) = sin(t) sin(x(t)), t > 0 ,
associé à la donnée initiale x(0) = 1.
[44]: from scipy.integrate import odeint
[45]: F= lambda t, x:np.sin(t)*np.sin(x) t0= tf= x0= n= t=np.linspace(t0,tf,n+1) # temps x=odeint(F,x0,t) # solution x(t)
plt.plot(t,x,'r') plt.title("solution de x'=sin(t)sin(x)") plt.xlabel('temps (t)') plt.ylabel('x(t)') plt.show()
