Chapitre 6 : Simulation physique
Alexandre Blondin Massé
Laboratoire d’informatique formelle
Université du Québec à Chicoutimi
12 juin 2014
Cours 8INF802
Département d’informatique et mathématique
A. Blondin Massé (UQAC) 12 juin 2014 1 / 51
Table des matières
1. Introduction
2. Mécanique des fluides
3. Pendule simple
4. Pendule double
5. Particules dans une boîte
6. Montagnes russes
A. Blondin Massé (UQAC) 12 juin 2014 2 / 51
Table des matières
1. Introduction
2. Mécanique des fluides
3. Pendule simple
4. Pendule double
5. Particules dans une boîte
6. Montagnes russes
A. Blondin Massé (UQAC) 12 juin 2014 3 / 51
Simulations physiques
I Dans certaines situations, la simulation repose sur des
principes physiques;
I Plusieurs outils de domaines mathématiques sont
fondamentaux :
I Algèbre linéaire;
I Calcul différentiel et intégral;
I Équations différentielles;
I Analyse numérique.
A. Blondin Massé (UQAC) 12 juin 2014 4 / 51
Bibliothèques
I Pour tous les langages de programmation populaire, il existe
évidemment des bibliothèques offrant ce genre d’outils;
I C++ : Box2D, Chipmunk, Havoc, Bullet, Vortex, etc.
I Java : JBox2D, dyn4j, JBullet, Falling2D, etc.
I C# : Farseer, Box2Dx, [Link], JigLib, Jitter,
Henge3D, etc.
I Python : PyODE, pymunk, Pygame, pybox2d, Panda3D,
Matplotlib, etc.
A. Blondin Massé (UQAC) 12 juin 2014 5 / 51
Scipy
I Scipy est un ensemble d’outils logiciels facilitant la
programmation scientifique;
I Il inclut plusieurs paquets, tels que
I Numpy, pour le calcul numérique;
I Matplotlib, pour créer des dessins et des animations;
I Sympy, pour le calcul symbolique,
I etc.
A. Blondin Massé (UQAC) 12 juin 2014 6 / 51
Numpy
I Numpy est le noyau de Scipy qui prend en charge les calculs
numériques;
I Il offre entre autres :
I Des tableaux multidimensionnels efficaces;
I Plusieurs fonctions pour manipuler les tableaux;
I La transformée de Fourier discrète;
I Des fonctions d’algèbre linéaire;
I Des fonctions mathématiques de base.
A. Blondin Massé (UQAC) 12 juin 2014 7 / 51
Matplotlib
I Matplotlib offre différents services de dessins et
d’animations;
I On peut entre autres créer
I Des graphiques de fonctions;
I Des histogrammes;
I Des graphiques en 3D;
I Des champs de vecteurs;
I Des diagrammes circulaires.
A. Blondin Massé (UQAC) 12 juin 2014 8 / 51
Exemple (1/3)
I Animation d’une vague : basic_animation.py;
I On importe les modules nécessaires :
# Imports
import numpy as np
from matplotlib import pyplot as plt
from matplotlib import animation
I On configure l’image :
# First set up the figure, the axis, and the plot element
# we want to animate
fig = [Link]()
ax = [Link](xlim=(0, 2), ylim=(-2, 2))
line, = [Link]([], [], ’b’, lw=2)
A. Blondin Massé (UQAC) 12 juin 2014 9 / 51
Exemple (2/3)
I Pour créer une animation, il faut une fonction
d’initialisation :
# initialization function: plot the background of each
# frame
def init():
line.set_data([], [])
return line,
I Puis une fonction principale pour l’animation :
# animation function. This is called sequentially
def animate(i):
x = [Link](0, 2, 1000)
y = [Link](2 * [Link] * (x - 0.01 * i))
line.set_data(x, y)
return line,
A. Blondin Massé (UQAC) 12 juin 2014 10 / 51
Exemple (3/3)
I On lance l’animation :
# call the animator. blit=True means only re-draw
# the parts that have changed.
anim = [Link](fig,
animate,
init_func=init,
frames=200,
interval=20,
blit=False)
# start animation
[Link]()
A. Blondin Massé (UQAC) 12 juin 2014 11 / 51
Simulations physiques simples
I En combinant
I Numpy,
I Scipy,
I Matplotlib,
on peut créer en très peu de temps des simulations
physiques complexes.
I On verra quelques modèles :
I Pendule simple;
I Pendule double;
I Des particules en collision;
I Une simulation de type montagne russe.
A. Blondin Massé (UQAC) 12 juin 2014 12 / 51
Table des matières
1. Introduction
2. Mécanique des fluides
3. Pendule simple
4. Pendule double
5. Particules dans une boîte
6. Montagnes russes
A. Blondin Massé (UQAC) 12 juin 2014 13 / 51
Équations de Navier-Stokes
Équation générale :
∂v
ρ + v · ∇v = −∇p + µ∇2 v.
∂t
où
I v est le champ de vitesse du fluide;
I ρ est la masse volumique du fluide;
I p est la pression;
I µ est la viscosité.
A. Blondin Massé (UQAC) 12 juin 2014 14 / 51
Applet
Voir [Link]
fluid_water_2/.
A. Blondin Massé (UQAC) 12 juin 2014 15 / 51
Table des matières
1. Introduction
2. Mécanique des fluides
3. Pendule simple
4. Pendule double
5. Particules dans une boîte
6. Montagnes russes
A. Blondin Massé (UQAC) 12 juin 2014 16 / 51
Simulation
Voir simple_pendulum.py
A. Blondin Massé (UQAC) 12 juin 2014 17 / 51
État du système
I Une extrémité est fixée;
I L’autre est libre;
I On suppose que la ficelle est toujours pleinement tendue;
I On suppose qu’il n’y a aucune friction.
I Pour modéliser le problème, il faudrait connaître θ(t).
A. Blondin Massé (UQAC) 12 juin 2014 18 / 51
Équation différentielle (1/4)
I Avant de se lancer dans la programmation, il faut d’abord
décrire l’équation différentielle qui modélise un pendule;
I Dans un premier temps, rappelons la seconde loi de
Newton :
F = ma.
A. Blondin Massé (UQAC) 12 juin 2014 19 / 51
Équation différentielle (2/4)
I Clairement, le vecteur vitesse de la boule est toujours
perpendiculaire à la ficelle;
I Par conséquent, la composante de la force de gravité
agissant de façon perpendiculaire à la ficelle est donnée par
F = −mg sin θ = ma =⇒ a = −g sin θ.
A. Blondin Massé (UQAC) 12 juin 2014 20 / 51
Équation différentielle (3/4)
I Il nous reste à calculer l’accélération angulaire;
I Or, le déplacement est donné par s = `θ;
I En dérivant deux fois, on trouve
d2 θ
a=` .
dt2
A. Blondin Massé (UQAC) 12 juin 2014 21 / 51
Équation différentielle (4/4)
I En recombinant, on obtient l’équation différentielle :
d2 θ g
+ sin θ = 0.
dt2 `
A. Blondin Massé (UQAC) 12 juin 2014 22 / 51
Résolution de l’équation
I L’équation
d2 θ g
+ sin θ = 0.
dt2 `
est difficile à résoudre de façon exacte;
I Par contre, elle devient beaucoup plus simple si on utilise
l’approximation sin(θ) ≈ θ;
I On obtient alors une équation du premier ordre linéaire :
d2 θ g
+ θ = 0.
dt2 `
I Avec conditions initiales θ(0) = θ0 et dθ/dt(0) = 0, la
solution particulière est
r
g
θ(t) = θ0 cos t .
`
A. Blondin Massé (UQAC) 12 juin 2014 23 / 51
La classe SimplePendulum (1/2)
class SimplePendulum:
"""Simple Pendulum Class"""
def __init__(self,
init_angle=50,
L=1.0, # length of pendulum in m
G=9.8, # acceleration due to gravity, in m/s^2
origin=(0,0)):
self.L = L
self.G = G
[Link] = origin
self.time_elapsed = 0
self.init_angle = init_angle * [Link] / 180.
[Link] = self.init_angle
A. Blondin Massé (UQAC) 12 juin 2014 24 / 51
La classe SimplePendulum (2/2)
def position(self):
"""compute the current x,y position of the pendulum"""
x = [Link]([[Link][0], self.L * [Link](self.
angle)])
y = [Link]([[Link][1], -self.L * [Link](self.
angle)])
return (x, y)
def step(self, dt):
"""execute one time step of length dt and update state
"""
self.time_elapsed += dt
[Link] = self.init_angle * [Link]([Link]( self.G /
self.L) * self.time_elapsed)
A. Blondin Massé (UQAC) 12 juin 2014 25 / 51
Initialisation
#------------------------------------------------------------
# set up initial state and global variables
pendulum = SimplePendulum(40)
dt = 1./30 # 30 fps
#------------------------------------------------------------
# set up figure and animation
fig = [Link]()
ax = fig.add_subplot(111, aspect=’equal’, autoscale_on=False,
xlim=(-2, 2), ylim=(-2, 2))
[Link]()
line, = [Link]([], [], ’o-’, lw=2)
time_text = [Link](0.02, 0.95, ’’, transform=[Link])
A. Blondin Massé (UQAC) 12 juin 2014 26 / 51
Fonctions d’animation
def init():
"""initialize animation"""
line.set_data([], [])
time_text.set_text(’’)
return line, time_text
def animate(i):
"""perform animation step"""
global pendulum, dt
[Link](dt)
line.set_data(*[Link]())
time_text.set_text(’time = %.1f’ % pendulum.time_elapsed)
return line, time_text
A. Blondin Massé (UQAC) 12 juin 2014 27 / 51
Lancement de l’animation
# choose interval based on dt and the time to animate one step
from time import time
t0 = time()
animate(0)
t1 = time()
interval = 1000 * dt - (t1 - t0)
ani = [Link](fig,
animate,
frames=300,
interval=interval,
blit=False,
init_func=init)
[Link]()
A. Blondin Massé (UQAC) 12 juin 2014 28 / 51
Table des matières
1. Introduction
2. Mécanique des fluides
3. Pendule simple
4. Pendule double
5. Particules dans une boîte
6. Montagnes russes
A. Blondin Massé (UQAC) 12 juin 2014 29 / 51
Simulation
Voir double_pendulum.py
A. Blondin Massé (UQAC) 12 juin 2014 30 / 51
État du système
I À tout moment, le système est
représenté par (0, 0)
I Les angles θ1 et θ2 ;
`1
I Les vitesses angulaires θ10 et
θ20 . θ1 (x1 , y1 )
I Les coordonnées sont obtenues
par : `2
x1 = `1 sin(θ1 ), θ2 (x2 , y2 )
y1 = `1 cos(θ1 ),
x2 = `1 sin(θ1 ) + `2 sin(θ2 ),
y2 = `1 cos(θ1 ) − `2 cos(θ2 ).
A. Blondin Massé (UQAC) 12 juin 2014 31 / 51
Paramètres du système
class DoublePendulum:
"""Double Pendulum Class
init_state is [theta1, omega1, theta2, omega2] in degrees,
where theta1, omega1 is the angular position and velocity of
the first
pendulum arm, and theta2, omega2 is that of the second
pendulum arm
"""
def __init__(self,
init_state = [120, 0, -20, 0],
L1=1.0, # length of pendulum 1 in m
L2=1.0, # length of pendulum 2 in m
M1=1.0, # mass of pendulum 1 in kg
M2=1.0, # mass of pendulum 2 in kg
G=9.8, # acceleration due to gravity, in m/s^2
origin=(0, 0)):
self.init_state = [Link](init_state, dtype=’float’)
[Link] = (L1, L2, M1, M2, G)
[Link] = origin
self.time_elapsed = 0
[Link] = self.init_state * [Link] / 180.
A. Blondin Massé (UQAC) 12 juin 2014 32 / 51
Position des pendules
def position(self):
"""compute the current x,y positions of the pendulum
arms"""
(L1, L2, M1, M2, G) = [Link]
x = [Link]([[Link][0],
L1 * sin([Link][0]),
L2 * sin([Link][2])])
y = [Link]([[Link][1],
-L1 * cos([Link][0]),
-L2 * cos([Link][2])])
return (x, y)
A. Blondin Massé (UQAC) 12 juin 2014 33 / 51
Énergie du système
def energy(self):
"""compute the energy of the current state"""
(L1, L2, M1, M2, G) = [Link]
x = [Link]([L1 * sin([Link][0]),
L2 * sin([Link][2])])
y = [Link]([-L1 * cos([Link][0]),
-L2 * cos([Link][2])])
vx = [Link]([L1 * [Link][1] * cos([Link][0]),
L2 * [Link][3] * cos([Link][2])
])
vy = [Link]([L1 * [Link][1] * sin([Link][0]),
L2 * [Link][3] * sin([Link][2])
])
U = G * (M1 * y[0] + M2 * y[1])
K = 0.5 * (M1 * [Link](vx, vx) + M2 * [Link](vy, vy))
return U + K
A. Blondin Massé (UQAC) 12 juin 2014 34 / 51
Équations différentielles (1/4)
I Les équations différentielles modélisant un système de deux
pendules sont considérablement plus complexes;
I Dans un premier temps, il faut calculer le lagrangien :
L = K − P.
où
m1 θ102 `21 + m2 [θ102 `21 + θ202 `22 + 2θ10 `1 θ20 `2 cos(θ1 − θ2 )]
K =
2
P = −(m1 + m2 )g`1 cos θ1 − m2 `2 g cos θ2
I Puis les équations différentielles sont données par
d ∂L ∂L
0 − = 0, i = 1, 2.
dt ∂θi ∂θi
A. Blondin Massé (UQAC) 12 juin 2014 35 / 51
Équations différentielles (2/4)
I On obtient alors le système suivant :
θ100 = − m2 `1 θ102 sin(θ1 − θ2 ) cos(θ1 − θ2 ) + gm2 sin θ2 cos(θ1 − θ2 )
−m2 `2 θ202 sin(θ1 − θ2 ) − (m1 + m2 )g sin θ1 `1 (m1 + m2 )
m2 `1 cos2 (θ1 − θ2 )
θ200 = m2 `2 θ202 sin(θ1 − θ2 ) cos(θ1 − θ2 ) + g sin θ1 cos(θ1 − θ2 )(m1 + m2 )
−`1 θ102 sin(θ1 − θ2 )(m1 + m2 ) − g sin θ2 (m1 + m2 ) `2 (m1 + m2 )
−m2 `2 cos2 (θ1 − θ2 ) .
I Il s’agit d’équations d’ordre deux, qu’on peut convertir en
un système d’ordre un.
A. Blondin Massé (UQAC) 12 juin 2014 36 / 51
Équations différentielles (3/4)
I Considérons le changement de variables :
z1 = θ1
z2 = θ2
z3 = θ10
z4 = θ10
I Alors
z10 = θ10
z20 = θ20
z30 = θ100
z40 = θ100
A. Blondin Massé (UQAC) 12 juin 2014 37 / 51
Équations différentielles (4/4)
On obtient ENFIN un système du premier ordre :
z10 = θ10
z20 = θ20
z30 = − m2 `1 z42 sin(z1 − z2 ) + gm2 sin z2 cos(z1 − z2 )
−m2 `2 z42 sin(z1 − z2 ) − (m1 + m2 )g sin z1
`1 (m1 + m2 ) − m2 `1 cos2 (z1 − z2 )
z40 = m2 `2 z42 sin(z1 − z2 ) cos(z1 − z2 ) + g sin(z1 ) cos(z1 − z2 )(m1 + m2 )
`1 z42 sin(z1 − z2 )(m1 + m2 ) − g sin z2 (m1 + m2 )
`2 (m1 + m2 ) − m2 `2 cos2 (z1 − z2 ) .
A. Blondin Massé (UQAC) 12 juin 2014 38 / 51
Écriture des équations en code
def dstate_dt(self, state, t):
"""compute the derivative of the given state"""
(M1, M2, L1, L2, G) = [Link]
dydx = np.zeros_like(state)
dydx[0] = state[1]
dydx[2] = state[3]
cos_delta = cos(state[2] - state[0])
sin_delta = sin(state[2] - state[0])
num1 = M2 * L1 * state[1] * state[1] * sin_delta *
cos_delta + M2 * G * sin(state[2]) * cos_delta + M2
* L2 * state[3] * state[3] * sin_delta - (M1 + M2)
* G * sin(state[0])
den1 = (M1 + M2) * L1 - M2 * L1 * cos_delta * cos_delta
dydx[1] = num1 / den1
num2 = (-M2 * L2 * state[3] * state[3] * sin_delta *
cos_delta + (M1 + M2) * G * sin(state[0]) *
cos_delta - (M1 + M2) * L1 * state[1] * state[1] *
sin_delta - (M1 + M2) * G * sin(state[2]))
den2 = (L2 / L1) * den1
dydx[3] = num2 / den2
return dydx
A. Blondin Massé (UQAC) 12 juin 2014 39 / 51
Table des matières
1. Introduction
2. Mécanique des fluides
3. Pendule simple
4. Pendule double
5. Particules dans une boîte
6. Montagnes russes
A. Blondin Massé (UQAC) 12 juin 2014 40 / 51
Simulation
I Voir particle_box.py
I Chaque particule est représentée par un quadruplet
(x, y, vx , vy )
indiquant la position et la vitesse;
I Il y a deux types de collisions à traiter :
I Collisions entre particules;
I Collisions avec mur.
A. Blondin Massé (UQAC) 12 juin 2014 41 / 51
Mise à jour positions/vitesses
I Avec la numpy, la mise à jour des positions est très simple :
# update positions
[Link][:, :2] += dt * [Link][:, 2:]
I La mise à jour des vitesses tient compte de la gravité :
# add gravity
[Link][:, 3] -= self.M * self.G * dt
I Il reste à prendre en compte les collisions.
A. Blondin Massé (UQAC) 12 juin 2014 42 / 51
Gestion des collisions avec les murs
I Il suffit de vérifier si les limites ont été dépassé;
I Ensuite, on inverse le signe de la composante du vecteur
vitesse s’il y a lieu :
# check for crossing boundary
crossed_x1 = ([Link][:, 0] < [Link][0] + [Link])
crossed_x2 = ([Link][:, 0] > [Link][1] - [Link])
crossed_y1 = ([Link][:, 1] < [Link][2] + [Link])
crossed_y2 = ([Link][:, 1] > [Link][3] - [Link])
[Link][crossed_x1, 0] = [Link][0] + [Link]
[Link][crossed_x2, 0] = [Link][1] - [Link]
[Link][crossed_y1, 1] = [Link][2] + [Link]
[Link][crossed_y2, 1] = [Link][3] - [Link]
[Link][crossed_x1 | crossed_x2, 2] *= -1
[Link][crossed_y1 | crossed_y2, 3] *= -1
A. Blondin Massé (UQAC) 12 juin 2014 43 / 51
Gestion des collisions entre particules (1/2)
I Lors d’une collision élastique, les moments sont préservés :
m1 v1 + m2 v2 = m1 v10 + m2 v20 .
I De la même façon, l’énergie cinétique est préservée :
m1 v12 m2 v22 m1 v102 m2 v202
+ = + .
2 2 2 2
I Ces équations sont vérifiées pour chaque dimension (ce sont
donc des scalaires).
A. Blondin Massé (UQAC) 12 juin 2014 44 / 51
Gestion des collisions entre particules (1/2)
I Par conséquent, en combinant la conservation du moment
et de l’énergie cinétique, on en déduit
m2 (v2 − v1 )
v10 = vmoy + ,
m1 + m2
m1 (v1 − v2 )
v20 = vmoy + .
m1 + m2
A. Blondin Massé (UQAC) 12 juin 2014 45 / 51
Table des matières
1. Introduction
2. Mécanique des fluides
3. Pendule simple
4. Pendule double
5. Particules dans une boîte
6. Montagnes russes
A. Blondin Massé (UQAC) 12 juin 2014 46 / 51
Simulation
Voir roller_coaster.py
A. Blondin Massé (UQAC) 12 juin 2014 47 / 51
Courbe paramétrée
t = −1.5 t = 1.5
t = −1 t=1
t=0
I Une paramétrisation possible est
~r(t) = (t, t2 ).
I Problème : la vitesse ne doit pas dépendre de la
paramétrisation.
A. Blondin Massé (UQAC) 12 juin 2014 48 / 51
Abscisse curviligne
` = −3.1521 ` = 3.1521
` = −1.4790 ` = 1.4790
`=0
I Il suffit de reparamétrer de sorte que le paramètre
corresponde à la longueur d’arc;
I Cette paramétrisation s’appelle l’abscisse curviligne.
A. Blondin Massé (UQAC) 12 juin 2014 49 / 51
Physique des montagnes russes
I À tout moment, le vecteur vitesse est tangent à la courbe;
I Si k(t) est la pente de la courbe au point t, alors la norme
du vecteur accélération est
−gk(t)
a= p .
1 + k(t)2
A. Blondin Massé (UQAC) 12 juin 2014 50 / 51
Mise à jour de l’état de la balle
Il suffit donc, à chaque étape, de mettre à jour la position et la
vitesse comme suit :
def step(self, dt):
r"""
Updates the position and velocity according to the
given time difference
"""
self.time_elapsed += dt
# update positions
self.p += dt * self.v
# update velocity
k = [Link].p_to_k(self.p)
if k is None:
self.v += -self.G * dt
else:
self.v += -self.G * k * dt / [Link](1 + k**2)
A. Blondin Massé (UQAC) 12 juin 2014 51 / 51