Résolution numérique d'une équation différentielle d'ordre 1

Nous verrons ici quelques techniques de résolution numérique d'équations différentielles pour des équations d'ordre 1 . Toutes les méthodes classiques sont implantées dans Maple, mais on peut également définir soi-même la procédure qui nous intéresse .

Nous résoudrons ici l'équation suivante y' = x + y avec la condition initiale y( 0 ) = 1

On cherche à déterminer y(1) ou à estimer une courbe solution quand x varie de 0 à 1.

Plutôt que d'utiliser les routines déjà implantées dans Maple, on programmera ici les procédures permettant d'utiliser les méthodes d'Euler, d'Euler amélioré et de Runge-Kutta.

Méthode d'Euler

La procédure suivante appliquera la méthode d'euler avec n étapes pour estimer la solution de l'équation

y' = f( x , y) en x =b en connaissant la valeur initiale y ( a ) = c . La variable ynew représente la nouvelle valeur de y à chaque itération.

> Euler:=proc(f,a,c,b,n)
local h,i,x,y,ynew;
  h:=(b-a)/n:
  x:=a: y:=c:
for i from 1 to n do
  ynew:=y+h*f(x,y):
  y:=ynew:
  x:=a+i*h:
end do:
[x,evalf(y)]:end;

Définissons la fonction f( x , y) à être utiliser par notre procédure et appliquons-la!

> f:=(x,y)->x+y;

> Euler(f,0,1,1,10);

Si on désire voir tous les couples de points [x , y] calculés par la méthode, on modifie la procédure comme suit:

> Euler_detail:=proc(f,a,c,b,n)
local h,i,x,y,ynew;
  h:=(b-a)/n:
  x:=a: y:=c: print(evalf(a),evalf(c)):
for i from 1 to n do
  ynew:=y+h*f(x,y):
  y:=ynew:
  x:=a+i*h:print(evalf(x),evalf(y)):
end do:
end;

> Euler_detail(f,0,1,1,10);

Méthode d'Euler amlioré

Voici le code pour afficher seulement la valeur finale cherchée, toujours avec le même exemple.

> EulerAmel:=proc(f,a,c,b,n)
local h,i,x,y,ynew,k1,k2;
  h:=(b-a)/n:
  x:=a: y:=c:
for i from 1 to n do
   k1:=f(x,y):
   k2:=f(x+h,y+h*k1):
  ynew:=y+h*(k1+k2)/2:
  y:=ynew:
  x:=a+i*h:
end do:
[x,evalf(y)]:end;

Exécutons cette procédure:

> EulerAmel(f,0,1,1,10);

Si on désire voir tous les couples de points [x , y] calculés avec la méthode d'Euler amélioré, on utilise:

> EulerA_detail:=proc(f,a,c,b,n)
local h,i,x,y,ynew,k1,k2;
  h:=(b-a)/n:
  x:=a: y:=c: print(evalf(a),evalf(c)):
for i from 1 to n do
   k1:=f(x,y):
   k2:=f(x+h,y+h*k1):
  ynew:=y+h*(k1+k2)/2:
  y:=ynew:
  x:=a+i*h:print(evalf(x),evalf(y)):
end do:
end;

Exécutons cette procédure:

> EulerA_detail(f,0,1,1,10);

Méthode de Runge-Kutta

Voici le code pour afficher seulement la valeur finale cherchée, toujours avec le même exemple.

> RungeK:=proc(f,a,c,b,n)
local h,i,x,y,ynew,k1,k2,k3,k4;
  h:=(b-a)/n:
  x:=a: y:=c:
for i from 1 to n do
   k1:=h*f(x,y):
   k2:=h*f(x+h/2,y+k1/2):
   k3:=h*f(x+h/2,y+k2/2):
   k4:=h*f(x+h,y+k3):
  ynew:=y+(k1+2*k2+2*k3+k4)/6:
  y:=ynew:
  x:=a+i*h:
end do:
[x,evalf(y)]:end;

Exécutons cette procédure:

> RungeK(f,0,1,1,10);

Si on désire voir tous les couples de points [x , y] calculés avec la méthode Runge-Kutta d'ordre 4, on utilise:

> RungeK_detail:=proc(f,a,c,b,n)
local h,i,x,y,ynew,k1,k2,k3,k4;
  h:=(b-a)/n:
  x:=a: y:=c: print(evalf(a),evalf(c)):
for i from 1 to n do    k1:=h*f(x,y):    k2:=h*f(x+h/2,y+k1/2):    k3:=h*f(x+h/2,y+k2/2):    k4:=h*f(x+h,y+k3):
  ynew:=y+(k1+2*k2+2*k3+k4)/6:
  y:=ynew:
  x:=a+i*h: print(evalf(x),evalf(y)):
end do:
end;

Exécutons cette procédure:

> RungeK_detail(f,0,1,1,10);

>