mouvement d'un double pendule [C] - C - Programmation
Marsh Posté le 25-02-2007 à 16:33:42
http://www.cpt.univ-mrs.fr/~grapperon/oldpdf/mip.pdf
Marsh Posté le 25-02-2007 à 18:22:48
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.
Marsh Posté le 25-02-2007 à 23:21:45
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.
Marsh Posté le 26-02-2007 à 15:38:43
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?
Marsh Posté le 26-02-2007 à 17:36:12
pedigus1 a écrit : |
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 : |
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 : |
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
Marsh Posté le 26-02-2007 à 20:29:31
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 :
|
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?
Marsh Posté le 26-02-2007 à 22:47:08
pedigus1 a écrit : |
Faut penser à initialiser x avant et à ne pas oublier les * des multiplications.
pedigus1 a écrit : |
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.
Marsh Posté le 27-02-2007 à 14:10:40
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.
Marsh Posté le 27-02-2007 à 20:10:29
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].
Marsh Posté le 27-02-2007 à 20:20:40
Voila l'ébauche du code:
Code :
|
Avec ces fonctions (qui d'ailleur devraient me permettre de supprimer mon paragraphe de /*à t=0*/ non?):
Code :
|
Marsh Posté le 27-02-2007 à 20:43:15
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
Marsh Posté le 27-02-2007 à 21:31:52
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
Marsh Posté le 28-02-2007 à 06:58:18
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? |
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.
Marsh Posté le 28-02-2007 à 14:08:07
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?
Marsh Posté le 28-02-2007 à 18:14:24
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 :
|
Avec le fichier des fonctions:
Code :
|
En tous cas, merci beaucoup pour ton aide !!!
Marsh Posté le 28-02-2007 à 18:24:46
pedigus1 a écrit : Je n'ai pas beaucoup de C derriere moi mais j'y travaille. |
Tu définis juste ta fonction
Code :
|
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.
Marsh Posté le 28-02-2007 à 18:28:25
ReplyMarsh Posté le 28-02-2007 à 18:30:01
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.
Marsh Posté le 28-02-2007 à 18:38:13
red faction a écrit : mon dieu que cest compliqué ! |
Tu as mieux à proposer? parce que c'est quand même un grand moment de solitude pour moi depuis le début, là !
Marsh Posté le 28-02-2007 à 18:40:49
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 )
Marsh Posté le 28-02-2007 à 18:50:21
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...
Marsh Posté le 28-02-2007 à 20:02:37
pedigus1 a écrit : ou y a t il de grosses erreurs ailleurs? |
pour l'instant ça calcule n'importe quoi, en tout cas.
Marsh Posté le 28-02-2007 à 21:48:10
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.
Marsh Posté le 28-02-2007 à 22:48:22
GrosBocdel a écrit :
|
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...
Marsh Posté le 28-02-2007 à 23:25:00
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 :
|
Avec les fonctions:
Code :
|
Marsh Posté le 01-03-2007 à 06:11:22
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).
Marsh Posté le 01-03-2007 à 06:59:25
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/
Marsh Posté le 01-03-2007 à 07:40:01
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 ...
Marsh Posté le 01-03-2007 à 07:43:56
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
Marsh Posté le 01-03-2007 à 15:03:47
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).
Marsh Posté le 01-03-2007 à 18:12:32
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 : |
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 : |
Demande toi où est ton pendule quand "il est à 50" et ce que prennent comme arguments les fonctions trigo.
Marsh Posté le 01-03-2007 à 18:29:18
GrosBocdel a écrit : C'est ton RK qui n'est pas bon. Compare ton programme avec celui tu site web. |
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
Marsh Posté le 23-03-2007 à 17:37:12
pedigus1 a écrit : |
T'es sûr? T'as changé quoi et ça donne quoi?
Marsh Posté le 23-03-2007 à 19:48:04
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
Marsh Posté le 24-03-2007 à 10:18:59
pedigus1 a écrit : Oula je pourrais plus te dire, sa remonte, mais on a beaucoup avancé depuis. |
Ok d'ac. Je vais regarder ça et aller jouer avec dans mon coin
Marsh Posté le 25-02-2007 à 15:32:26
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