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:
Les exemples:
Les fonctions sage suivantes sont abordés
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?
Par exemple: $y'+3y=\mathrm{e}^{x}$.
desolve(diff(y,t)+3*y==exp(t),y,show_method=True)
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)
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()
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})$.
r,t=var('r,t')
x = function('x')(t)
R_malthus(z)=r*z; R_malthus
Malthus=(diff(x,t)-R_malthus(x)==0); Malthus
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:
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))
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))
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.
r,K=var('r,K')
x = function('x')(t)
R_verhulst(z)=r*z-r*(z**2)/K; R_verhulst
Verhulst=(diff(x,t)-R_verhulst(x)==0);R_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]
Chercher les solutions d'états initiales respectifs:
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:
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)
On rappel la fonction de verhust et on verifie sa formule
R_verhulst
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
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
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()
t=var('t')
x=function('x')(t)
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
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
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?
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.
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:
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.