Discussion :

[lucienGr]

Bonsoir,

Dans le cadre d'un projet concernant l'analyse numérique des équations différentielles ordinaires, il m'est demandé de résoudre le système de Lotka-Volterra (système proie-prédateur) et ensuite de représenter la fonction solution en lui appliquant différents schémas numériques (Euler progressif, Euler rétrograde, Crank-Nicolson, Runge-Kutta d'ordre 4...)

J'ai donc calculé la solution qui se trouve être de la forme
H(t)=kx(t)-kpln(x(t))+y(t)-ln(y(t)) (avec k et p deux constantes strictement positives)

Il me faut donc maintenant tracer cette fonction grâce aux différents schémas numériques. Le problème c'est que je n'ai découvert Matlab que cette année et que je suis vraiment débutant. Je pense donc qu'il y a des fonctions existantes qui doivent m'aider à résoudre mon problème mais je ne les connais pas et ne sait pas comment les appliquer.

En vous remerciant d'avance pour votre aide

[Dam2227]

Salut

Tu peux regarder du côté des fonctions ode. Il suffit de taper ode dans le help de matlab. Il y des exemples qui t'expliqueront comment les utiliser.

++

[lucienGr]

Merci à toi de ta réponse

En cherchant un peu, j'ai trouvé les éléments nécessaires pour calculer la solution d'une équation différentielle par les méthodes numériques d'Euler. Je fais comme suit avec l'équation y'(t)=-t*y(t)^2

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
f=’-t*y.ˆ2’;
h=0.2; t0=0; tf=4; y0=2;
[t_EP, y_EP] = euler_prog(t0,y0,tf,h,f)
 
df=’-2*t*y’;
[t_ER, y_ER] = euler_retr(t0,y0,tf,h,f,df)
 
t=[0:0.05:4]; sol_ex=’2./(1+t.ˆ2)’;
plot(t,eval(sol_ex),’b’,t_EP,y_EP,’ro-’,t_ER, y_ER’,’go-’)
legend(’sol-ex’,’Euler prog.’,’Euler retr’)

Je peux donc ainsi comparer sur un même graphe la solution exacte et les solutions calculées.

Mon problème est que le système de Lotka-Volterra est un système de deux équations différentielles croisées.
Quelqu'un pourrait-il donc me fournir les commandes nécessaires à la résolution d'un système d'équations différentielles croisées car j'ai beau chercher, je ne trouve pas les commandes nécessaires

En vous remerciant d'avance

[Dam2227]

J'imagine que ce que tu entends par équations différentielles croisées, c'est en fait que tes équations sont non-linéaires : si tu as une variable x et une variable y, alors dans tes équations, tu as un truc du style x*y. En gros, tu ne peux pas résoudre séparément x et y chacun de leur côté puisqu'ils dépendent l'un de l'autre. Si ce n'est pas ça, qu'est-ce que tu entends par équations différentielles croisées?

Sinon, pour répondre à ta question, l'exemple traité dans Matlab inclu exactement ce que tu cherches (enfin ce que je pense que tu cherches^^). Donc tape ode dans le help de maltab et regarde les exemples.

++

[lucienGr]

Bonsoir, merci de tes réponses.

Je me suis donc mis à la tache (je suis en fait passé à Scilab car ma licence Matlab a expiré et il était plus facile de passer sous Scilab)

Bref j'arrive à appliquer les méthodes de Euler et de Runge-Kutta

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
 
xbasc();
 
// Définition du pas de temps
 
h = 0.01;
Tmax=100;
Taille = 10*Tmax;
Nbintervalle=Tmax/h;
 
// Vecteur des abscisses
 
T = [1:Taille];
 
T = 10*h*T;
T = T';
 
// Initialisation des vecteurs des populations
 
for I=1:Taille,
 X(I) = 0;
 Y(I) = 0;
end
 
// Constantes de dynamique intrinsèque et populations initiales
 
alpha1=1;
alpha2=1;
x=2;
y=2;
I=0;
 
X(1)=x;
Y(1)=y;
 
// Calcul de l'évolution des populations
 
for I=1:Nbintervalle-1,
        I = I + 1;
        dx = alpha1*x - x*y;
        dy = - alpha2*y + x*y;
        x = x + h*dx;
        y = y + h*dy;
 
if modulo(I,10)==0 then
         X(I/10) = x;
        Y(I/10) = y;
end
 
end
 
// Affichage
 
xbasc();
 
plot2d([T T],[X Y],[2,3],"111","Proies@Prédateurs",[0,0,Tmax,4],[2,10,2,10]);
 
halt();
xbasc();
 
// Affichage de la trajectoire dans le plan de phase
 
plot2d([X],[Y],[2],"111","Trajectoire des phases",[0,0,4,4],[2,10,2,10]);

Après vérification, ce code fonctionne très bien

En revanche, lorsque je veux appliquer la méthode de Crank-Nicolson, cette dernière n'étant pas explicites mais implicites, il me faut calculer le terme f(n+1) pour l'injecter dans ma formule et la je ne sais pas du tout comment faire.

Quelqu'un pourrait-il m'aider?

En vous remerciant

[Dam2227]

Salut, je ne sais pas si tu as vu cette page internet :
http://en.wikipedia.org/wiki/Crank%E...icolson_method
C'est en anglais, mais bon. Il y a un exemple clair et précis qui accompagne l'explication...

++

[salseropom]

En revanche, lorsque je veux appliquer la méthode de Crank-Nicolson, cette dernière n'étant pas explicites mais implicites, il me faut calculer le terme f(n+1) pour l'injecter dans ma formule et la je ne sais pas du tout comment faire.

Quelqu'un pourrait-il m'aider?

Salut, l'intérêt des schéma numériques implicites par rapport aux schémas explicites c'est qu'ils sont plus stables. Ils ont un inconvénient : il faut résoudre à un moment donné une équation du type F(y) = 0. En multi-d tu n'as qu'une seule solution : la méthode de Newton (et toutes ses dérivées)

exemple pour Euler implicte :

y(n+1) = h_n * f(t(n+1),y(n+1))

soit Y = y(n+1) et F(Y) = Y - h_n * f(t(n+1),Y)

tu résous F(Y) = 0 et c'est gagné