TP d'introduction à la résolution et à l'étude des équations différentielles¶

Auteur: slimane ben miled. Date: 9 février 2017

On ce propose dans ce TP d'introduire les utils de résolution et d'étude des EDO

On traite les notions:
- Solution d'équations différentielles
- Points d'équilibres
- Etude de stabilité
- Champs de vecteurs
Les exemples:
- Equation de Malthus, Verhulst
Les fonctions sage suivantes sont abordés
1. Rappel
  - plot
2. Nouvelles fonctions
  - *diff
  - desolve
  - desolve_rk4*
  - boucle: for

1. Introduction¶

On commence par faire un test pour voir si tout marche bien

1+2

On défini une variable $t$ et une fonction $y$ dépendant de cette variable:

t = var('t')
y = function('y')(t)

On utilise ensuite:

desolve(equation,variable,ics=...,ivar=...,show_method=... contrib_ode=...)

Faire un help sur desolve: desolve?

desolve?

Équations linéaires: Il s'agit d'équations du type: $$y'+P(x)y=Q(x)$$ avec $P$ et $Q$ des fonctions continues sur des intervalles donnés.

Par exemple: $y'+3y=\mathrm{e}^{x}$.

desolve(diff(y,t)+3*y==exp(t),y,show_method=True)

Équations à variables séparables: Il s'agit d'équations du type: $$P(x)=y'Q(y)$$ avec $P$ et $Q$ des fonctions continues sur des intervalles donnés.

Par exemple: $yy'=x$.

desolve(y*diff(y,t)==t,y,show_method=True)

Attention: Sage parfois ne reconnaît pas les équations à variables séparables et les traite comme des équations exactes.

Par exemple: $y'=\mathrm{e}^{x+y}$.

 desolve(diff(y,t)==exp(t+y),y,show_method=True)

2. Equations linéaires:¶

Résolvons par exemple $$y'+2y=t^2-2t+3$$

 desolve(diff(y,t)+2*y==t**2-2*t+3,y)

Ordonnons un peu tout ça avec la commande expand:

 desolve(diff(y,t)+2*y==t**2-2*t+3,y).expand()

Comment connaitre la mêthode utilisée par sage?

 desolve(diff(y,t)+2*y==t**2-2*t+3,y,show_method=True)[1]

Ajoutons des conditions initiales, par exemple $x(0)=1$:

 desolve(diff(y,t)+2*y==t**2-2*t+3,y,ics=[0,1]).expand()

3. Résolution de l'équation de Malthus¶

Le taux relatif de croissance d'une population est constant.

Pour l'étudier, on est amené à résoudre une équation du type~:

$$x'=rx$$

avec $r$ le taux d'accroissement de la population c'est un paramètre réel positif.

Remarque: A partir de maintenant la solution des équations différentielle se note $x\in C^1(\mathbb{R})$.

Définition des variables et des fonctions¶

r,t=var('r,t')

x = function('x')(t)

Definir la reponse fonctionelle de Malthus:¶

R_malthus(z)=r*z; R_malthus

Definir l'équation de malthus:¶

Malthus=(diff(x,t)-R_malthus(x)==0); Malthus

Recherche les solutions analytiques de l'équation différentielle:¶

On utilisera la fonction desolve. Les solutions sont stoquées dans une variable: sol

sol=desolve(Malthus,[x,t]); sol

Chercher les solutions d'états initiales:

$x(0)=1$
$x(0)=0$
$x(0)=2$

Utiliser l'option ics=$[t_0,x_0]$ de la fonction dsolve

sol1=desolve(Malthus,[x,t],ics=[0,1]);print(sol1)
sol0=desolve(Malthus,[x,t],ics=[0,0]);print(sol0)
sol2=desolve(Malthus,[x,t],ics=[0,2]);print(sol2)

Dessine la solution de l'équation de Malthus d'état initial $y(0)=1$ et pour $r=0.5$ et $t\in[0,1]$. Vous utiliserai la fonction plot.

plot?

sol=desolve(Malthus,[x,t],ics=[0,1]); sol
plot(sol.substitute(r==0.5) ,(0,1))

Etude de l'effet de $r$¶

Dessiner les solutions de l'équation de Malthus pour $r=0.5,1$ et $1.5$. Vous prendrer $x(0)=1$ comme etat initial pour $t\in [0,1]$.

plot([sol.substitute(r==rr) for rr in [-0.1..0.1, step=0.01]],(0,2))

4. Résolution du modèle de Verhust:¶

Le taux relatif de croissance d'une population est une fonction linéairement décroissante de la population.

Pour l'étudier, on peut être amené à résoudre une équation du type~:

$$x'=rx(1-x/K)$$

avec $r$ et $K$ des paramètres réels positifs.

Définir ce qui est une variable et ce qui est une fonction¶

r,K=var('r,K')

x = function('x')(t)

Definir la réponse fonctionnelle de Verlhust¶

R_verhulst(z)=r*z-r*(z**2)/K; R_verhulst

Verhulst=(diff(x,t)-R_verhulst(x)==0);R_verhulst

Résolution de analytique de l'EDO de Verhulst¶

solGen=desolve(Verhulst,[x,t]);solGen

Remarque: La solGen est une equation en $x(t)$. Essayons de la resoudre avec solve, la fonction qui résoud les équations algèbrique.

Commençons par la simplifier.

solGen.simplify_log()

Multiplions de part et autre par $r$

eq2=r*solGen.simplify_log();eq2

Résolvons maintenant l'équation eq2. Pour cela on utilisera la fonction solve:

eq3=solve(eq2,x(t));eq3

Attention Solve renvois un tableau avec les solutions

eq3[0]

Calcul des solutions avec différents états initiaux¶

Chercher les solutions d'états initiales respectifs:

$x(0)=1$
$x(0)=K$

sol1=desolve(Verhulst,[x,t],ics=[0,1]);print(sol1)
eq1=solve(r*sol1.simplify_log(),x(t));print(eq1)

Remarque: La solution est sous la forme d'un log(x(t)) et donc l'etat initiale $x(0)=0$ ne peux pas marcher ni même $x(0)=K$ (devinez pourqoi ?)

sol0=desolve(Verhulst,[x,t],ics=[0,0]);

solK=desolve(Verhulst,[x,t],ics=[0,K]); print(solK) 
#eqK=solve(r*solK.simplify_log(),x(t));print(eqK)

Dessiner les solutions de l'équation de Verhulst pour les états initiales:

$y(0)=1$

avec K=10 et $r=2$ pour $t\in [0,5]$

 plot(eq1[0].substitute(r==2,K==10).rhs() ,(0,10))

Dessiner les solutions de l'équation de Verhulst pour $r=-0.5,0,0.5$ avec $K=5,7,10$. Vous prendrer $x(0)=1$ comme etat initial pour $t\in [0,5]$.

g = Graphics()
for rr in [-0.5,0,0.5]:
    for kk in [5,7,10]:
         g+=plot(eq1[0].substitute(r==rr,K==kk).rhs(),(0,15))
show(g)

Calcul des points d'équilibres et étude de leurs stabilités¶

On rappel la fonction de verhust et on verifie sa formule

R_verhulst

Calcul des points fixes¶

Pour calculer les solutions de manière analytique on utilisera la fonction sage solve

x_eq=solve(R_verhulst(x)==x,x);x_eq

x_equilibre1=x_eq[0].right_hand_side();  x_equilibre1

x_equilibre2=x_eq[1].right_hand_side();  x_equilibre2

Etude de la stabilité des deux points d'équilibres¶

Pour cela on dois calculer la dérivé de R_verhulst, R_verhulst' en chaque point d'équilibre, x_equilibre1 et x_equilibre2 et le comparé sa valeur absolue à 1, pour toutes valeurs de $r$

Deriv_R_verhulst(z)=diff(R_verhulst(z),z);Deriv_R_verhulst

On doit construire l'application $r\mapsto $R_verhulst'(x_equilibre)

Deriv_R_verhulst(x_equilibre1)

Posons $f_1:r\mapsto $R_verhulst'(x_equilibre1)

f1(r)=Deriv_R_verhulst(x_equilibre1);f1

Faite la même chôse avec x_equilibre2.

Posons $f_2:r\mapsto $R_verhulst'(x_equilibre2)

f2(r)=Deriv_R_verhulst(x_equilibre2);f2

Dessiner $f_1$ et $f_2$ sur $\mathbb{R}$

G=plot(f1(r),(r,-1,2)) + plot(f2(r),(r,-1,2),color="red",title='$f_1(r)$ et $f_2(r)$')
G+=plot(1,(r,-1,2),thickness=3, color="green")
show(G)

Cherchons les valeurs de $r$ pour lesquels x_equilibre1 est stable (i.e. solutions $f_1(r)<1$)

EnsembleStable1=solve(f1(r)<1,r); EnsembleStable1

Même question avec x_equilibre2 (i.e. solutions $f_2(r)<1$)

EnsembleStable2=solve(f2(r)<1,r); EnsembleStable2

5. Champ de vecteur et solution numérique¶

Vous utiliserai la fonction desolve_rk4 qui permet de resoudre une équation différentielle en utilisant la mêthode de RK d'ordre 4. Vous utiliserai aussi l'option output='slop_filed' qui permet de dessiner le champs de vecteurs associé.

Refaire le même travail avec un solveur numérique, par exemple on prendra: desolve_rk4()

Vous dessinerez les solutions pour r=2 et K=10
Vous dessinerez le champs de vecteur associé.

Définition des variables¶

t=var('t')
x=function('x')(t)

Définition de la réponse fonctionelle de Verhulst¶

Remarque: Cette partie est inutile car la réponse fonctionelle a déja été définie plus haut. Cependant, on veut imposser a ce que $z\in \mathbb{R}$

R_verhulst(z)=(2.*z*(1-z/10));R_verhulst

Définition de EDO de Verhulst¶

Remarque: Cette partie est inutile car la réponse fonctionelle a déja été définie plus haut.

Verhulst_eq = (diff(x, t) == R_verhulst(x));Verhulst_eq

Résolution numérique de l'équation de Verhulst¶

On utilisera l'integrateur numérique desolve_rk4 basé sur la mêthode numérique Range Kutta d'ordre 4 (voir votre cours d'analyse numérique).

desolve_rk4(Verhulst_eq,x, ics=[0,10],step=1, end_points=40)

tx,ty = var('tx,ty')
desolve_rk4(R_verhulst(ty),ty, ics=[0,12],ivar=tx,output='slope_field',end_points=[0,6],thickness=3)

Exercice 1: Modèle avec pêche/chasse:

Modifiez le modèle de compétition intra-specifique de Verlhust pour inclure l'idée suivante: A chaque pas de temps, on introduit un individu. Résoudre le modèle, imprimer le graphique et décrire ce qui se passe avec la capacité limite.
En utilisant l'idée précédente, combien d'individu peuvent-ils être récoltés à chaque instant avant extinction?
Comment changer le modèle, si le nombre d'individu pêchés/chassés est proportionnel au nombre total d'individu présents? (exemple: on pêchera/chassera 25% des individu à chaque moment.)

Exercice 2: La NASA s'intéresse à la dynamique de population d'une base lunaire. Répondez aux questions ci-dessous en vous basant sur un modèle approprié pour la croissance démographique. Sauvegardez vos réponses en utilisant les graphiques appropriés.

Supposons que la population maximale durable sur la base est de 100 000 personnes, le taux de natalité dans la base par an est de 11/10 000 et la population initiale est de 30 000.

Comment la population change-t-elle au cours des 15 premières années d'activité? Quand la population at-elle atteint 90% de ca capacité limite ?
La NASA prévoit d'envoyer environ 500 «immigrants» par an. Comment pourriez-vous modeler cela? Comment cela affecte-t-il la population? À l'aide de votre modèle, quant est ce que 90% de la population maximale sera atteinte ?

Exercice 3: Une espèce de grue des sables d'Amérique du Nord (Grus canadensis) a été déclarée espèce en voie de disparition en 1916. Ces oiseaux causent des dommages considérables aux cultures lorsque leur nombre augmente, de sorte que les permis de chasse ont été autorisés à partir de 1961. La principale question d'intérêt est: Combien de permis devront être vendus par an? Nous avons les données suivantes de R. Miller et D. Botkin sur la base de leur étude démographique de 1974:

population maximale viable est de $194 600$
population actuelle est de $190 000$
Increasing rate $5.07\ 10^{-7}$ (par unité de temps [année])

Nous voulons construire un modèle qui intègre ces données avec un terme supplémentaire pour la limite de chasse (en supposant qu'un nombre constant sont tués chaque année). Nous voulons choisir ce nombre pour que l'espèce ne disparaisse pas. Trouvez un bon nombre de permis de chasse.