Runge-Kutta à une variable?

Version imprimable

08/11/2006, 18h33
PadawanDuDelphi

Runge-Kutta à une variable?

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.

@+.
08/11/2006, 19h02
millie

Code:

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.
08/11/2006, 19h07
j.p.mignot

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

Citation:

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.

Citation:

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
09/11/2006, 09h51
PadawanDuDelphi

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.

Citation:

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:

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:

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

Encore merci...

@+.
09/11/2006, 10h50
millie

Ici, tu as :

Citation:

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)).
09/11/2006, 15h19
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...)
09/11/2006, 16h16
millie

Citation:

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. ;)
09/11/2006, 16h41
PadawanDuDelphi

Merci bocoup,

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

@+.
09/11/2006, 19h06
j.p.mignot

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.
10/11/2006, 11h09
PadawanDuDelphi

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.

@+.

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:

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

@+.

21/11/2006, 10h24
millie

Code:

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.

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

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.

21/11/2006, 13h02
millie

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

Citation:

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 ?
21/11/2006, 14h02
PadawanDuDelphi

Citation:

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.

@+.

Voilà pour me faire pardonner.

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

Je résout y' = ay.

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
#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:

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;
}

21/11/2006, 16h39
PadawanDuDelphi

Merci pour tout,

Ton code fonctionne nickel...

@+.
21/11/2006, 18h00
j.p.mignot

Ne pas oublier, si cela est nécessaire, les pas adaptatifs. Cela peut considérablement réduire les temps de calculs!
22/11/2006, 09h20
PadawanDuDelphi

Citation:

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 :D )

En tout cas merci pour tout.