Forum |  HardWare.fr | News | Articles | PC | S'identifier | S'inscrire | Shop Recherche
2179 connectés 

  FORUM HardWare.fr
  Programmation
  C

  [C] mouvement d'un double pendule

 
 


 Mot :   Pseudo :  
 
Bas de page
Auteur Sujet :

[C] mouvement d'un double pendule

n°1519640
pedigus1
Posté le 25-02-2007 à 15:32:26  profilanswer
 

Salut à tous,
j'aimerais representer à l'écran (avec SDL) le mouvement d'un double pendule attaché a un plafond par une barre indeformable.Les equations differentielles du mouvement sont là: http://fr.wikipedia.org/wiki/Double_pendule et je pense qu'utiliser la methode d'integration de runge kutta d'ordre 4 permettrai d'avoir quelquechose de stable. Mais j'arrive pas a traduire cette methode en code. Si vous pouviez m'expliquer clairement le principe ce serait cool (même si ce n'est que pour un pendule simple). J'ai à peine 6 mois de C derriere moi.
Merci

mood
Publicité
Posté le 25-02-2007 à 15:32:26  profilanswer
 

n°1519659
Sve@r
Posté le 25-02-2007 à 16:33:42  profilanswer
 

http://www.cpt.univ-mrs.fr/~grapperon/oldpdf/mip.pdf


---------------
Vous ne pouvez pas apporter la prospérité au pauvre en la retirant au riche.
n°1519691
GrosBocdel
Posté le 25-02-2007 à 18:22:48  profilanswer
 

Pour utiliser ta méthode de Runge Kutta, il faut dans un premier temps que tu te ramènes à sa situation d'application, c'est à dire dx/dt=f(x)
f est ici très fortement non linéaire, sinon ce ne serait pas drôle.
Tu te prends un double x[4];
x[0]=theta1
x[1]=theta2
x[2]=dtheta1/dt
x[3]=dtheta2/dt
 
à t=0, tu sais le remplir, ce sont tes conditions initiales.
 
Tu te prends un double y[4], c'est ton dx/dt.
y[0]=dtheta1/dt
y[1]=dtheta2/dt
y[2]=d2theta1/dt2
y[3]=d2theta2/dt2
 
Pour y[0] et y[1], pas trop difficile, y[0]=x[2] et y[1]=x[3]. (pas certain qu'on s'en serve...)
Pour y[2] et y[3], tu mets tes équations de mouvement en fonction de x[0], x[1], etc...
 
Ca, c'était pour t=0. Tu te choisis un pas de calcul h ad hoc. S'il est trop grand, le calcul part dans le décors.
Tu te prends un double k[4] et tu calcules les k1,k2,k3,k4 de la méthode de Runge Kutta. http://en.wikipedia.org/wiki/Runge-Kutta. Ici les équations de mouvement ne dépendent pas du temps, donc sur la page wiki f(t,y) tu zappes et c'est f(y). Avec ça, tu obtiens tes nouvelles positions et vitesses à t0+h, que tu réinjectes dans x[] et tu boucles.


Message édité par GrosBocdel le 25-02-2007 à 18:25:51
n°1519758
pedigus1
Posté le 25-02-2007 à 23:21:45  profilanswer
 

merci de vos réponses.
Donc d'après ce que j'ai compris, il faut que je ramene d'abord mes equa diff a une forme:    d²théta1/dt²=f(théta2,théta2) et d²théta2/dt²=f(théta2,théta1)      
je vais ensuite me servir de cette fonction f pour calculer mes k1,k2...  
Je vais obtenir mes nouveaux thétas, mais est ce que ma vitesse sera la dérivée de mon théta? En plus j'ai du cos et sin dans mes equa diff, il existe pas des bibliotheques qui ont deja ces fonctions?
Je crois avoir compris le principe pour une equa diff du 1e degré à une seule inconnue mais la avec ces deux equa diff du second degré liée c'est plus dur.

Message cité 1 fois
Message édité par pedigus1 le 25-02-2007 à 23:23:21
n°1520206
pedigus1
Posté le 26-02-2007 à 15:38:43  profilanswer
 

Pour y[2] et y[3], je dois mettre mes équations de mouvement en fonction de x[0], x[1], etc... c'est la que je bloque.Est-ce que je dois isoler dans chaque equation mon terme qui est dérivé 2 fois? Et comment est ce que je peut calculer mes K1,K2,K3,K4 de runge kutta avec deux equa diff en même temps?

Message cité 1 fois
n°1520342
GrosBocdel
Posté le 26-02-2007 à 17:36:12  profilanswer
 

pedigus1 a écrit :


Je vais obtenir mes nouveaux thétas, mais est ce que ma vitesse sera la dérivée de mon théta? En plus j'ai du cos et sin dans mes equa diff, il existe pas des bibliotheques qui ont deja ces fonctions?


 
Oui et non. Ce qu'on veut ici c'est intégrer les équations du mouvement, pour ça, on prend un pas de calcul h, et on suppose un truc du style :
d2theta/dt2 * h = dtheta/dt
donc v(t+h)=v(t)+a(t)*h. Les k1, k2 etc calculent juste une pente.
sin et cos sont dans math.h.  L'intérêt de cette méthode c'est de ne pas avoir à linéariser les sinus, en ayant une solution exacte (même si elle n'est pas analytique)
 

pedigus1 a écrit :


Pour y[2] et y[3], je dois mettre mes équations de mouvement en fonction de x[0], x[1], etc... c'est la que je bloque.Est-ce que je dois isoler dans chaque equation mon terme qui est dérivé 2 fois?


Oui, tu mets tes theta'' à gauche, le reste à droite. Faut juste recopier les équations que tu as, en remplaçant theta1' par x[2], theta2' par x[3], etc...
 

pedigus1 a écrit :


 Et comment est ce que je peut calculer mes K1,K2,K3,K4 de runge kutta avec deux equa diff en même temps?


Arrivé à ce stade du calcul, tu n'as plus 2 équa diff du second ordre, mais une seule du premier ordre, dx/dt=f(x), avec x qui est un vecteur à 4 valeurs. Tu as dû voir dans tes cours un truc qui s'appelle l'espace des phases.
Je me suis mal exprimé sur cette partie : Les K1, K2, K3, K4 sont chacun des double[4]. K1 te donne une estimation de la variation de chacun de tes paramètres, K2 une autre un peu plus loin, etc... La somme pondérée de chacun des K te donne une estimation plus précise de la variation de chacun des paras.
Sinon, tu peux faire une recherche anglophone sur le pendule double, il y a un code en c qui traine. Comme ce que tu fais a l'air d'un exo, je ne te le donne pas  ;)  

n°1520532
pedigus1
Posté le 26-02-2007 à 20:29:31  profilanswer
 

Déjà merci de m'aider c'est très sympa
En fait c'est pour un projet info pour mon école.J'ai trouvé le code in english mais bon je comprend pas tout et j'ai pas vraiment envi de recopier.
Voici ce que j'ai commencé à coder:  

Code :
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <math.h>
  5. int main()
  6. {
  7. double x[4],y[4],K1[4],K2[4],K3[4],K4[4];
  8. int m1,m2,l1,l2;  /*m1 et m2 masses des pendules, l1 et l2 longeurs des fils*/
  9. float h,g=9.81;            /*le pas et la constante de gravitation universelle*/
  11. /* à t=0 */
  12. y[0]=x[2];
  13. y[1]=x[3];
  14. y[2]=(1/((m1+m2)l1))*(-m2*l2*y[3]*cos(x[0]-x[1])-m2*l2*x[3]*x[3]*sin(x[0]-x[1])-(m1+m2)*g*sin(x[0]));
  15. y[3]=(1/l2)(-l1*y[2]*cos(x[0]-x[1])+l1*x[2]*sin(x[0]-x[1])-g*sin(x[1]));
  18. /*Runge-kutta*/
  19. K1[0]=y[0];
  20. K1[1]=y[1];
  21. K1[2]=Y[2];
  22. K1[3]=y[3];
  24. K2[0]=y[0]+(h/2)*K1[0];
  25. K2[1]=y[1]+(h/2)*K1[1];
  26. K2[2]=y[2]+(h/2)*K1[2];
  27. K2[3]=y[3]+(h/2)*K1[3];


 
En fait mon equa serait devenue y=f(x) (avec y et x mes doubles)
Mais est-ce normal que je me retrouve avec du y[3] dans mon y[2]? et j'arrive pas bien à identifier mon f pour le y[2] et y[3].De même pour le calcul des K2 sa devrait être f(y[]+(h/2)K1[]) ,pareil, quel est mon f?

Message cité 1 fois
Message édité par pedigus1 le 28-02-2007 à 21:37:41
n°1520602
GrosBocdel
Posté le 26-02-2007 à 22:47:08  profilanswer
 

pedigus1 a écrit :


 y[0]=x[2];
 y[1]=x[3];
 y[2]=(1/((m1+m2)l1))*(-m2*l2*y[3]*cos(x[0]-x[1])-m2*l2*x[3]*x[3]*sin(x[0]-x[1])-(m1+m2)*g*sin(x[0]));
 y[3]=(1/l2)(-l1*y[2]*cos(x[0]-x[1])+l1*x[2]*sin(x[0]-x[1])-g*sin(x[1]));
 


Faut penser à initialiser x avant et à ne pas oublier les * des multiplications.
 

pedigus1 a écrit :


Mais est-ce normal que je me retrouve avec du y[3] dans mon y[2]? et j'arrive pas bien à identifier mon f pour le y[2] et y[3].De même pour le calcul des K2 sa devrait être f(y[]+(h/2)K1[]) ,pareil, quel est mon f?


Non, ce n'est pas normal. Il y a encore un peu de travail à effectuer sur les équations 1 et 2 de la page wikipedia. Deux équations, deux inconnues (theta1'' et theta2'') à isoler. Sinon la fonction f est la même partout dans le programme. C'est celle qui transforme ton x[] en ton y[]. Alors autant la coder séparément.


Message édité par GrosBocdel le 26-02-2007 à 22:51:10
n°1520970
pedigus1
Posté le 27-02-2007 à 14:10:40  profilanswer
 

J'ai corriger les * des multiplications et les x je verrais sa plus tard.  
En fait il faut que j'arrive a ramener mon système de deux équations a 2 inconnues liées à un système de 4 equations du 1e ordre.c'est sa? Et pour la fonction f, je n'en aurais pas quatres differentes? (une pour y[0] une autre pour y[1] etc..).
Et sinon pour le travail a effectuer sur les equations 1 et 2 de wikipedia, est-ce que je dois remplacer mes y[2] de l'equation y[3] par le y[2] trouvé juste avant?
Merci de ton aide , je suis assez perdu là-dessus.

n°1521256
GrosBocdel
Posté le 27-02-2007 à 20:10:29  profilanswer
 

Les équations de wikipedia, dans un premier temps, il faut juste les mettre en forme. Elles n'ont juste pas la tête qu'il faut. La résolution des équa diff, c'est bien Runge Kutta qui s'en occupe.
 
aX+bY=c
dX+eY=d
Je pense que tu sais tirer X et Y d'un système d'équations comme ça. Là c'est pareil, sauf qu'il faut le faire avec theta1'' et theta2''. Juste pour mettre en forme comme il faut.
Une fois que tu as les expressions de theta1'' (sans theta2'' dedans mais avec tout ce que tu veux de theta1' theta2', theta1,theta2) et theta2''(sans theta1'' dedans) c'est tout c'est fini!  Tu remplis y[2] et y[3] avec ça, comme tu avais commencé.
Oui, si tu veux, tu as 4 f différentes. Tu veux voir ça comme ça. T'étais bien parti, ton prog a la bonne forme. Tes y[0] et y[1] sont bons, tu as juste à changer y[2] et y[3].

mood
Publicité
Posté le 27-02-2007 à 20:10:29  profilanswer
 

n°1521260
pedigus1
Posté le 27-02-2007 à 20:20:40  profilanswer
 

Voila l'ébauche du code:

Code :
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <math.h>
  5. int main()
  6. {
  7. double x[4],y[4],x1[4],y1[4];
  8. double K1[4],K2[4],K3[4],K4[4];
  9. int m1,m2,l1,l2; 
  10. float h,g=9.81;
  12. /* à t=0 */
  13. y[0]=x[2];
  14. y[1]=x[3];
  15. y[2]=-(1/((m1+m2*sin(x[0]-x[1])*sin(x[0]-x[1]))*l1))*(m2*(-l1*x[2]*sin(x[0]-x[1])+g*sin(x[1]))*cos(x[0]-x[1])+m2*l2*x[3]*x[3]*sin(x[0]-x[1])+(m1+m2)*g*sin(x[0]));
  16. y[3]=-(1/l2)*(l1*(-(1/((m1+m2*sin(x[0]-x[1])*sin(x[0]-x[1]))*l1))*(m2*(-l1*x[2]*sin(x[0]-x[1])+g*sin(x[1]))*cos(x[0]-x[1])+m2*l2*x[3]*x[3]*sin(x[0]-x[1])+(m1+m2)*g*sin(x[0])))*cos(x[0]-x[1])-l1*x[2]*sin(x[0]-x[1])+g*sin(x[1]));
  20. /*Runge-kutta*/
  21.     /*pour eq0*/
  22. K1[0]=fonction0(x[0]);
  23. K2[0]=fonction0(x[0]+(h/2)*K1[0]);
  24. K3[0]=fonction0(x[0]+(h/2)*K2[0]);
  25. K4[0]=fonction0(x[0]+h*K2[0]);
  27. x1[0]=x[0]+(h/6)(K1[0]+2*K2[0]+2*K3[0]+K4[0]);
  28.  
  29.    /*pour eq1*/
  30. K1[1]=fonction1(x[1]);
  31. K2[1]=fonction1(x[1]+(h/2)*K1[1]);
  32. K3[1]=fonction1(x[1]+(h/2)*K2[1]);
  33. K4[1]=fonction1(x[1]+h*K2[1]);
  35. x1[1]=x[1]+(h/6)(K1[1]+2*K2[1]+2*K3[1]+K4[1]);
  37.    /*pour eq2*/
  38. K1[2]=fonctionf(x[2]);
  39. K2[2]=fonctionf(x[2]+(h/2)*K1[2]);
  40. K3[2]=fonctionf(x[2]+(h/2)*K2[2]);
  41. K4[2]=fonctionf(x[2]+h*K2[2]);
  43. x1[2]=x[2]+(h/6)(K1[2]+2*K2[2]+2*K3[2]+K4[2]);
  45.    /*pour eq3*/     
  46. K1[3]=fonctiong(x[3]);
  47. K2[3]=fonctiong(x[3]+(h/2)*K1[3]);
  48. K3[3]=fonctiong(x[3]+(h/2)*K2[3]);
  49. K4[3]=fonctiong(x[3]+h*K2[3]);
  51. x1[3]=x[3]+(h/6)(K1[3]+2*K2[3]+2*K3[3]+K4[3]);
  52.  
  53.    /*injection des nouvelles valeurs dans les anciennes*/
  54. x[0]=x1[0];
  55. x[1]=x1[1];
  56. x[2]=x1[2];
  57. x[3]=x1[3];


 
 
Avec ces fonctions (qui d'ailleur devraient me permettre de supprimer mon paragraphe de /*à t=0*/ non?):

Code :
  1. double fonction0(double t[0])
  2. {
  3. double f[0];
  4. f[0]=t[0];
  5. return f[0];
  6. }
  8. double fonction1(double s[1])
  9. {
  10. double j[1];
  11. j[1]=s[1];
  12. return j[1];
  13. }
  15. double fonctionf(double v[2])
  16. {
  17. double F[2];
  18. F[2]=-(1/((m1+m2*sin(x[0]-x[1])*sin(x[0]-x[1]))*l1))*(m2*(-l1*v[2]*sin(x[0]-x[1])+g*sin(x[1]))*cos(x[0]-x[1])+m2*l2*x[3]*x[3]*sin(x[0]-x[1])+(m1+m2)*g*sin(x[0]));
  19. return F[2];
  20. }
  22. double fonctiong(double w[3])
  23. {
  24. double G[3];
  25. G[3]=-(1/l2)*(l1*(-(1/((m1+m2*sin(x[0]-x[1])*sin(x[0]-x[1]))*l1))*(m2*(-l1*x[2]*sin(x[0]-x[1])+g*sin(x[1]))*cos(x[0]-x[1])+m2*l2*w[3]*w[3]*sin(x[0]-x[1])+(m1+m2)*g*sin(x[0])))*cos(x[0]-x[1])-l1*x[2]*sin(x[0]-x[1])+g*sin(x[1]));
  26. return G[3];
  27. }


Message édité par pedigus1 le 28-02-2007 à 21:37:17
n°1521269
GrosBocdel
Posté le 27-02-2007 à 20:43:15  profilanswer
 

Bon c'est pas super optimisé, mais si tu ne t'es pas gouré dans les équations (trouves en des toutes faites sur le net, là les tiennes sont simplifiables), si tu mets tout ça dans une boucle for et que tu initialises tes x avant, ça doit pouvoir te sortir des résultats. Peut-être.
Oui tu peux virer ton paragraphe à t=0


Message édité par GrosBocdel le 27-02-2007 à 20:56:00
n°1521286
pedigus1
Posté le 27-02-2007 à 21:31:52  profilanswer
 

Comment est-ce que je pourrais optimiser plus?Et est-ce que mes fonctions sont correctes?J'ai mis des tableux en argument, je devrai pas rajouter une étoile devant?est est-ce que je ne pourrais pas remplacer ces tableux par juste des doubles?
Sinon j'ai chercher sur le net,et je n'ai pas trouver q'autres equations.
Merci de ton aide

Message cité 1 fois
n°1521351
GrosBocdel
Posté le 28-02-2007 à 06:58:18  profilanswer
 

pedigus1 a écrit :

Comment est-ce que je pourrais optimiser plus?Et est-ce que mes fonctions sont correctes?J'ai mis des tableux en argument, je devrai pas rajouter une étoile devant?est est-ce que je ne pourrais pas remplacer ces tableux par juste des doubles?
Sinon j'ai chercher sur le net,et je n'ai pas trouver q'autres equations.
Merci de ton aide


 
Bon, j'ai parlé trop vite. Des fois avec des trucs comme ça, sur un heureux hasard, ça peut marcher mais là, non.
Dans un premier temps, ouvre vite un bouquin de c. Là tu n'y coupes pas.
Après, je t'ai dis qu'une fonction f, il n'y en a qu'une. Donc à toi de coder une void f(double x[], double y[]), qui te prend ton x et te calcule y.
Mais là, il faut que tu ouvres un bouquin de c. C'est obligatoire.
 
Pour te donner une idée : actuellement, ta fonction0 prend un double en paramètre, ta fonction1 en prend deux, fonctionf en prend trois, fonctiong en prend quatre ! Et aucune ne peut savoir qui est x avec lequel elles travaillent.
Donc un bouquin de c et après on en rediscute si tu rebloques sur des trucs pas expliqués dans le bouquin.
Sinon c'est un début, hein, pas de problème.


Message édité par GrosBocdel le 28-02-2007 à 07:08:11
n°1521672
pedigus1
Posté le 28-02-2007 à 14:08:07  profilanswer
 

Je n'ai pas beaucoup de C derriere moi mais j'y travaille.
Mais je ne vois pas comment je pourrais avoir qu'une seule fonction f , j'ai en tout quatre equa diff qui sont differentes et il faut que j'applique runge-kutta a chacune d'elles.Pour moi sa fait quatre fonctions.
J'ai changé mes paramètres de fonctions pour des doubles simples (et non plus des tableaux comme j'avais mis). Est- ce que sa merde juste pour les fonctions? ou y a t il de grosses erreurs ailleurs?

Message cité 2 fois
Message édité par pedigus1 le 28-02-2007 à 14:11:25
n°1521905
pedigus1
Posté le 28-02-2007 à 18:14:24  profilanswer
 

Bon j'ai trvaillé sur le code (avec l'aide du site du zero) et je suis arrivé à un truc qui marche(compile et execute sans problèmes). Mais comment est-ce que je pourrais verifier mes resultats? j'imprime mes théta dans un fchier txt.
Voila le code:
 

Code :
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <math.h>
  4. #include "entete.h"
  5. #include <SDL/SDL.h>
  7. double x[4];
  8. int m1,m2,l1,l2;
  10. int main(int argc, char *argv[])
  11. {
  12. double y[4],x1[4],y1[4];
  13. double K1[4],K2[4],K3[4],K4[4]; 
  14. float h;
  15. FILE *fi;
  16. int duree,tpsfin=0,tpsactuel=0;
  18. /*conditions intiales*/
  19. x[0]=10;
  20. x[1]=20;
  21. x[2]=30;
  22. x[3]=10;
  23. y[0]=10;
  24. y[1]=20;
  25. y[2]=30;
  26. y[3]=10;
  27. duree=100;
  28. m1=20;
  29. m2=10;
  30. l1=3;
  31. l2=6;
  32. h=0.05;
  34. fi=fopen("donnee.txt","w" );
  35. tpsfin = SDL_GetTicks()+duree;
  37. while(tpsactuel<tpsfin)
  38. {
  39. /*Runge-kutta*/
  40. fprintf(fi,"%f %f\n",x[0],x[1]);
  41.    /*pour eq0*/     
  42. K1[0]=fonction0(x[0]);
  43. K2[0]=fonction0(x[0]+(h/2)*K1[0]);
  44. K3[0]=fonction0(x[0]+(h/2)*K2[0]);
  45. K4[0]=fonction0(x[0]+h*K2[0]);
  47. x1[0]=x[0]+(h/6)*(K1[0]+2*K2[0]+2*K3[0]+K4[0]);
  48.  
  49.   /*pour eq0*/     
  50. K1[1]=fonction1(x[1]);
  51. K2[1]=fonction1(x[1]+(h/2)*K1[1]);
  52. K3[1]=fonction1(x[1]+(h/2)*K2[1]);
  53. K4[1]=fonction1(x[1]+h*K2[1]);
  55. x1[1]=x[1]+(h/6)*(K1[1]+2*K2[1]+2*K3[1]+K4[1]);
  57.    /*pour eq2*/     
  58. K1[2]=fonctionf(x[2]);
  59. K2[2]=fonctionf(x[2]+(h/2)*K1[2]);
  60. K3[2]=fonctionf(x[2]+(h/2)*K2[2]);
  61. K4[2]=fonctionf(x[2]+h*K2[2]);
  63. x1[2]=x[2]+(h/6)*(K1[2]+2*K2[2]+2*K3[2]+K4[2]);
  65.    /*pour eq3*/     
  66. K1[3]=fonctiong(x[3]);
  67. K2[3]=fonctiong(x[3]+(h/2)*K1[3]);
  68. K3[3]=fonctiong(x[3]+(h/2)*K2[3]);
  69. K4[3]=fonctiong(x[3]+h*K2[3]);
  71. x1[3]=x[3]+(h/6)*(K1[3]+2*K2[3]+2*K3[3]+K4[3]);
  72.  
  73.    /*injection des nouvelles valeurs dans les anciennes*/
  74. x[0]=x1[0];
  75. x[1]=x1[1];
  76. x[2]=x1[2];
  77. x[3]=x1[3];
  79. tpsactuel=SDL_GetTicks();
  80. }
  82. system("PAUSE" ); 
  83. return 0;
  84. }


 
Avec le fichier des fonctions:

Code :
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <math.h>
  5. double x[4];
  6. int m1,m2,l1,l2;
  7. float g=9.81;
  9. double fonction0(double t)
  10. {
  11. double f;
  12. f=t;
  13. return f;
  14. }
  16. double fonction1(double s)
  17. {
  18. double j;
  19. j=s;
  20. return j;
  21. }
  23. double fonctionf(double v)
  24. {
  25. double F;
  26. F=-(1/((m1+m2*sin(x[0]-x[1])*sin(x[0]-x[1]))*l1))*(m2*(-l1*v*sin(x[0]-x[1])+g*sin(x[1]))*cos(x[0]-x[1])+m2*l2*x[3]*x[3]*sin(x[0]-x[1])+(m1+m2)*g*sin(x[0]));
  27. return F;
  28. }
  30. double fonctiong(double w)
  31. {
  32. double G;
  33. G=-(1/l2)*(l1*(-(1/((m1+m2*sin(x[0]-x[1])*sin(x[0]-x[1]))*l1))*(m2*(-l1*x[2]*sin(x[0]-x[1])+g*sin(x[1]))*cos(x[0]-x[1])+m2*l2*w*w*sin(x[0]-x[1])+(m1+m2)*g*sin(x[0])))*cos(x[0]-x[1])-l1*x[2]*sin(x[0]-x[1])+g*sin(x[1]));
  34. return G;
  35. }


 
En tous cas, merci beaucoup pour ton aide !!!

n°1521908
GrosBocdel
Posté le 28-02-2007 à 18:24:46  profilanswer
 

pedigus1 a écrit :

Je n'ai pas beaucoup de C derriere moi mais j'y travaille.
Mais je ne vois pas comment je pourrais avoir qu'une seule fonction f , j'ai en tout quatre equa diff qui sont differentes et il faut que j'applique runge-kutta a chacune d'elles.Pour moi sa fait quatre fonctions.
J'ai changé mes paramètres de fonctions pour des doubles simples (et non plus des tableaux comme j'avais mis). Est- ce que sa merde juste pour les fonctions? ou y a t il de grosses erreurs ailleurs?


 
Tu définis juste ta fonction  

Code :
  1. void f(double x[], double y[])
  2. {
  3. y[0]=x[2];
  4. y[1]=x[3];
  5. y[2]=-(1/((m1+m2*sin(x[0]-x[1])*sin(x[0]-x[1]))*l1))*(m2*(-l1*x[2]*sin(x[0]-x[1])+g*sin(x[1]))*cos(x[0]-x[1])+m2*l2*x[3]*x[3]*sin(x[0]-x[1])+(m1+m2)*g*sin(x[0]));
  6. y[3]=-(1/l2)*(l1*(-(1/((m1+m2*sin(x[0]-x[1])*sin(x[0]-x[1]))*l1))*(m2*(-l1*x[2]*sin(x[0]-x[1])+g*sin(x[1]))*cos(x[0]-x[1])+m2*l2*x[3]*x[3]*sin(x[0]-x[1])+(m1+m2)*g*sin(x[0])))*cos(x[0]-x[1])-l1*x[2]*sin(x[0]-x[1])+g*sin(x[1]));
  7. }


 
Ca ne fait bien qu'une fonction. Le faire de cette manière va t'obliger à réécrire ton prog un peu différemment. Note que je ne vérifie pas la liste des sinus et cos qu'il y a dans tes expressions. Il y a peut-être une erreur, je ne sais pas.

Message cité 1 fois
n°1521911
red factio​n
Posté le 28-02-2007 à 18:28:25  profilanswer
 

mon dieu que cest compliqué !
tout ca pour trouver deux angles x[0] et x[1]!???

Message cité 1 fois
n°1521912
pedigus1
Posté le 28-02-2007 à 18:30:01  profilanswer
 

Effectivement c'est plus pratique, mais j'avais mal interpreté ce que tu voulais dire par fonction. Je pensais que tu parlais des fonctions mathematiques ( par exemple d²théta1/dt²=f(théta1,etc..)   et théta2/dt²=g(théta2,etc..) et f et g étaient differentes).
Mais je pense qu'il y quelquechose dans la méthode de runge kutta que je n'avais pas compris (je me fait aussi aider sur un forum de maths http://les-mathematiques.u-strasbg [...] ?15,357592 ).Donc je vais devoir changer pas mal de choses ce soir.


Message édité par pedigus1 le 28-02-2007 à 18:39:20
n°1521916
GrosBocdel
Posté le 28-02-2007 à 18:38:13  profilanswer
 

red faction a écrit :

mon dieu que cest compliqué !
tout ca pour trouver deux angles x[0] et x[1]!???


 
Tu as mieux à proposer? parce que c'est quand même un grand moment de solitude pour moi depuis le début, là !  :)

Message cité 1 fois
n°1521917
pedigus1
Posté le 28-02-2007 à 18:40:49  profilanswer
 

GrosBocdel a écrit :

Tu as mieux à proposer? parce que c'est quand même un grand moment de solitude pour moi depuis le début, là !  :)


 
Je suis d'accord (mais je m'en plains pas :) )

n°1521922
red factio​n
Posté le 28-02-2007 à 18:50:21  profilanswer
 

Non je n'ai rien de mieux , je pensait juste que cetait bcp plus simple a realiser
 
Je ne sais pas si tu doit faire ca pour lecole ou autre, mais peut etre que tu pourrait obtenir des resultats semblables avec un 3 particules (dont une fixe) connectées par des strings (tres tendues et avec un pas tres fin)
 
mais au final ca doit surement consommer plus de cpu , et donner quelque chose de moins precis que tes equations...

n°1521944
GrosBocdel
Posté le 28-02-2007 à 20:02:37  profilanswer
 

pedigus1 a écrit :

ou y a t il de grosses erreurs ailleurs?


pour l'instant ça calcule n'importe quoi, en tout cas.

Message cité 1 fois
n°1521977
el muchach​o
Comfortably Numb
Posté le 28-02-2007 à 21:48:10  profilanswer
 

GrosBocdel a écrit :

pour l'instant ça calcule n'importe quoi, en tout cas.


Comment le sais-tu ? Tu as visualisé les résultats ? Le mouvement du pendule est aléatoire, donc en lisant une suite de chiffres, je doute que tu puisses en conclure quoi que ce soit. A part visualiser l'animation, je ne vois pas trop comment on peut vérifier les résultats. Enfin si des points successifs sont disposés un peu partout quelque soit la taille de h, alors oui, c'est n'importe quoi. Il faut au moins qu'ils décrivent une trajectoire, mais après, pour vérifier cette trajectoire, sans visualisation, je ne vois pas trop. Peut-être qu'il y a moyen de vérifier si certains invariants sont conservés.

Message cité 2 fois
Message édité par el muchacho le 01-03-2007 à 06:08:09
n°1522002
pedigus1
Posté le 28-02-2007 à 22:48:22  profilanswer
 

GrosBocdel a écrit :


Code :
  1. void f(double x[], double y[])
  2. {
  3. y[0]=x[2];
  4. y[1]=x[3];
  5. y[2]=-(1/((m1+m2*sin(x[0]-x[1])*sin(x[0]-x[1]))*l1))*(m2*(-l1*x[2]*sin(x[0]-x[1])+g*sin(x[1]))*cos(x[0]-x[1])+m2*l2*x[3]*x[3]*sin(x[0]-x[1])+(m1+m2)*g*sin(x[0]));
  6. y[3]=-(1/l2)*(l1*(-(1/((m1+m2*sin(x[0]-x[1])*sin(x[0]-x[1]))*l1))*(m2*(-l1*x[2]*sin(x[0]-x[1])+g*sin(x[1]))*cos(x[0]-x[1])+m2*l2*x[3]*x[3]*sin(x[0]-x[1])+(m1+m2)*g*sin(x[0])))*cos(x[0]-x[1])-l1*x[2]*sin(x[0]-x[1])+g*sin(x[1]));
  7. }


 


 Je cherche depuis un moment mais je n'arrive pas à utiliser cette fonction.Je pense que sa pourrait beaucoup simplifier le code.
Je pense commencer par là: .
K1[0]=f(x[0],y[0]);
mais je vois pas la suite...

n°1522015
pedigus1
Posté le 28-02-2007 à 23:25:00  profilanswer
 

J'ai quand même retapé le code (sans utilisé cette fonction).Je pense que cette fois c'est la bonne méthode (enfin j'espere).Sa donne des resultats plus plausibles,à part pour ceux de x[1] qui me semblent bizarre.Mais est-ce que tu pourrais m'expliquer comment implementer la fonction dans mon code.Je poste le nouveau code:

Code :
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <math.h>
  4. #include "entetes.h"
  5. #include <SDL/SDL.h>
  7. double x[4];
  8. int m1,m2,l1,l2;
  10. int main(int argc, char *argv[])
  11. {
  12. double y[4],x1[4],y1[4];
  13. double K[4],L[4],M[4],N[4]; 
  14. float h;
  15. FILE *fi;
  16. int duree,tpsfin=0,tpsactuel=0;
  18. /*conditions intiales*/
  19. x[0]=10;
  20. x[1]=20;
  21. x[2]=30;
  22. x[3]=10;
  23. y[0]=10;
  24. y[1]=20;
  25. y[2]=30;
  26. y[3]=10;
  27. duree=100;
  28. m1=20;
  29. m2=10;
  30. l1=3;
  31. l2=6;
  32. h=0.05;
  34. fi=fopen("donnee.txt","w" );
  35. tpsfin = SDL_GetTicks()+duree;
  37. while(tpsactuel<tpsfin)
  38. {
  39. /*Runge-kutta*/
  40. fprintf(fi,"%f %f\n",x[0],x[1]);
  41.  
  42.   K[0]=fonctionf(x[0],x[1],x[2],x[3]);
  43.   L[0]=fonctiong(x[0],x[1],x[2],x[3]);
  44.   M[0]=fonctionm(x[2]);
  45.   N[0]=fonctionn(x[3]);
  46.  
  47.   K[1]=fonctionf(x[0]+(h/2)*M[0],x[1]+(h/2)*N[0],x[2]+(h/2)*K[0],x[3]+(h/2)*L[0]);
  48.   L[1]=fonctiong(x[0]+(h/2)*M[0],x[1]+(h/2)*N[0],x[2]+(h/2)*K[0],x[3]+(h/2)*L[0]);
  49.   M[1]=fonctionm(x[2]+(h/2)*K[0]);
  50.   N[1]=fonctionn(x[3]+(h/2)*L[0]);
  51.  
  52.   K[2]=fonctionf(x[0]+(h/2)*M[1],x[1]+(h/2)*N[1],x[2]+(h/2)*K[1],x[3]+(h/2)*L[1]);
  53.   L[2]=fonctiong(x[0]+(h/2)*M[1],x[1]+(h/2)*N[1],x[2]+(h/2)*K[1],x[3]+(h/2)*L[1]);
  54.   M[2]=fonctionm(x[2]+(h/2)*K[1]);
  55.   N[2]=fonctionn(x[3]+(h/2)*L[1]);
  56.  
  57.   K[3]=fonctionf(x[0]+h*M[2],x[1]+h*N[2],x[2]+h*K[2],x[3]+h*L[2]);
  58.   L[3]=fonctiong(x[0]+h*M[2],x[1]+h*N[2],x[2]+h*K[2],x[3]+h*L[2]);
  59.   M[3]=fonctionm(x[2]+h*K[2]);
  60.   N[3]=fonctionn(x[3]+h*L[2]);
  61.  
  62.   x1[2]=x[2]+(h/6)*(K[0]+2*K[1]+2*K[2]+K[3]);
  63.   x1[3]=x[3]+(h/6)*(L[0]+2*L[1]+2*L[2]+L[3]);
  64.   x1[0]=x[0]+(h/6)*(M[0]+2*M[1]+2*M[2]+M[3]);
  65.   x1[1]=x[1]+(h/6)*(N[0]+2*N[1]+2*N[2]+N[3]);
  66.  
  67.  
  68. /*injection des nouvelles valeurs dans les anciennes*/
  69. x[0]=x1[0];
  70. x[1]=x1[1];
  71. x[2]=x1[2];
  72. x[3]=x1[3];
  74. tpsactuel=SDL_GetTicks();
  75. }
  77. system("PAUSE" ); 
  78. return 0;
  79. }


 
Avec les fonctions:

Code :
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <math.h>
  5. double x[4];
  6. int m1,m2,l1,l2;
  7. float g=9.81;
  9. double fonctionm(double t)
  10. {
  11. double f;
  12. f=t;
  13. return f;
  14. }
  16. double fonctionn(double s)
  17. {
  18. double j;
  19. j=s;
  20. return j;
  21. }
  23. double fonctionf(double a,double b,double c,double d)
  24. {
  25. double F;
  26. F=-(1/((m1+m2*sin(a-b)*sin(a-b))*l1))*(m2*(-l1*c*sin(a-b)+g*sin(b))*cos(a-b)+m2*l2*d*d*sin(a-b)+(m1+m2)*g*sin(a));
  27. return F;
  28. }
  30. double fonctiong(double a,double b ,double c,double d)
  31. {
  32. double G;
  33. G=-(1/l2)*(l1*(-(1/((m1+m2*sin(a-b)*sin(a-b))*l1))*(m2*(-l1*c*sin(a-b)+g*sin(b))*cos(a-b)+m2*l2*d*d*sin(a-b)+(m1+m2)*g*sin(a)))*cos(a-b)-l1*c*sin(a-b)+g*sin(b));
  34. return G;
  35. }


Message édité par pedigus1 le 28-02-2007 à 23:31:27
n°1522053
el muchach​o
Comfortably Numb
Posté le 01-03-2007 à 06:11:22  profilanswer
 

Dans un premier temps, tu peux déjà tester ton intégrateur RK sur un pendule simple, avec plusieurs conditions initiales différentes.
Les équations du mvt sont plus simples.
http://fr.wikipedia.org/wiki/Pendule_pesant
http://www.chimix.com/an4/term4/pendule.htm

 

Tu devrais obtenir des mouvements sur un arc de cercle ou un cercle entier.

 

Ah, et puis découpe le calcul de tes fonctions pour qu'elles tiennent sur plusieurs lignes, histoire de les rendre plus digestes, avec des noms qui ont un sens.
Définis des variables temporaires pour sin(theta1-theta2), cos(theta1-theta2), sin(theta1) et sin(theta2), ça évite les recalculs (enfin je pense que le compilo optimise, mais c'est mieux quand même).

 

Ensuite, si ça marche, tu peux essayer les méthodes avec pas adaptatif ou l'algo de Beeman http://en.wikipedia.org/wiki/Beeman%27s_algorithm
(tiens je ne savais pas que l'intégration de Verlet venait de Loup Verlet, un de mes profs à Orsay).


Message édité par el muchacho le 01-03-2007 à 06:37:52
n°1522055
GrosBocdel
Posté le 01-03-2007 à 06:59:25  profilanswer
 

el muchacho a écrit :

Comment le sais-tu ? Tu as visualisé les résultats ? Le mouvement du pendule est aléatoire, donc en lisant une suite de chiffres, je doute que tu puisses en conclure quoi que ce soit. A part visualiser l'animation, je ne vois pas trop comment on peut vérifier les résultats. Enfin si des points successifs sont disposés un peu partout quelque soit la taille de h, alors oui, c'est n'importe quoi. Il faut au moins qu'ils décrivent une trajectoire, mais après, pour vérifier cette trajectoire, sans visualisation, je ne vois pas trop. Peut-être qu'il y a moyen de vérifier si certains invariants sont conservés.


 
Juste que les deux angles ne restent pas entre 0 et 2Pi, et que ça diverge aux temps longs...
Je propose qu'on arrète de faire comme si on n'avait pas déjà ce programme
http://www.physics.usyd.edu.au/~wheat/dpend_html/
 
 

n°1522060
esox_ch
Posté le 01-03-2007 à 07:40:01  profilanswer
 

el muchacho a écrit :

Comment le sais-tu ? Tu as visualisé les résultats ? Le mouvement du pendule est aléatoire, donc en lisant une suite de chiffres, je doute que tu puisses en conclure quoi que ce soit. A part visualiser l'animation, je ne vois pas trop comment on peut vérifier les résultats. Enfin si des points successifs sont disposés un peu partout quelque soit la taille de h, alors oui, c'est n'importe quoi. Il faut au moins qu'ils décrivent une trajectoire, mais après, pour vérifier cette trajectoire, sans visualisation, je ne vois pas trop. Peut-être qu'il y a moyen de vérifier si certains invariants sont conservés.


 
Dans ce système c'est l'énéergie mécanique qui est conservée (suffi de voir la forme de l'hamiltonien pour s'en convaincre), par contre verifier la conservation d'énergie à partir de points ... Je vois pas trop comment faire ...

Message cité 1 fois
Message édité par esox_ch le 01-03-2007 à 07:40:19

---------------
Si la vérité est découverte par quelqu'un d'autre,elle perd toujours un peu d'attrait
n°1522062
GrosBocdel
Posté le 01-03-2007 à 07:43:56  profilanswer
 

esox_ch a écrit :

Dans ce système c'est l'énéergie mécanique qui est conservée (suffi de voir la forme de l'hamiltonien pour s'en convaincre), par contre verifier la conservation d'énergie à partir de points ... Je vois pas trop comment faire ...


 
Suffit de vérifier que les deux équa diff restent égales à zéro pendant le calcul

n°1522299
pedigus1
Posté le 01-03-2007 à 15:03:47  profilanswer
 

Sympa ce site GrosBoc ,mais je n'ai pris que les equations simplifiées.
Pour le fait que sa dépasse 2Pi , je pense que c'est normal parceque je n'ai pas encore converti mes angles en radian.
Sinon j'avais une autre question, pour mes conditions initiales, je dois specifier x[0] et x[1] (les deux théta), mais est-ce que je dois mettre les x[2] et [3] et tous les y[] nuls (car ils sont les dérivées de mes angles de départs qui sont constants)?
J'ai fait par exemple un test avec x[0]=50 et x[1]=40 et tout le reste nul, et je vois que x[0] augmente dès le début.Ce n'est pas normal, le pendule devrait descendre.Et si j'initialise tous les x[] et les y[] à 0 ,mon fichier de sortie n'est composé que de zeros (normal).

Message cité 1 fois
Message édité par pedigus1 le 01-03-2007 à 15:09:03
n°1522451
GrosBocdel
Posté le 01-03-2007 à 18:12:32  profilanswer
 

pedigus1 a écrit :

Sympa ce site GrosBoc ,mais je n'ai pris que les equations simplifiées.


C'est ton RK qui n'est pas bon. Compare ton programme avec celui tu site web.
 

pedigus1 a écrit :

Pour le fait que sa dépasse 2Pi , je pense que c'est normal parceque je n'ai pas encore converti mes angles en radian.


 
J'ai converti tes angles en radian en compilant ton programme  
 

pedigus1 a écrit :


Sinon j'avais une autre question, pour mes conditions initiales, je dois specifier x[0] et x[1] (les deux théta), mais est-ce que je dois mettre les x[2] et [3] et tous les y[] nuls (car ils sont les dérivées de mes angles de départs qui sont constants)?


Tes y ne sont utilisés nulle part. Sinon tes valeurs initiales, tu en fais ce que tu veux. Tu peux très bien avoir envie de pousser ton pendule au démarrage.
 

pedigus1 a écrit :


J'ai fait par exemple un test avec x[0]=50 et x[1]=40 et tout le reste nul, et je vois que x[0] augmente dès le début.Ce n'est pas normal, le pendule devrait descendre.Et si j'initialise tous les x[] et les y[] à 0 ,mon fichier de sortie n'est composé que de zeros (normal).


Demande toi où est ton pendule quand "il est à 50" et ce que prennent comme arguments les fonctions trigo.
 

Message cité 1 fois
n°1522458
pedigus1
Posté le 01-03-2007 à 18:29:18  profilanswer
 

GrosBocdel a écrit :

C'est ton RK qui n'est pas bon. Compare ton programme avec celui tu site web.
 
J'ai converti tes angles en radian en compilant ton programme  
 
Demande toi où est ton pendule quand "il est à 50" et ce que prennent comme arguments les fonctions trigo.


 
D'accord, je vois pour mon 50 il faut que je le passe en radian avant.
 
Mais pour ce qui est du code du site web, j'ai passé du temps dessus et je vois pas vraiment comment il procède.il utilise des choses que je n'ai pas encore vues.Tu pourrais me dire ou est mon erreur pour runge kutta(ou alors si tout est faux).(le code du site web ne marche pas chez moi, donne une erreur à l'execution du programme).
EDIT: j'ai corrigé quelques trucs dans le code et je trouve des valeurs coherentes (j'ai testé avec le pendule au repos qui reste bien au repos et pour les autres conditions intiales j'ai quelquechose de continu).Je pars pour une semaine, je finirais le code en rentrant.
Merci de votre aide
 
EDIT:
C'est bon sa marche maintenant.
Merci de votre aide

Message cité 1 fois
Message édité par pedigus1 le 14-03-2007 à 12:59:01
n°1532947
GrosBocdel
Posté le 23-03-2007 à 17:37:12  profilanswer
 

pedigus1 a écrit :


EDIT:
C'est bon sa marche maintenant.
Merci de votre aide


 
T'es sûr? T'as changé quoi et ça donne quoi?

Message cité 1 fois
n°1532982
pedigus1
Posté le 23-03-2007 à 19:48:04  profilanswer
 

GrosBocdel a écrit :

T'es sûr? T'as changé quoi et ça donne quoi?


 
Oula je pourrais plus te dire, sa remonte, mais on a beaucoup avancé depuis.
Je poste un lien vers le projet complet (qui marche) (il faut avoir la SDL et SDL-ttf d'installé ,ils sont inclus dans le zip)
 
http://rapidshare.com/files/224374 [...] _.zip.html

Message cité 1 fois
n°1533111
GrosBocdel
Posté le 24-03-2007 à 10:18:59  profilanswer
 

pedigus1 a écrit :

Oula je pourrais plus te dire, sa remonte, mais on a beaucoup avancé depuis.
Je poste un lien vers le projet complet (qui marche) (il faut avoir la SDL et SDL-ttf d'installé ,ils sont inclus dans le zip)
 
http://rapidshare.com/files/224374 [...] _.zip.html


 
Ok d'ac. Je vais regarder ça et aller jouer avec dans mon coin  :)  
 

n°1533281
pedigus1
Posté le 24-03-2007 à 22:37:42  profilanswer
 

lol ok

mood
Publicité
Posté le   profilanswer
 


Aller à :
Ajouter une réponse
  FORUM HardWare.fr
  Programmation
  C

  [C] mouvement d'un double pendule

 

Sujets relatifs
incrementer un calendrier avec une precision double.page executee en double et redirection
Double tri dans un Array pour un classementprobleme d'overflow avec variables type Double [résolu]
[C#] [Résolu] Comment communiquer avec une dll c++[C] Liste chaibee Double Probleme
mettre un null à la place d'un double[C] Des accolades "just pour le fun" ?
Variable et double menu deroulantConditions multiples ( double WHERE...)
Plus de sujets relatifs à : [C] mouvement d'un double pendule

Forum MesDiscussions.Net, Version 2010.2
(c) 2000-2011 Doctissimo
Page générée en 0.126 secondes

Copyright © 1997-2025 Groupe LDLC (Signaler un contenu illicite / Données personnelles)