CoCalc -- Systeme.sagews

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 $\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é $x_2$ de la première équation puis remplacé, et ainsi obtenu une équation d'ordre 2, linéaire à coefficients constants en $x_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

[

\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)}

\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 $x_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,x_1(t), x_2(t))$ décrit une courbe dans l'espace. À $t=0$ elle commence au point correspondant aux conditions initiales, puis, à mesure que $t$ tend vers $\infty$ , les composantes $x_1$ et $x_2$ tendent vers $0$

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 $x_1x_2$ , on laisse tomber, dans un certain sens, la variable $t$ . On peut aussi d'intéresser au graphique dans le plan de phase, c'est à dire le plan $x_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 $(x_1,x_2)$ du plan on dessine un vecteur collinéaire à $(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$'])

(

\displaystyle x_{1}

\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)$ désigne le nombre de proies, et $y(t)$ le nombre de prédateurs, le modèle de Lotka - Volterra s'écrit: $\left\{ \begin{array}{ccc}x' &= &ax-bxy\\ y' & = & dxy-cy\end{array}\right.$ où $a,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)$ et $y(t)$ . Nous avons montré que $\frac{y^a}{e^{by}}\frac{x^c}{e^{dx}} = K$ où $K$ est donnée par les conditions initiales. Cette équation décrit une courbe. Ci-bas on la trace en prennant $x_0 = 2.2$ et $y(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)$ et $y(t)$ en fonction de $t$ , on peut faire une courbe paramétrique $\{(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)]$ . En manipulant ces listes on peut produire les courbes $(t,x(t))$ , $(t,y(t))$ , ou encore, dans le plan de phase, la courbe $(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︠