Algorithme de Runge Kutta -- besoin d'exercice

**kromartien** · 24/12/2007, 18h38

Bonjour, j'ai étudié l'algorithme de Runge Kutta de résolution d'équations différentielles, et j'ai trouvé que :

Soit

f(t,y)=y', l'équation différentielle régissant y.
y(0)=y0, la condition initiale sur y.
f(t,y0)=y'0, la condition initiale sur y'.

On peut écrire l'algorithme de Runge Kutta. Est-ce bien certain qu'il faille connaître toutes ces données pour pouvoir écrire l'algorithme de Runge Kutta ?
En particulier, je m'interroge sur la fonction dérivée fonction de y et de t.

Avez-vous des exemples concrets d'équations différentielles existant dans la nature pouvant être calculée d'après la méthode de Runge Kutta ?
J'en voudrai bien une ou deux pour pouvoir écrire un code, afin de mieux en comprendre le fonctionnement. Sont-ce bien des équations différentielles du premier ordre non linéaires, par exemple du type y' = t * y(t) ?

Merci beaucoup.

**kromartien** · 24/12/2007, 23h10

Alors j'ai réalisé l'exercice pour l'équation différentielle
y' + y/x = 0, qui correspond bien sûr à la fonction y=1/x, et j'ai obtenu les résultats suivants

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
 
#Temps			y(t)
1.000000		1.000000
1.000001		1.000001
1.000002		0.999998
1.000004		0.999996
1.000008		0.999992
1.000016		0.999984
1.000032		0.999968
1.000064		0.999936
1.000128		0.999872
1.000256		0.999744
1.000512		0.999488
1.001024		0.998977
1.002048		0.997955
1.004096		0.995917
1.008192		0.991858
1.016384		0.983816
1.032768		0.968024
1.065536		0.937582
1.131072		0.881008
1.262144		0.783254
1.524288		0.636937
2.048576		0.470839
3.097152		0.352208
5.194304		0.278881
9.388608		0.244652
17.777216		0.209791

avec l'algorithme suivant (langage C)

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
64
65
66
67
68
69
70
71
72
73
74
75
76
 
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
 
 
double calc_next_step(long, double, double, double);
double calc_y(double, double);
 
int main(void)
{
    long indice;
    double min_pas;
    double max_time;
    double t_0,y_0,y_prim_0;
    double t_i,y_i,y_prim_i;
    FILE * fp_out;
    fp_out = fopen("data_out.dat","w");
    if(fp_out == NULL)
    {
        puts("erreur d'ouverture de fichier");
        return 1;
    }
    max_time = 10;
	min_pas = 0.000001;
	/*************************
	Conditions initiales
	*************************/
	t_0 = 1;
	y_0 = 1;
	y_prim_0 = 1;
	/************************
	initialisation des variables
	************************/
	y_prim_i = y_prim_0;
	t_i = t_0;
	indice=0;
	fprintf(fp_out,"#Temps\t\ty(t)\n");
	fprintf(fp_out,"%lf\t\t%lf\n",t_0,y_0);
	/******************
        boucle d'iteration
	******************/
	while(t_i<max_time)
	{
	    /*calcul du y_i*/
	    y_i=calc_next_step(indice, min_pas, y_prim_i, y_0);
	    /*incrementation du temps courant*/
            t_i = t_0 + min_pas * pow(2,indice);
	    /*calcul de y_prim_i a partir de y_i et de t_i*/
	    y_prim_i = calc_y_prim(y_i,t_i);
            /*ecriture dans le fichier de sortie*/
	    fprintf(fp_out,"%lf\t\t%lf\n",t_i,y_i);
	    /*incrementation de l'indice de boucle*/
	    indice = indice+1;
	}
	fclose(fp_out),fp_out = NULL;
	puts("fin du programme");
	puts("appuyez sur une touche pour continuer");
	getchar();
    return 0;
}
 
double calc_y_prim(double y_i, double t_i)
{
    double y_prim_i;
    y_prim_i = -(1.0/t_i) * y_i ;
    return y_prim_i;
}
 
 
double calc_y(long i, double min_pas, double y_prim_i, double y_0)
{
    double y_i;
    y_i = y_0 + min_pas * pow(2, i) * y_prim_i;
    return y_i;
}

Pouvez-vous me dire si ce code correspond à une bonne utilisation de la méhode de Runge Kutta ? Merci beaucoup.

**FR119492** · 25/12/2007, 19h40

Salut !

Il existe diverses méthodes pour intégrer les équations différentielles ordinaires avec conditions initiales, et même plusieurs méthodes de Runge-Kutta. La plus utilisée est celle du 4ème ordre, appelée aussi RK4. Elle s'applique aussi bien à un vecteur d'état à plusieurs composantes qu'à une seule variable d'état.

Il n'est pas nécessaire de reprogrammer la méthode elle-même pour chaque application: on la stocke sous la forme d'un sous-programme RK4 dans une bibliothèque, et pour chaque nouvelle application, on écrit un sous-programme qui, connaissant le temps t et les valeurs y des diverses variables d'état, calcule les valeurs de leurs dérivées, et un programme principal qui fixe les conditions initiales et parcourt une boucle dans laquelle le sous-programme RK4 est appelé à chaque passage.

Il n'est donc pas nécessaire d'avoir les valeurs initiales des dérivées. Je te joins un sous-programme appliquant la méthode RK4:

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
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
 
      Subroutine D006(SM,N,T,DT,Y)
C
C***********************************************************************
C
C     Bibliothèque JMBPACK
C
C     Série D: Equations différentielles
C
C***********************************************************************
C
C     Méthode de Runge-Kutta d'ordre 4 (1/6 + 2/6 + 2/6 + 1/6)
C
C***********************************************************************
C
C     Version 3.0: Jean-Marc Blanc, juillet 2004
C
C **********************************************************************
C
C     SM     Nom d'un sous-programme fourni par l'utilisateur pour le
C            calcul des dérivées à partir des valeurs de la variable
C            indépendante et des variables d'état. Ce sous-programme
C            doit être de la forme
C              Subroutine SM(N,T,Y,DY)
C                DY(1)= ...
C                DY(2)= ...
C                 .
C                 .
C                DY(N)= ...
C                Return
C              End
C            et avoir été déclaré dans une instruction External.
C
C     N      Ordre du système différentiel.
C
C     T      Variable indépendante.
C
C     DT     Pas d'intégration.
C
C     Y      Vecteur des variables d'état.
C
C **********************************************************************
C
      Implicit None
C
      Integer N
      Real*8 DT,T,Y(N)
C
      Integer I
      Real*8 DY(N),K1(N),K2(N),K3(N),K4(N),Z(N)
C
      Call SM(N,T,Y,DY)
      Do I=1,N
        K1(I)=DT*DY(I)
        Z(I)=Y(I)+K1(I)/2.d0
      End Do
C
      Call SM(N,T+DT/2.d0,Z,DY)
      Do I=1,N
        K2(I)=DT*DY(I)
        Z(I)=Y(I)+K2(I)/2.d0
      End Do
C
      Call SM(N,T+DT/2.d0,Z,DY)
      Do I=1,N
        K3(I)=DT*DY(I)
        Z(I)=Y(I)+K3(I)
      End Do
C
      Call SM(N,T+DT,Z,DY)
      Do I=1,N
        K4(I)=DT*DY(I)
        Y(I)=Y(I)+(K1(I)+2.d0*K2(I)+2.d0*K3(I)+K4(I))/6.d0
      End Do
C
      T=T+DT
      Return
C
      End

Comme exercices, je te suggère

y'=a*y y(0)=1 a>0

y'=a*y y(0)=1 a<0

y1'=y2 y1(0)=0
y2'=-y1 y2(0)=1

Jean-Marc Blanc

**kromartien** · 25/12/2007, 20h35

Bonjour. Merci beaucoup pour la réponse et vote code. En réalité, j'ai vu la méthode dite RK4 sur Wikipédia, mais je voulais savoir programmer RK1 pour mieux en comprendre le fonctionnement. Merci malgré tout.

D'après vos exemples, a est une constante, on a donc y'(x) qui est fonction uniquement de y, ce sont des équations différentielles à coefficients constants. À ce moment, la méthode d'Euler par propagation de proche en proche semble tout à fait indiquée.

La méthode de Runge Kutta n'est-elle pas utile seulement pour des équations différentielles à coefficients non constants ?

**FR119492** · 25/12/2007, 21h01

Salut !

Si tu veux comprendre le fonctionnement d'Euler ou RK1, plutôt que de programmer, tu devrais prendre une grande feuille de papier millimétré et faire 2 ou 3 pas à la main, ce qui te visualisera bien mieux la méthode.

Jean-Marc Blanc