Appliquer Runge-Kutta(RK4) aux équations du mouvement d'un corps rigide

**johnstanp** · 07/05/2008, 01h31

Cela fait un petit moment que je cherche à comprendre comment implémenter Runge-Kutta pour le système d'équations de Newton.
Voici le système d'équations :

dx(t)/dt = v(t)
dR(t)/dt = w(t)*R(t)
dP(t)/dt = F(t)
dL(t)/dt = T(t)

x(t) : position du centre de masse
R(t) : matrice d'orientation
P(t) : moment linéaire
L(t) : moment angulaire

v(t) : vitesse du centre de masse
w(t) : vitesse angulaire
F(t) : somme des forces appliquées sur le corps
T(t) : somme des couples générés par les forces appliquées

Les grandeurs auxilliaires se calculent comme suit:

v(t) = P(t)/M
I(t) = R(t) Ib R(t)^T
w(t) = I(t)^-1 L(t)

I(t) = tenseur d'inertie dans le repère fixe choisi
Ib = tenseur d'inertie dans le repère local
R(t)^T = transposée de la matrice d'orientation
I(t)^-1 = inverse du tenseur d'inertie exprimé dans le repère fixe

En fait , je reprends ces équations du papier de Baraff "An introduction to physically based modeling" disponible à l'adresse suivante :

http://www.cs.cmu.edu/~baraff/sigcourse/notesd1.pdf

J'ai écrit une ébauche de moteur physique utilisant la méthode d'Euler(pour sa trop grande simplicité) : il s'avère que je n'obtiens pas de résultats satisfaisants lorsque je simule freinage et virage simultanément pour un mobile...
J'aimerais donc utiliser une méthode plus précise d'où celle de Runge-Kutta d'ordre 4. Le seul problème vient du fait que dans celle-ci , on doit disposer d'une fonction dépendant explicitement du temps et des variables que l'on veut intégrer.Il s'avère que dans l'équation donnée , je ne dispose pas d'une fonction me permettant de calculer le membre de droite en fonction du temps et des variables x(t) , R(t) , P(t) et L(t) , puisque je ne connais pas à priori ce qu'elles valent aux instants t = t0 + dt/2 et t = t0 + dt.
Les seules données dont je dispose sont les forces et couples générés à l'instant t0( que je peux supposer constantes entre t0 et t0 + dt ou faire varier linéairement , exponentiellement , etc ). Il n'empêche que je ne peux calculer le membre de droite en un instant autre que l'instant t0.

Quelqu'un pourrait-il me dire comment appliquer Runge-Kutta dans ce cas-ci?
Merci pour vos réponses.

**millie** · 07/05/2008, 12h01

Normalement, runge kutta est juste fait pour résoudre des équations du type :

dy/dx = f(x,y) (où f est constant !)

C'est à dire qu'il faut que tu essayes par exemple d'avoir une équation sous la forme :
dx(t)/dt = f(x(t), t) (avec f constant)

Mais apparemment ce n'est pas vraiment fait pour ton cas car ça ne semble pas possible a priori

(je n'ai pas plus creusé que ça cela dit)

**FR119492** · 07/05/2008, 15h19

Salut!
En fait, je ne vois pas où est ton problème: au début de chaque pas d'intégration, tu connais x(t), R(t), P(t) et L(t), la première fois comme conditions initiales, et par la suite comme résultat du pas précédent. A partir de cela, tu as tout ce qu'il te faut pour calculer les termes de droite dans tes équations différentielles.
Voici un code Fortran de la méthode RK4 que j'ai retrouvé dans mes vieilleries. Tu peux sans autres le traduire dans le langage de ton choix.

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

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
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
      Subroutine D006(SM,N,T,DT,Y)
C
C***********************************************************************
C
C     Bibliothèque JMBPACK
C
C     Série D: Equations différentielles
C
C***********************************************************************
C
C     Méthode de Runge-Kutta d'ordre 4 (1/6 + 2/6 + 2/6 + 1/6)
C
C***********************************************************************
C
C     Version 3.0: Jean-Marc Blanc, juillet 2004
C
C **********************************************************************
C
C     SM     Nom d'un sous-programme fourni par l'utilisateur pour le
C            calcul des dérivées à partir des valeurs de la variable
C            indépendante et des variables d'état. Ce sous-programme
C            doit être de la forme
C              Subroutine SM(N,T,Y,DY)
C                DY(1)= ...
C                DY(2)= ...
C                 .
C                 .
C                DY(N)= ...
C                Return
C              End
C            et avoir été déclaré dans une instruction External.
C
C     N      Ordre du système différentiel.
C
C     T      Variable indépendante.
C
C     DT     Pas d'intégration.
C
C     Y      Vecteur des variables d'état.
C
C **********************************************************************
C
      Implicit None
C
      Integer N
      Real*8 DT,T,Y(N)
C
      Integer I
      Real*8 DY(N),K1(N),K2(N),K3(N),K4(N),Z(N)
C
      Call SM(N,T,Y,DY)
      Do I=1,N
        K1(I)=DT*DY(I)
        Z(I)=Y(I)+K1(I)/2.d0
      End Do
C
      Call SM(N,T+DT/2.d0,Z,DY)
      Do I=1,N
        K2(I)=DT*DY(I)
        Z(I)=Y(I)+K2(I)/2.d0
      End Do
C
      Call SM(N,T+DT/2.d0,Z,DY)
      Do I=1,N
        K3(I)=DT*DY(I)
        Z(I)=Y(I)+K3(I)
      End Do
C
      Call SM(N,T+DT,Z,DY)
      Do I=1,N
        K4(I)=DT*DY(I)
        Y(I)=Y(I)+(K1(I)+2.d0*K2(I)+2.d0*K3(I)+K4(I))/6.d0
      End Do
C
      T=T+DT
      Return
C
      End

Jean-Marc Blanc

**millie** · 07/05/2008, 15h24

Voici une version Java http://www.developpez.net/forums/sho....php?p=2528269
Et C++ http://www.developpez.net/forums/sho...77&postcount=2

mais je travaillais avec une équation du type y'(x)= f(x,y(x)) avec f constant.

**acx01b** · 07/05/2008, 19h40

salut je n'ai pas bien compris, car si tes forces sont constantes par morceaux alors c'est facile, pas besoin d'équa diff, pour un moteur physique ça suffira à priori non ?

tes forces sont constantes par morceaux ?

le problème sera pour les collisions, mais pas pour les équations newton