CoCalc Public FilesSysteme.sagewsOpen with one click!
Author: Juan Carlos Bustamante
Views : 68
Compute Environment: Ubuntu 18.04 (Deprecated)

Deux réservoirs couplés.

Le problème traité en cours mène à un système de deux équations différentielles linéaires à coefficients constants, à savoir {x1=110x1+340x2x2=110x115x2\left\{\begin{array}{ccr}x'_1 & = & -\frac{1}{10}x_1 +\frac{3}{40}x_2\\ x'_2 & = &\frac{1}{10}x_1 -\frac{1}{5}x_2\end{array} \right.

Nous avons résolu le système par substitution, nous avons isolé x2x_2 de la première équation puis remplacé, et ainsi obtenu une équation d'ordre 2, linéaire à coefficients constants en x1x_1. Cette équation se résoud facilement (fait en classe).

À l'ordinateur on peut résoudre le système comme suit :

typeset_mode(True) t = var('t') x1 = function('x1')(t) x2 = function('x2')(t) de1 = diff(x1,t) + 1/10*x1-3/40*x2 == 0 de2 = diff(x2,t) - 1/10*x1 +1/5*x2 == 0 Sols = desolve_system([de1, de2], [x1,x2], ics=[0,-17,-21]) Sols
[x1(t)=1658e(120t)+298e(14t)\displaystyle x_{1}\left(t\right) = -\frac{165}{8} \, e^{\left(-\frac{1}{20} \, t\right)} + \frac{29}{8} \, e^{\left(-\frac{1}{4} \, t\right)}, x2(t)=554e(120t)294e(14t)\displaystyle x_{2}\left(t\right) = -\frac{55}{4} \, e^{\left(-\frac{1}{20} \, t\right)} - \frac{29}{4} \, e^{\left(-\frac{1}{4} \, t\right)}]

**Note : ** Sols est une liste, le zéroième élément est x1x_1, qui est une expression, pour accéder à la fonction elle-même, il faut prendre le coté droit de l'expression, ce qui se fait avec rhs

x1 = Sols[0].rhs() x2 = Sols[1].rhs() C1 = plot(x1,t,0,50,legend_label="$x_1$") C2 = plot(x2,t,0,50, legend_label="$x_2$",color="red") show(C1+C2, axes_labels=['$t$','$x$'])

On peut penser que (t,x1(t),x2(t))(t,x_1(t), x_2(t)) décrit une courbe dans l'espace. À t=0t=0 elle commence au point correspondant aux conditions initiales, puis, à mesure que tt tend vers \infty, les composantes x1x_1 et x2x_2 tendent vers 00

Courbe = parametric_plot3d([t,x1(t),x2(t)],(t,0,100), color='red') show(Courbe, aspect_ratio=[1,5,5])
3D rendering not yet implemented

Le plan de phase

Il s'agit du plan x1x2x_1x_2, on laisse tomber, dans un certain sens, la variable tt. On peut aussi d'intéresser au graphique dans le plan de phase, c'est à dire le plan x1x2x_1x_2. En somme, il s'agit de la projection de la courbe ci-haut.

C = parametric_plot([x1(t),x2(t)], (t,0,50), color="blue") show(C, aspect_ratio=1/5, axes_labels=['$x_1$','$x_2$'])

Plus illustratif encore : on peut tracer le champ de vecteurs donné par le système. En chaque point (x1,x2)(x_1,x_2) du plan on dessine un vecteur collinéaire à (x1,x2)(x'_1, x'_2). Il indique la direction de la courbe tangente en chaque point.

var("x1,x2"); F=plot_vector_field([-1/10*x1 +3/40*x2,1/10*x1-1/5*x2],(x1,-20,0),(x2,-20,0), color="red") show(F, axes_labels=['$x_1$','$x_2$'])
(x1\displaystyle x_{1}, x2\displaystyle x_{2})

Mais ce qu'il y a de mieux, c'est de considérer la courbe solution avec le champ.

show(F+C, axes_labels=['$x_1$','$x_2$'], aspect_ratio=1/2)

Des requins et des sardines.

On a considéré le modèle qui décrit un système proie - prédateur. Si x(t)x(t) désigne le nombre de proies, et y(t)y(t) le nombre de prédateurs, le modèle de Lotka - Volterra s'écrit: {x=axbxyy=dxycy\left\{ \begin{array}{ccc}x' &= &ax-bxy\\ y' & = & dxy-cy\end{array}\right.a,b,c,d>0a,b,c,d>0.

Il s'agit d'un système non-linéaire, mais on peut quand même dire quelque chose au sujet de la relatin entre x(t)x(t) et y(t)y(t). Nous avons montré que yaebyxcedx=K\frac{y^a}{e^{by}}\frac{x^c}{e^{dx}} = KKK est donnée par les conditions initiales. Cette équation décrit une courbe. Ci-bas on la trace en prennant x0=2.2x_0 = 2.2 et y(0)=5.5y(0) = 5.5

a=3 b=1 c=1 d=0.1

Le plan de phase

Plutôt que dessiner les deux fonctions x(t)x(t) et y(t)y(t) en fonction de tt, on peut faire une courbe paramétrique {(x(t),y(t))tR}\{(x(t),y(t))| t\in \mathbb{R}\}. C'est le diagramme de phase.

var('x,y') C= implicit_plot((y^a)*(x^c)*e^(-d*x-b*y)-K, (x,0,35),(y,0,8),color="blue") show(C,aspect_ratio =3,figsize=7, axes_labels=['Proies','Predateurs'] )
(x, y)

Comme avant, on peut aussi tracer le champ de vecteurs, ça nous dit en quel sens la courbe est parcourue.

F=plot_vector_field([a*x - b *x*y, -c*y + d*x*y],(x,0,35),(y,0,8), color="red") show(F+C, aspect_ratio =3,figsize=7, axes_labels=['Proies','Predateurs'] )

Le système est non linéaire, de sorte qu'on ne peut pas le résoudre (facilement). Par contre, on peut demander à SAGE de le faire numériquement. La commande desolve_system_rk4 fait el boulot. Ici rk4 dit que c'est la méthode de Runge-Kutta d'ordre 4 qui est utilisée, chose à voir dans vos cours de méthodes numériques.

Cette commande produit une liste qui, dans notre cas contient des triplets [t,x(t),y(t)][t,x(t),y(t)]. En manipulant ces listes on peut produire les courbes (t,x(t))(t,x(t)), (t,y(t))(t,y(t)), ou encore, dans le plan de phase, la courbe (x(t),y(t))(x(t),y(t)).

from sage.calculus.desolvers import desolve_system_rk4 x,y,t=var('x y t') P=desolve_system_rk4([a*x - b *x*y, -c*y + d*x*y],[x,y],ics=[0,0.5,2],ivar=t,step=0.05,end_points=10) X=[ [i,j] for i,j,k in P] LX=list_plot(X, plotjoined=True, color="red",legend_label="proies") Y=[ [i,k] for i,j,k in P] LY=list_plot(Y, plotjoined=True,color="green", legend_label="predateurs", axes_labels=['$t$',' ']) show(LX + LY)
Phase = [ [j,k] for i,j,k in P] LP = list_plot(Phase, plotjoined=True) show(LP)
8034ee0f-9496-4a50-936c-e2aa9302315a a18a4415-fecc-4eb2-b154-e020eab93078