Runge-Kutta à une variable?

**PadawanDuDelphi** · 08/11/2006, 18h33

Bonjour,

J'ai un problème avec Runge-Kutta...Je pensais avoir compris mais là ça coince (sur un truc simple, j'en suis sûr). J'ai une équation du type:

dx/dt = f(x) uniquement.
Or Runge-Kutta me demande de calculer des valeurs K1=deltat * f(x,y);

J'ai bien lu le numeric recipices et les autres liens données sur le forum, mais je coince toujours. Quelqu'un aurait-il un exemple simple pour moi?

Merci bocoup.

@+.

**millie** · 08/11/2006, 19h02

Code :

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

dx/dt = f(x) uniquement.

Ca me semble bizarre : x est ici une fonction (tu la dérives par rapport à t). Donc, il faudrait au moins que f(x) dépende de t. A mon avis, tu as une équations du type (comme Runge-Kutta) du type : dx/dt(t) = f(x,t)

Par exemple :

dx/dt(t) = cos(x(t)) avec f : (x,y) -> cos(x(t)) où x est une fonction.

**j.p.mignot** · 08/11/2006, 19h07

Voir par exemple
http://www.ac-nancy-metz.fr/enseign/...doc/RK4_02.pdf
bien entendu la connaissance de

dx/dt(t) = f(x,t)

est insuffisante.
la connaissance de X(t=0) est absolument necessaire.
la solution peut changer drastiquement avec cette initialisation.

Ca me semble bizarre

non ; il faut trouver X(t) tel que dX(t)/dt = f(X(t),t)
si dans un cas particulier on peut exprimer dX/dt = f(x(t)) on a le temps via x(t)
comme par exemple avec
dx/dt = x => d(log(x)) = dt => log(x) = t +K => x = K1.e^t

Note
J'avais déjà débatu sur ce sujet et sur ce forum: voir
http://www.developpez.net/forums/sho...d.php?t=104554

**PadawanDuDelphi** · 09/11/2006, 09h51

Tout d'abord merci pour vos réponse.

J'avais déjà regarder le topic, mais les équations différentielles sont loin derrière moi et j'ai un peu de mal.

Ca me semble bizarre : x est ici une fonction (tu la dérives par rapport à t). Donc, il faudrait au moins que f(x) dépende de t. A mon avis, tu as une équations du type (comme Runge-Kutta) du type : dx/dt(t) = f(x,t)

Par exemple :

dx/dt(t) = cos(x(t)) avec f : (x,y) -> cos(x(t)) où x est une fonction.

Ce dois être ça. En fait le plus simple c'est peut-être que je mette l'équation:

Code :

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

dO/dt=A*O(t)+B*O(t)^2+...+DO(t)^4

et on a O0 pour t=0.
Je vois pas avec ça comment trouver, par exemple

Code :

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

K2=deltat*f(O0+deltat/2,t0+K1/2)

Encore merci...

@+.

**millie** · 09/11/2006, 10h50

Ici, tu as :

dO/dt=A*O(t)+B*O(t)^2+...+DO(t)^4

Donc : f(x,t) = a * x(t) + b*x(t)²+...d*x(t)^4

Donc :

K2=deltat*f(O0+deltat/2,t0+K1/2) = deltat * (a * (0 + deltat/2) + b * (....))

Car f(x,t) t n'intervient pas seul dans son expression (c'est toujours x(t)).

**PadawanDuDelphi** · 09/11/2006, 15h19

Ok,

Alors si j'ai bien compris mon calcul de K1, K2 et K3 ne seront pas récursifs? (c.a.d que K3 ne dépendra pas de K2 qui ne dépendra pas de K1...)

**millie** · 09/11/2006, 16h16

Envoyé par PadawanDuDelphi

Ok,

Alors si j'ai bien compris mon calcul de K1, K2 et K3 ne seront pas récursifs? (c.a.d que K3 ne dépendra pas de K2 qui ne dépendra pas de K1...)

Lors d'une étape, effectivement, les Ki ne dépendront pas l'un de l'autre vu que t n'apparait pas explicitement seul dans l'équation diff que tu as.

**PadawanDuDelphi** · 09/11/2006, 16h41

Merci bocoup,

Je pense avoir compris. Je vais implémenter ça et je vous tiens au courant.

@+.

**j.p.mignot** · 09/11/2006, 19h06

On trouve déjà de noubreux exemples de code sur le WEB
comme
http://www.physics.orst.edu/~rubin/n...k_c/srk_c.html
ou
http://www.physics.orst.edu/~rubin/n...k_c/srk_f.html
ou
http://www.library.cornell.edu/nr/bookcpdf/c16-1.pdf
ou
http://beige.ucs.indiana.edu/B673/node53.html

...
Si vous avez quelques problèmes, je pourrais vous fournir un extrait de mes codes perso fait pour delphi ( Pascal ) ou C: je ne sais plus!

La méthode est facile à implémenter ( finalement 4 X Euler ). Cela devient + délicat avec des pas adaptatifs.

**PadawanDuDelphi** · 10/11/2006, 11h09

Merci bocoup,

J'ai vu qu'il y avait des exemples en Fortran, et c'est exactement le langage que je souhaite utilisé.

Encore merci à tous les deux pour vos exemples, vos conseils et votre patience.

@+.

**PadawanDuDelphi** · 21/11/2006, 10h12

Bonjour,

Je relance le sujet car il semblerair que ma méthode d'implémentation ne fonctionne pas. Mon problème c'est lorsque j'utilise cette algorithme pour résoudre des équations simples, et que je les compare avec une méthode de résolution manuelle, alors les résultats sont différents dans la partie transitoire.
Je place donc mon algorithme dans l'espoir que quelqu'un puisse y trouver une erreur.
Merci d'avance.

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
f = dX/dt = K1*X(t) + K2 + K3*X(t)^4  // version simplifiee de l'equation
C.I : X(t=0) = 20.
/******************************/
pas=0
deltaT = 1
 
Xinitial =  f(20)
 
For pas=0 to Fin_pas do:
{
K1 = f(Xinitial) * deltaT;
K2 = f(Xinitial + K1/2) * deltaT
K3 = f(Xinitial + K2/2) * deltaT
K4 = f(Xinitial + K3) * deltaT
 
deltaX = (1/6) * [K1 + 2K2 + 2K3 +K4]
Xinitial = Xinitial + deltaX
afficher (Xinitial)
pas++

Encore merci

@+.

**millie** · 21/11/2006, 10h24

Code :

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

f = dX/dt = K1*X(t) + K2 + K3*X(t)^4

Je ne comprends pas ce que tes Ki viennent faire ici. Ils ne devraient pas être dans la définition de ta fonction.

Note que f:t -> aX(t) + b X(t)^4 + c

EDIT : Attention, note que sur leur page, leur x, c'est ton t, et ton x c'est leur y.

**PadawanDuDelphi** · 21/11/2006, 12h05

En fait je m'était inspiré de cet exemple, qui semblait correspondre au mien, pour l'implémentation:

Code :

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

1
2
3
4
5
6
7
Exercice RK4_1 : décharge d’un condensateur dans une résistance.

du/dt =  - (1/RC) * U(t)
 
K1 =  - (1/RC * Ui) dt
K2 = -(1/RC * [Ui+K1/2])dt
.....

Merci.

**millie** · 21/11/2006, 13h02

Arf, désolé pour le post précédent, j'ai craqué, c'est moi qu'ait tout confondu

mon problème c'est lorsque j'utilise cette algorithme pour résoudre des équations simples"]on problème c'est lorsque j'utilise cette algorithme pour résoudre des équations simples

Quel genre d'équation tu as cherché à résoudre. Comment appelles-tu la fonction f ?

**PadawanDuDelphi** · 21/11/2006, 14h02

Envoyé par millie

Arf, désolé pour le post précédent, j'ai craqué, c'est moi qu'ait tout confondu

C'est pas grave, ce qui compte c'est sympa à toi de t'intérésser à mon problème.

En fait je cherche à résoudre une equation thermique; sa dérivée est simplifiable en un polynôme de degré 4...
Et pour moi f c'est la derivee de ma fonction (=dx/dt).

Sinon j'avais essayer de résoudre dx/dt = A * X(t), avec X(t=0) = 1.
La solution devrait donc être X(t) = exp(At).

Mais je retrouve pas ma courbe exponentielle exacte (en fait, ma courbe calculée avec Runge-Kutta atteint la valeur limite bien plus vite que la fonction exp(At).

Encore merci.

@+.

**millie** · 21/11/2006, 14h38

Voilà pour me faire pardonner.

Tiens, voici un code que j'ai écris et testé et ça marche.

Je résout y' = ay.

Code cpp :

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
#include <iostream>
#include <cmath>
 
 
void rungeK(double dt, double y0,  double (*f)(double), double (*solution) (double))
{
  double k1;
  double k2;
  double k3;
  double k4;
 
  double y = y0;
 
  for(int i = 0; i<20;i++)
  {
 
      k1 = f(y) * dt;
      k2 = f(y+k1/2.0) *dt;
      k3 = f(y+k2/2.0) *dt;
      k4 = f(y+ k3) *dt;
 
      y += 1.0/6.0 * (k1 + 2.0*k2 + 2.0*k3 + k4);
      std::cout<<"Y : "<<y<<" Vrai solution: "<<solution(x0 + dt *(i+1))<<std::endl;
  }
 
}
 
double expA(double x)
{
 double a = 2.0;
 return exp(a* x);
}
 
double IdentiteA(double x)
{
  double a= 2.0;
  return a * x;
}
 
int main(void)
{
  double x0 = 0.0;
  rungeK(x0, 0.001, expA(x0), IdentiteA, expA);
 
  return EXIT_SUCCESS;
}

Et donc pour toi, cela devrait être de la forme :

Code cpp :

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
#include <iostream>
#include <cmath>
 
 
void rungeK(double x0, double dt, double y0,  double (*f)(double))
{
  double k1;
  double k2;
  double k3;
  double k4;
 
  double y = y0;
 
  for(int i = 0; i<20;i++)
  {
 
      k1 = f(y) * dt;
      k2 = f(y+k1/2.0) *dt;
      k3 = f(y+k2/2.0) *dt;
      k4 = f(y+ k3) *dt;
 
      y += 1.0/6.0 * (k1 + 2.0*k2 + 2.0*k3 + k4);
      //traiter y
  }
 
}
 
double eqnFun(double x)
{
  double a= 2.0;
  double b = 1.0f;
  double c = 1.0f;
 
  return(a * pow(x,4.0) + b *x + c);
}
 
int main(void)
{
  double y0 = 0.0;
  rungeK(0.001, y0, eqnFun);
 
  return EXIT_SUCCESS;
}

**PadawanDuDelphi** · 21/11/2006, 16h39

Merci pour tout,

Ton code fonctionne nickel...

@+.

**j.p.mignot** · 21/11/2006, 18h00

Ne pas oublier, si cela est nécessaire, les pas adaptatifs. Cela peut considérablement réduire les temps de calculs!

**PadawanDuDelphi** · 22/11/2006, 09h20

Ne pas oublier, si cela est nécessaire, les pas adaptatifs. Cela peut considérablement réduire les temps de calculs!

Effectivement, j'avais lu ça dans un article, mais ça m'a paru un peu compliqué pour mon niveau (surtout que ce calcul est utilisé dans un modèle avec deux autre calculs qui utilisent un très grand nombre d'inversion de matrices 150*150...Donc la vitesse n'ai pas un problème apparemment

)

En tout cas merci pour tout.

Runge-Kutta à une variable?

Algorithmes et structures de données

Discussions similaires

Partager

Partager