Méthode Runge Kutta d'ordre 4 en PASCAL

[Boo-Bouh]

Bonsoir à tous !
Depuis quelques temps je me suis mis aux méthodes de résolutions d'equo diff de façon numérique et evidemment je suis tombé sur la méthode de Runge kutta d'ordre 4 qui permet une approximation assez précise de la solution a l'aide de la formule :

yn+1=yn+ h/6(k1+2k2+2k3+4)

Je cherche actuellement un algorithme en PASCAL qui permettrait d'avoir en entrée une equo diff du premier ordre et en sortie un tableau avec comme premiere colonne le temps qui doit suivre la relation tn+1=t+h et comme autre colonne les images calculees a l 'aide de l'algo(en gros les images de la fonction solution).

Connaissez vous donc un site ou je pourrais trouver cet algorithme?
Mes connaissances en Pascal sont assez faibles,malgré la simplificité du programme j'eprouve beaucoup de difficultés a écrire cet algo.


Merci d'avance.

[FR119492]

Salut!
Je l'ai fait en Fortran et non en Pascal, mais je pense que ça peut te donner des idées. Le programme se composera toujours de trois morceaux:

1. Un programme principal que tu devras écrire toi-même pour chaque nouveau problème, qui comportera une boucle pour calculer les valeurs successives en appelant à chaque pas le sous-programme D006.
2. Le sous-programme (ici D006) qui contient la méthode de Runge-Kutta d'ordre 4, qui est programmé une fois pour toutes et que tu pourras réutiliser chaque fois que tu voudras intégrer des équations différentielles.
3. Un sous-programme que tu devras écrire de cas en cas et qui contiendra les équations différentielles.

Code :

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
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
C
      Subroutine D006(SM,N,T,DT,Y,KErr)
C
C***********************************************************************
C
C     Bibliothèque JMBPACK
C
C     Série D: Equations différentielles
C
C***********************************************************************
C
C     Routine D006: Méthode de Runge-Kutta d'ordre 4
C                   (1/6 + 2/6 + 2/6 + 1/6)
C
C     Version 1.0, Jean-Marc Blanc, Novembre 2011
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     A l'entrée:

C     N      Ordre du système différentiel.
C
C     T      Variable indépendante au début du pas.
C
C     DT     Pas d'intégration.
C
C     Y      Vecteur des variables d'état au début du pas.
C
C     A la sortie:
C
C     T      Variable indépendante à la fin du pas.
C
C     Y      Vecteur des variables d'état à la fin du pas.
C
C     KErr   Message d'erreur:
C              0  L'intégration s'est effectuée correctement.
C              1  N non strictement positif.
C              2  DT non strictement positif.
C
C **********************************************************************
C
      Implicit None
C
      External SM
      Integer N,KErr
      Real*8 DT,T,Y(N)
C
      Integer I
      Real*8 DY(N),K1(N),K2(N),K3(N),K4(N),Z(N)
C
C **********************************************************************
C
      KErr=0
      If (N.Le.0) Then
        KErr=1
        Call U004('N non strictement positif')
        Return
      End If
C
      If (DT.Le.0.d0) Then
        KErr=2
        Call U004('DT non strictement positif')
        Return
      End If
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

Jean-Marc Blanc

[Boo-Bouh]

Tout d'abord merci pour cette réponse rapide.
J'ai analysé votre algo bien que je ne sache pas du tout coder en Fortran (je commence a peine la programmation donc je connais que le pseudo-code et un peu de pascal).
Est-il compliqué de passer du frotran au pascal? Parcequ'a vu d'oeil ça me dit vraiment rien tout ça.

Merci de votre patience

[Boo-Bouh]

Le code en Pascal est il construit de la même façon?
Parceque moi je voyais plutot:
-Declaration des tableaux de valeurs (le temps et les images)
-Une boucle qui calcule les valeurs du temps (avec tn+1=t+h)
-Le corps du programme

[pseudocode]

Envoyé par Boo-Bouh

Le code en Pascal est il construit de la même façon?
Parceque moi je voyais plutot:
-Declaration des tableaux de valeurs (le temps et les images)
-Une boucle qui calcule les valeurs du temps (avec tn+1=t+h)
-Le corps du programme

il y a une version Fortran plus compact sur netlib.

Sinon il y a des exemples en C++ et en Java disponibles dans ce forum.

[Boo-Bouh]

Merci.
Je vais donc essayer d'ecrire l'algo en Pascal mais ce qui me fait défaut c'est plutôt la syntaxe a utiliser.

[Boo-Bouh]

Bonsoir, j'ai fait un algo sur la méthode d'euleur que j'ai trouvé beaucoup simple.
Comment faire maintenant pour passer à la méthode RK4?
Il suffit que je rajoute des nouvelles variables mais je bloque un peu.
Voici le code:

Code :

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
program euler
type Tableau= array [0..1000] of Real;    //création type tableau
var a, b, n,i: Integer;                   // équation sous la forme ay'+by=g(t)
y0,g,h,u,t1,ti: real;                     //initialisation variables
TableauY,Tableaut,TableauG : Tableau;             //initialisation type tableau

begin
	a:=1;
	b:=3;
	y0:=0;
	
	readln(h);
	readln(n);
	
	Tableaut[0]:=0;
	TableauY[0]:=y0;
	
	
	For i:=1 to n do
	begin
		ti=Tableaut[i-1];
		Tableau[i]:=ti+h;
		t1:=Tableau[i];
		TableauG[i]:=sin(t1);
		g:=TableauG[i];
		u:=TableauY[i-1];
		TableauY[i]:=u+h*((g-b*u)/a);
		
		writeln(Tableaut[i]:12:10,'',TableauY[i]:12:10);
	end;
end.