Runge Kutta 4

Version imprimable

Runge Kutta 4

Bonjour!
Je cherche à faire fonctionner runge kutta 4 sous c++ pour modéliser la trajectoire d'une planète autour du soleil. J'ai beau rentrer correctement le code à base des équations de Newton, impossible de trouver des résultats convenables.. Mon x ne bouge pas de position tandis que y diverge complètement...
J'ai tout essayé, du coup j'appelle à l'aide maintenant. Voici mon code , merci à tous ceux qui prendront le temps de m'aider ! :)

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
double f(double x,double y){
    double G=6.674*pow(10,-11);
    double Ms=1.989*pow(10,30);
    double Aso=13.34*pow(10,19);
    double F;
    double r=sqrt(x*x+y*y);
    F=-((Aso*x)/(r*r*r));
    return F;
 
}    
double g(double x,double y){
    double G=6.674*pow(10,-11);
    double Ms=1.989*pow(10,30);
    double Aso=13.34*pow(10,19);
    double r=sqrt(x*x+y*y);
    double g;
    g=-((Aso*y)/(r*r*r));
    return g;
}    
 
 
 
double RK4_PLANETE(double x, double y, double u, double v, double n){
 
    double h=0.1; 
   double kx1,kx2,kx3,kx4,A,B; 
   double ky1,ky2,ky3,ky4; 
   double ku1,ku2,ku3,ku4; 
   double kv1,kv2,kv3,kv4;
   ofstream traj("traj.txt");
 
    for (float t=0 ; t<1 ; t+=h){
 
    kx1=u*h;
    ky1=v*h;
    ku1=f(x,y)*h;
    kv1=g(x,y)*h;
 
    kx2=(u+ku1/2.0)*h;
    ky2=(v+kv1/2.0)*h;
    ku2=(f(x+kx1/2.0,y+ky1/2.0))*h;
    kv2=(g(x+kx1/2.0,y+ky1/2.0))*h;
 
    kx3=(u+ku2/2.0)*h;
    ky3=(v+kv2/2.0)*h;
    ku3=(f(x+kx2/2.0,y+ky2/2))*h;
    kv3=(g(x+kx2/2.0,y+ky2/2))*h;
 
    kx4=(u+ku3)*h;
    ky4=(v+kv3)*h;
    ku4=(f(x+kx3,y+ky3))*h;
    ku4=(g(x+kx3,y+ky3))*h;
 
    x=x+((1/6.0)*(kx1+2*kx2+2*kx3+kx4));
    y=y+((1/6.0)*(ky1+ky2*2+ky3*2+ky4));
    u=u+((1/6.0)*(ku1+ku2*2+ku3*2+ku4));
    v=v+((1/6.0)*(kv1+kv2*2+kv3*2+kv4));
 
    traj<<x<<"  "<<y<<endl;     
 
    }; 
 
}

05/05/2018, 22h29
dourouc05

:salut:

Déjà, ta manière d'encoder tes constantes est assez particulière : tu pourrais écrire tout simplement 6.674e-11, par exemple. Aussi, au lieu de définir une variable g, d'y assigner une valeur, puis de la retourner, tu pourrais simplement retourner -((Aso*y)/(r*r*r)).

As-tu vérifié le résultat de f et g ? (D'ailleurs, que sont f et g, très exactement ?)
Ton pas h est-il assez petit ?
05/05/2018, 23h09
aureliaa

Bonjour, Oui je conçois mais c'est pour que je sache à peu près ce que je fais ! :)

f et g correspondent aux valeurs de la Force de gravitation pour résoudre l'équation de Newton .

Pour h j'ai tout essayé, ca ne change rien ..
05/05/2018, 23h16
aureliaa

1 pièce(s) jointe(s)

Pièce jointe 377080

Voici la figure que j'obtiens, quand je code avec un pas de1 pour 365 itérations..
07/05/2018, 08h30
wolinn

Je n'ai pas vérifié le détail des formules utilisées, juste qu'elles ont l'air bien symétriques en x, y, il ne parait pas y avoir de fautes de frappe à ce niveau.
Ca peut être un problème de précision numérique.
Par exemple, lorsque tu calcules kx2 = (u + ku1/2)*h, que vaut le rapport ku1/u ?
Si ce rapport est inférieur à 1e-15, le bruit numérique dû aux erreurs d'arrondis devient dominant, avec la précision du type 'double'.
Tu pourrais ajouter un test et afficher une alerte si cela se produit.
Si au contraire il est trop grand, les dérivées sont peut-être mal calculées.
Quelles sont les conditions de convergence de la méthode et comment s'accumulent les erreurs ?
Par ailleurs, est-ce que les coordonnées cartésiennes sont bien le meilleur choix et est-ce que ça ne marcherait pas mieux en polaire ?
(pour une trajectoire de planète, r est une fonction lentement variable, et theta une fonction à dérivée lentement variable, ce qui est plutôt intéressant pour une résolution numérique stable d'équation différentielle).

Ce sont plus des questions d'algorithmique que de langage de programmation...