Appliquer Runge-Kutta 4 à un pendule élastique.

neokiller007 · 21/01/2010 - 12h05

Bonjour,

J'ai crée un programme C qui me permet de simuler l'évolution du pendule élastique (pendule avec un ressort au lieu d'une corde).

Pour cela j'ai utilisé la méthode d'Euler:
$f(a+h)=f(a)+hf'(a)+o(h)$
Donc on approxime en enlevant le o(h)

Du coup on a:
$x(t_0+\Delta t)=x(t_0)+\Delta t v(t_0)$
$v(t_0+\Delta t)=v(t_0)+\Delta t a(t_0)$

Et j'applique l'algorithme suivant:
On connait x(t) et v(t)(conditions initiales), on calcul a(t) par le PFD (pour cela on a besoin de x(t) et v(t) à cause de la force de frottement). Puis on calcul $x(t+\Delta t)$ grace à euler, puis $v(t+\Delta t)$ et on recommence.

Mais le problème c'est que les erreur numériques sont trop grandes et amplifient le mouvement quand on ne met pas de frottements.

Donc j'aimerais utiliser une autre méthode.
J'ai pensé à la méthode de Verlet qui nous donne:
$x(t_0+\Delta t)=x(t_0)+\Delta t v(t_0)+\frac12 (\Delta t)^2 a(t_0)$
$v(t_0+\Delta t)=v(t_0)+\Delta t \frac{a(t_0+\Delta t)+a(t_0)}{2}$

Mais le problème c'est que pour calculer $v(t_0+\Delta t)$ j'ai besoin de $a(t_0+\Delta t)$ . Et pour calculer $a(t_0+\Delta t)$ j'ai besoin de $x(t_0+\Delta t)$ ça pas de problème je le calcul avec x(t0) et v(t0), mais j'ai aussi besoin de $v(t_0 + \Delta t)$ ! (à cause des forces de frottements) donc le problème est insoluble

J'ai donc pensé à la méthode de runge-kutta d'ordre4:
Sous wikipédia on donne les formules ici http://fr.wikipedia.org/wiki/M%C3%A9...de_Runge-Kutta
Mais je ne comprend pas à quoi correspond $y_n$ , $y_{n+1}$ et $f(t_n, y_n)$ .
f est une fonction à deux variables alors que moi je n'ai que des fonctinos d'une variable donc je comprend pas.

Est-ce que quelqu'un pourrait m'expliciter ça en terme de f(a+h)=quelque chose ou alors en terme de x(t), v(t) et $\Delta t$

Merci.

**LPFR** · 21/01/2010 - 12h39

Bonjour.
Je viens de relire ma bible de calcul numérique "Numerical recipes".
Et c'est plus clair que wikipedia (rien de surprenant!).
En fait la fonction f(x, y) de wikipedia est écrite dans mon bouquin f'(x, y). C'est bêtement la dérivée de la fonction au point x.
Les k1, .., k4 sont les dérivées aux points indiqués.
Au revoir.

neokiller007 · 21/01/2010 - 12h52

Donc si je résume avec des notations compréhensible ça nous donne:

$f(a+h)=f(a)+\frac{h}{6}(k_1+2k _2+2k_3+k_4)$
avec:
$k_1=f'(a)$
$k_2=f'(a+\frac{h}{2}k_1)$
$k_3=f'(a+\frac{h}{2}k_2)$
$k_4=f'(a+hk_3)$

Est-ce correcte ?

**LPFR** · 21/01/2010 - 12h55

Re.
Oui, c'est ça.
A+

philou21 · 21/01/2010 - 14h29

Bonjour
Verlet c'est pas :
$x(t+\Delta t)=2x(t)-x(t-\Delta t)+ a(t)(\Delta t)^2$
qui laisse une erreur en $(\Delta t)^4$ ?

neokiller007 · 21/01/2010 - 14h58

Envoyé par philou21

Bonjour
Verlet c'est pas :
$x(t+\Delta t)=2x(t)-x(t-\Delta t)+ a(t)(\Delta t)^2$
qui laisse une erreur en $(\Delta t)^4$ ?

Je ne sais pas, j'ai trouvé cette formule sur internet.
En fait c'est la même chose qu'Euler mais on fait un DL à l'ordre 2.
Puis on approxime a'(t0) par son taux d'accroissement en $\Delta t$ et on obtient cette expression.

Mais si tu me dis comment trouver ton expression je serais ravis de l'utiliser.

philou21 · 21/01/2010 - 15h32

Oui, facile :
tu developpes au 3ème ordre d'un coté :
$x(t + h) = x(t) + v(t)h + {\textstyle{1 \over 2}}a(t){h^2} + {\textstyle{1 \over 6}}{x^{'''}}(t){h^3} + O({h^4})$

et de l'autre (en remontant le temps) :
$x(t - h) = x(t) - v(t)h + {\textstyle{1 \over 2}}a(t){h^2} - {\textstyle{1 \over 6}}{x^{'''}}\left( t \right){h^3} + O\left( {{h^4}} \right)$

tu additionnes les deux expressions et tu obtiens :
$x(t + h) = 2x(t) - x(t - h) + a(t){h^2} + O({h^4})$

les vitesses et le troisième ordre s'annulent. Tu as donc une précision jusqu'au 3ème ordre.

Mais il faut que tu calcules le premier point $x(t_0+h)$ pour pouvoir enclencher le processus.
Tu peux le faire avec :
$x(t_0 + h) = x(t_0) + v(t_0)h + {\textstyle{1 \over 2}}a(t_0){h^2}$

Le seul problème est que tu n'as pas les vitesses, ça peut être embêtant si tu as des frottements. Tu peux toujours l'estimer par :
$v(t) = (x(t) - x(t - h))/h$

neokiller007 · 21/01/2010 - 15h51

Merci pour la démo.

Effectivement c'est embêtant pour la vitesse.
On calcul avec précision x mais on fait une erreur grossière sur la vitesse :/

philou21 · 21/01/2010 - 16h20

Envoyé par neokiller007

Merci pour la démo.

Effectivement c'est embêtant pour la vitesse.
On calcul avec précision x mais on fait une erreur grossière sur la vitesse :/

Oui et non ... l'erreur ( en O(h²)) ne s'accumule pas trop puisque la vitesse est calculée à chaque fois à partir de valeurs assez précises (les x).
essaye...
C'est facile à programmer et ça donne en général de bon résultats

remarque : de plus, souvent les forces de frottements dépendent de v de façon un peu empirique, je ne suis pas persuadé qu'une valeur très précise soit absolument indispensable.
si tu n'as pas besoin de v(t) pour calculer a(t) mais juste pour information il vaudra mieux, bien sûr, calculer la vitesse a posteriori par :
$v(t) = (x(t + h) - x(t - h)/2h$

neokiller007 · 22/01/2010 - 10h06

Effectivement,c'est quand même beaucoup plus stable qu'avec euler.

philou21 · 22/01/2010 - 10h12

Envoyé par neokiller007

Effectivement,c'est quand même beaucoup plus stable qu'avec euler.

Ça ne m'étonne pas... (de toute façon il faut régler Δt de manière à avoir une stabilité)
Tu as mis un peu de dissipation ? sous quelle forme ?

neokiller007 · 22/01/2010 - 11h40

Frottements visqueux en régime turbulent (car dans l'air on atteins très vite le régime turbulent)

philou21 · 22/01/2010 - 12h40

Envoyé par neokiller007

Frottements visqueux en régime turbulent (car dans l'air on atteins très vite le régime turbulent)

de toute façon, je pense que tu résous les équations pour l et thêta et que les accélérations comportent forcement des termes de vitesses.

neokiller007 · 22/01/2010 - 13h42

Envoyé par philou21

de toute façon, je pense que tu résous les équations pour l et thêta et que les accélérations comportent forcement des termes de vitesses.

Non non, sans frottement y a pas de vitesse dans l'accélération.

J'ai une question, pourquoi la plupart des logiciels scientifiques sont codés en Fortran ?
Je trouve ce langage horrible et bien moins pratique que le C/C++. Et niveau performance c'est le C/C++ qui remporte.
De plus le C/C++ sont très utilisés donc plus facile de trouvé de la doc et des bibliothèques (openGL, SDL que j'utilise, SFML...) notamment des GUI comme Qt etc...

philou21 · 22/01/2010 - 13h52

Envoyé par neokiller007

Non non, sans frottement y a pas de vitesse dans l'accélération.

J'ai une question, pourquoi la plupart des logiciels scientifiques sont codés en Fortran ?
Je trouve ce langage horrible et bien moins pratique que le C/C++. Et niveau performance c'est le C/C++ qui remporte.
De plus le C/C++ sont très utilisés donc plus facile de trouvé de la doc et des bibliothèques (openGL, SDL que j'utilise, SFML...) notamment des GUI comme Qt etc...

A bon ?
j'avais souvenir que dans $\ddot l$ il y avait un terme en ${\dot \theta }$
et pareil pour ${\ddot \theta }$

neokiller007 · 22/01/2010 - 14h02

Envoyé par philou21

A bon ?
j'avais souvenir que dans $\ddot l$ il y avait un terme en ${\dot \theta }$
et pareil pour ${\ddot \theta }$

Un petit schéma:

- $\vec{P}=mg\vec{j}$
- $\vec{F}=-k \Delta l \vec{Tu}=-k \Delta l (sin \theta \vec{i}+ cos \theta \vec{j})$ avec: $l_0$ la longueur du ressort à vide et $\Delta l=l-l_0$ l'allongement du ressort par rapport à la longueur à vide.
- $\vec{f}=-C||\vec{v}|| \vec{v}$
PFD:
$m\vec{a}=\vec{P}+\vec{F}+\vec{ f}$ avec: $\vec{a}=(a_x,a_y)=\frac{d \vec{v}}{dt}=\frac{d^2u}{dt^2}$
Donc
$m(a_x \vec{i} + a_y \vec{j})=mg\vec{j}-k \Delta l (sin \theta \vec{i}+ cos \theta \vec{j})-Cv(v_x \vec{i}+ v_y \vec{j})$ $\Leftrightarrow \left\{ \begin{array}{ll} a_x=-\frac{k}{m} \Delta l sin \theta - C \frac{v v_x}{m} \\ a_y=g-\frac{k}{m} \Delta l cos \theta - C \frac{v v_y}{m} \end{array} \right$

$\mbox{Avec: } \left\{ \begin{array}{llll} sin \theta=\frac{u_x-T_x}{l} \\ cos \theta=\frac{u_y-T_y}{l} \\ l=\sqrt{(u_x-T_x)^2+(u_y-T_y)^2}\\ v=\sqrt{v_x^2+v_y^2} \end{array} \right$

On voit bien que les vitesse ne sont là qu'a cause des frottements.

philou21 · 22/01/2010 - 14h06

En ce qui concerne le fortran on va pas relancer ce vieux débat...

Toujours est-il qu'ils existent des millions de lignes codés en fortran et on ne va pas recoder à chaque fois que sort un nouveau langage...
Contrairement à ce que tu penses, le fortran génère un code très efficace (depuis le temps, le compilateur est plutôt performant) de plus les bibliothèques sont extrêmement bien fournies.
Tu sais, ce qu'exécute une machine c'est du code machine, c'est pas du C ou du fortran ou autres...

philou21 · 22/01/2010 - 14h12

Envoyé par neokiller007

...On voit bien que les vitesse ne sont là qu'a cause des frottements.

Oui, oui, OK
je parlais de résolution en ${\ddot l}$ et en ${\ddot \theta }$
pas en ${\ddot x}$ et ${\ddot y}$

neokiller007 · 22/01/2010 - 14h17

Envoyé par philou21

Oui, oui, OK
je parlais de résolution en ${\ddot l}$ et en ${\ddot \theta }$
pas en ${\ddot x}$ et ${\ddot y}$

Ok, mais je vois pas comment coder le truc si tu résouds en ${\ddot l}$ et en ${\ddot \theta }$

philou21 · 22/01/2010 - 14h31

Envoyé par neokiller007

Ok, mais je vois pas comment coder le truc si tu résouds en ${\ddot l}$ et en ${\ddot \theta }$

Non, si tu veux inclure le frottement, je pense qu'il vaut mieux faire en cartésien comme tu as fait.
(autrement c'est lagrangien...)

neokiller007 · 22/01/2010 - 14h58

intéressant je savais pas que c'était possible.

Appliquer Runge-Kutta 4 à un pendule élastique.

Appliquer Runge-Kutta 4 à un pendule élastique.

Re : Appliquer Runge-Kutta 4 à un pendule élastique.

Re : Appliquer Runge-Kutta 4 à un pendule élastique.

Re : Appliquer Runge-Kutta 4 à un pendule élastique.

Re : Appliquer Runge-Kutta 4 à un pendule élastique.

Re : Appliquer Runge-Kutta 4 à un pendule élastique.

Re : Appliquer Runge-Kutta 4 à un pendule élastique.

Re : Appliquer Runge-Kutta 4 à un pendule élastique.

Re : Appliquer Runge-Kutta 4 à un pendule élastique.

Re : Appliquer Runge-Kutta 4 à un pendule élastique.

Re : Appliquer Runge-Kutta 4 à un pendule élastique.

Re : Appliquer Runge-Kutta 4 à un pendule élastique.

Re : Appliquer Runge-Kutta 4 à un pendule élastique.

Re : Appliquer Runge-Kutta 4 à un pendule élastique.

Re : Appliquer Runge-Kutta 4 à un pendule élastique.

Re : Appliquer Runge-Kutta 4 à un pendule élastique.

Re : Appliquer Runge-Kutta 4 à un pendule élastique.

Re : Appliquer Runge-Kutta 4 à un pendule élastique.

Re : Appliquer Runge-Kutta 4 à un pendule élastique.

Re : Appliquer Runge-Kutta 4 à un pendule élastique.

Re : Appliquer Runge-Kutta 4 à un pendule élastique.

Discussions similaires

Runge Kutta

Méthode Runge-Kutta heeeelp !

Runge Kutta en méca céleste

Runge kutta pour la RFD

Runge Kutta (et Fortran)