PROJET PHYSIQUE 330

Projet A: Évacuation d’une salle & déplacement d’une foule dans une rue
Le comportement d’une foule est un problème aux applications multiples: évacuation d’une salle,
couloir du métro aux heures de pointes, manifestations… On peut en imaginer des modèles simples.
P. ex., on peut décrire chaque individu par sa position, sa vitesse, et comme étant soumis à des «
forces »:
• Une force qui spécifie la direction dans laquelle l’individu veut se déplacer,
fdir=(v0−v(t))/τ
, où v0 est la direction et la vitesse que la personne veut atteindre, v sa vitesse actuelle, et τ
• un temps caractéristique d’ajustement.
• Une force qui l’oblige à éviter des obstacles qui peuvent être fixes (un mur, un massif de
fleurs, …), ou qui peuvent être les autres individus eux-mêmes. On pourra essayer
fobs(d)=aexp(−d/d0)
, où d est la distance entre le piéton et l’obstacle, d0 la « portée » de la force, et a son amplitude.
On pourra varier les différents paramètres apparaissant ci-dessus, tout en leur donnant une
interprétation physique réelle, et étudier leur influence dans des situations concrètes. P. ex., à quelle
vitesse, en fonction de v0 et de la densité de piétons, se déplace une foule contrainte à avancer dans
un couloir si chaque individu veut maintenir une vitesse v0 ?
Comment s’organise l’évacuation d’une salle initialement uniformément peuplée, avec une ou
plusieurs sorties, et en la présence éventuels d’obstacles ?
Il est également possible d’essayer d’autres expressions pour les forces.
Il existe une littérature conséquente sur le sujet, que l’on pourra explorer si besoin

Fig. 13.6 Figure: Un obstacle aide à l’évacuation (DR).

Projet B: États de diffusion pour l’équation de Schrödinger 1D stationnaire
On s’intéresse à la diffusion d’une particule de masse m à travers un potentiel carré défini par
V(x)=V0 pour 0≤x≤a , et 0 sinon.
Les solutions de cette équation en dehors de la région où règne le potentiel sont connues. Les
paramètres d’intégration de ces fonctions d’onde peuvent se déterminer par les relations de
continuité aux frontières avec la région où règne le potentiel. En résolvant l’équation différentielle
dans la région du potentiel pour x allant de a à 0 on peut obtenir une autre valeur pour ces paramètre
d’intégration. Il faut ensuite appliquer un algorithme de minimisation pour déterminer les
constantes d’intégration.
Les objectifs de ce projet sont :

• Écrire un programme qui résolve l’équation de Schrödinger.

• Représenter les fonctions d’ondes des trois régions avec la barrière de potentiel
• En déduire et représenter les fonctions de transmission et de réflexion en fonction de la
largeur de la barrière puis en fonction de l'Énergie

Projet C: Solitons
On considère un câble sous tension auquel sont rigidement et régulièrement attachés des pendules.
Les pendules sont couplés grâce au câble à travers sa constante de torsion. Dans un tel système on
peut observer une large gamme de phénomènes ondulatoires. Le but de cet projet est d’étudier une
solution très particulière : le soliton.
Imaginons qu’une des extrémités du câble est attachée à une manivelle qui peut tourner librement. Il
est alors possible de donner une impulsion au système en faisant un tour rapide ce qui déclenche la
propagation d’un soliton. Dans ce projet, on considérera les pendules individuellement. Il n’est pas
demandé de passer au modèle continu et de résoudre l’équation obtenue.
Pour chaque pendule n dont la position est décrite par θn
, l’équation d’évolution s’écrit:

sont des paramètres physiques. On résoudra numériquement cette équation pour chaque pendule. En
donnant un « tour de manivelle numérique », on essayera d’obtenir la solution soliton. On cherchera
en particulier à ajuster la solution par une équation du type

où a,b,n0
sont des paramètres à déterminer.
De très nombreuses questions se posent (il ne vous est pas demandé de répondre à chacune d’entre
elle):
 • Est-il toujours possible d’obtenir un soliton?
• Sa vitesse est-elle constante?
• Le soliton conserve-t-il sa forme?
• Que se passe-t-il avec des pendules plus lourds? ou plus rapprochés? avec un câble plus
rigide? avec un frottement?
• Comment le soliton se réfléchit-il si l’extrémité du câble est rigidement fixée? et si elle
tourne librement?
• Dans ce système, le soliton est chiral. En effet, on peut tourner la manivelle à gauche ou à
droite. Un anti-soliton a-t-il les mêmes propriétés (taille, vitesse, énergie) qu’un soliton?
• Si on place une manivelle à chaque extrémité, on peut faire se collisionner des solitons.
Cette étude est très intéressante et pleine de surprises. Que se passe-t-il lors de la collision de
deux solitons? Entre un soliton et un anti-soliton?

Projet D. Modélisation de la chute libre

Nous lâchons une balle d'une hauteur donnée. Une fois que la balle a quitté la main, la seule force
qui s'applique sur elle est la force de gravité (nous négligeons les frottements).
On peut alors écrire le principe fondamental de la dynamique :

En projettant sur l'axe vertical, on obtient :

La masse de l'objet n'intervient pas dans la vitesse de chute.

Méthode d'Euler :
Nous allons maintenant calculer à chaque instant la position et la vitesse de la balle à l'aide d'une
méthode numérique de résolution de l'équation différentielle basée sur le schéma d'Euler, décrit ci-
dessous.
À partir de la connaissance d'une valeur de la vitesse à un instant (v n à l'instant tn), on calcule la
valeur suivante (vn+1 à l'instant tn+1) à l'aide de l'approximation de la dérivée :
Principe de la méthode d'Euler. Le cas présenté correspond à la méthode d'Euler explicite : le calcul
se fait par approximation sur la tangente au point initial (v n).

Un choix simple est de prendre un intervalle de temps régulier, c'est à dire de fixer la
valeur de tn+1 - tn. Mais on pourrait tout-à-fait choisir un pas de temps variable en
fonction de l'optimisation et de la précision du calcul pour des moments choisis.
Nous présentons ici la méthode d'Euler, il existe de nombreux autres algorithmes
d'intégration plus ou moins compliqués à mettre en œuvre : Runge Kutta...
On remplace donc la dérivée par la valeur de la pente. Une telle approximation n'est raisonnable
que si la variation entre les deux instants est toute petite.
Finalement, il suffit de répéter l'itération afin de connaitre à chaque instant la valeur de la vitesse le
long de la chute :
• Ecrire un programme python pour  écrire la chute de la balle en représentant la position et
la vitesse de la balle en fonction du temps sur deux courbes séparées

•  Faire une résolution analytique et représenter les solutions sur la même courbe que cette
des solutions numérique

• Utiliser la fonction animation() de python pour montrer la dynamique de chute de la balle.

Représenter la balle comme “*” et le sol comme une droite.
 Projet PHY330: 2021–2022

Projet E

Modélisation d’une épidémie avec le modèle SIR

1 Présentation du Sujet
Le modèle SIR représente un modèle simpliste, qui permet la modélisation du cours d’une épidémie dans une
population. Créé en 1927 par A.G. McKendrick et W.O. Kermack [1], ce modèle épidémiologique est capable
de décrire qualitativement l’évolution du nombre des individus infectieux dans une population en fonction du
temps : Dans le modèle SIR, un individu commence sain, peut devenir infecté puis se remettre de sa maladie
avec une immunisation. Le modèle peut être facilement adapté pour simuler l’effet des mesures d’intervention
([Link]. un confinement) et avec un paramétrage sophistiqué on peut même faire des prédictions réalistes du taux
d’infection pour l’épidémie Covid-19 actuelle.

F IGURE 1 – (A) Schéma du modèle SIR. (B) Nouveaux cas Covid-19 confirmés par jour au début 2020 en
Allemagne et au Royaume-Uni. (C) Dynamique typique du modèle SIR.

2 Modèle et Méthode
La population est divisée en trois compartiments qui correspondent à différents états vis-à-vis d’une ma-
ladie : Les individus susceptibles à la maladie S, les individus infectieux I et les individus éliminés R (guéris
et donc immunisés, ainsi que les individus morts). Le nombre d’individus dans chaque compartiment varie en
fonction du temps t, on les désigne par S(t), I(t) et R(t).
Pour spécifier les taux d’individus passant pendant un incrément de temps dt entre les compartiments, on
définit un ensemble de règles : Le nombre d’individus sains est diminué par les nouvelles infections, ce qui
dépend du nombre d’individus pouvant transmettre la maladie (I) et du nombre d’individus disponibles pour
l’infection (S). Nous désignons le paramètre de transmission par β , de sorte que S varie comme suit

dS(t)
= −β .I(t).S(t) (R1)
dt
Parmi les individus infectés, une partie des individus se retire, soit parce qu’ils ont récupéré de la maladie et

1
sont devenus immunisés, soit parce qu’ils sont morts. Ce taux de transition de I à R est noté par γ, soit
dI(t)
= β .I(t).S(t) −γ.I(t) (R2)
dt
Par conséquence le nombre d’individus retirés augmente en fonction du temps comme
dR(t)
= γ.I(t) (R3)
dt
Les règles (R1)-(R3) forment un système d’équations différentielles, qu’on peut résoudre numériquement.

3 Objectif
Dans la première étape du projet, il est nécessaire de coder et analyser l’évolution des compartiments en
fonction du temps selon le modèle SIR : Quel paramétrage mène à une infection de toute la population ? Dans
quelles conditions la maladie continue-t-elle de gagner du terrain ? Quel est le nombre de reproduction de base
R0 et le nombre de reproduction à l’instant t, R(t), dans les deux cas ? [2]
Dans une seconde partie, on pourra modifier et étendre le modèle pour le rendre plus réaliste.

4 Études possibles d’intérêt

— Taux de transition dépendants du temps : Pour simuler un changement du virus ou des circonstances
concomitantes, on pourra modifier les taux de transition en fonction du temps. Cela permettra d’étudier,
par exemple, l’effet d’un confinement ou d’une mutation du virus sur l’évolution de la pandémie.
— Ajout d’autres compartiments : Le modèle SIR est seulement le modèle de base. On pourra ajouter
plusieurs aspects en introduisant des compartiments supplémentaires : un compartiment V pour simu-
ler l’effet de l’immunité protectrice par la vaccination ; un compartiment Q pour étudier l’effet d’une
quarantaine ; un compartiment E pour les individus exposés (infecté, mais pas encore infectieux) ; etc.
— Ajustement des paramètres : Même si le modèle SIR simplifie fortement la dynamique infectieuse
d’une pandémie, on pourra essayer d’ajuster ses paramètres pour décrire aussi fidèlement que possible
une épidémie réelle. À cet égard, il est essentiel de simuler une épidémie bien documentée et pour
laquelle il existe suffisamment de données fiables (Covid-19 ; la peste à Madagascar 2017 ; etc.).
— Fluctuations statistiques - l’algorithme de Gillespie : On pourra en fait considérer les taux de transition
entre les différents compartiments comme différentes réactions chimiques possibles avec des taux de
réaction donnés. Pour un tel système d’équations stochastiques, D. Gillespie a développé en 1976 un
algorithme Monte-Carlo : On tire deux nombres aléatoires pour calculer le temps τ jusqu’à la prochaine
réaction et pour déterminer quelle réaction a effectivement lieu au temps t + τ, voir la référence [3].
— Discrétisation du modèle sur réseau : Si l’on veut simuler la propagation d’un virus à Paris, on pourra
représenter chaque arrondissement par un site du réseau et ne relier que les arrondissements voisins
par un lien. En définissant alors les compartiments du modèle SIR pour chaque site et en ajoutant une
mobilité des individus entre les sites, nous obtenons une version discrétisée du modèle SIR. Pour un
exemple d’une telle discrétisation pour une épidémie de zombies aux Etats-Unis, voir la référence [4].

Références
[1] A.G. McKendrick et W.O. Kermack, "A Contribution to the Mathematical Theory of Epidemics." Proc. Roy. Soc. A
115, 700–721 (1927).
[2] N. Bacaër, "Un modèle mathématique des débuts de l’épidémie de coronavirus en France." Mathematical Modelling
of Natural Phenomena 15, 29 (2020).
[3] D.T. Gillespie, A. Hellander et L.R. Petzold, "Perspective : Stochastic algorithms for chemical kinetics." J. Chem.
Phys. 138, 170901 (2013).
[4] A.A. Alemie, M. Bierbaum, C.R. Myers et J.P. Sethna, "You can run, you can hide : The epidemiology and statistical
mechanics of zombies." Phys. Rev. E 92, 052801 (2015).

2
LU3PY124 PROJET F SU

Projet: Dynamique Quantique

La possibilité de décrire l’évolution temporelle d’un système quantique, c’est à dire

l’évolution au cours du temps de sa fonction d’onde Ψ(�r, t), est un des résultats majeurs
de la physique quantique.
Pour cela, il est nécessaire de résoudre l’équation de Schrödinger dépendante du temps
∂
ih̄ Ψ(�r, t) = Ĥ(t)Ψ(�r, t) (1)
∂t
où Ĥ(t) est le hamiltonien du système et où l’on suppose connue la fonction d’onde Ψ(�r, t =
0) à un instant initial t = 0. Hormis quelques cas particuliers où des solutions analytiques
peuvent êtres trouvées (souvent au prix d’approximations), on a recourt en général à des
méthodes d’intégration numériques.

Figure 1: Exemple d’évolution temporelle d’un système quantique, en l’occurrence

l’observation des oscillations de Bloch à l’aide d’un "simulateur quantique" constitué par 5
qubits supraconducteurs. Première ligne : résultats expérimentaux. Deuxième ligne simula-
tions numériques. Figure extraite de l’article Guo, XY., Ge, ZY., Li, H. et al. Observation
of Bloch oscillations and Wannier-Stark localization on a superconducting quantum proces-
sor. npj Quantum Inf 7, 51 (2021).

Dans ce projet, on mettra en œuvre ces méthodes numériques dans le cas d’une particule
de masse m soumis à l’action d’un potentiel unidimensionnel (1D) indépendant du temps
V (x) : � �
∂ h̄2 ∂ 2
ih̄ Ψ(x, t) = − + V (x) Ψ(x, t). (2)
∂t 2m ∂x2

1/4
LU3PY124 DQ SU

Par commodité, on prendra un système d’unité dans lequel l’équation (2) peut se mettre
sous la forme: � �
∂ 1 ∂2
i Ψ(x, t) = − + V (x) Ψ(x, t) (3)
∂t 2 ∂x2

Pour pouvoir effectuer une résolution numérique du problème, on restreindra les valeurs
de x à un intervalle de largeur L autour de l’origine. Cela revient à exclure la présence de
la particule pour les valeurs de |x| > L/2 (on suppose donc V (x) = ∞ pour |x| > L/2).
De même, l’évolution temporelle sera limitée à une durée T à déterminer en fonction des
paramètres du problème.
De plus, on effectuera une discrétisation du problème à la fois pour les positions x et
pour le temps t, c’est-à-dire que l’on fait l’approximation de remplacer la variable continue
x par N variables discrètes {xi } = {x0 , x1 , . . . xN −1 } espacées de δx et la variable continue t
par M variables discrètes tn = {t0 , t1 , . . . tM −1 } espacés de δt. Ceci permettra de remplacer
les dérivées partielles par rapport à x et t par des différences finies entre valeurs de Ψ(xi , tn ).

1 États stationnaires
Dans un premier temps, on pourra considérer uniquement les états stationnaires du prob-
lème, ce qui permettra de traiter d’abord la discrétisation spatiale. Il s’agit alors de
déterminer numériquement les énergies propres Ek et les fonctions d’ondes φ( k) associées,
obtenues comme solutions de
Ĥφk (x) = Ek φk (x) (4)
.
Ceci revient à calculer les valeurs propres et les vecteurs propres associés d’une matrice
H de dimension (N × N ) finie, représentant Ĥ.
Du fait de la discrétisation spatiale, on pourra considérer la fonction d’onde comme un
vecteur composé des N valeurs ψi = ψ(xi ) (de même, on peut noter Vi = V (xi ). Pour
construire H, on fera l’approximation d’ordre 2 en δx,
∂ 2 ψ(x) ψi+1 + ψi−1 − 2ψi
2
|x=xi ≈ (5)
∂x δx2
qui conduit à une matrice tridiagonale. Il ne sera pas forcément nécessaire d’utiliser une
écriture matricielle pour H, les termes non-diagonaux étant constants et égaux.
Pour trouver les solutions de l’équation aux valeurs propres obtenues, on pourra se
servir de bibliothèques déjà existantes et testées. Pour vérifier les solutions obtenues, on
les comparera avec les solutions analytiques dans les cas particuliers du puits de potentiel
infini V (x) = 0 pour −L/2 < x < L/2 et du potentiel harmonique V (x) = ω 2 x̂2 /2. Ceci
permettra également de discuter l’effet de la discrétisation spatiale et des conditions aux
limites x = ∓L/2 sur le spectre et les fonctions d’ondes des états stationnaires.

2/4
LU3PY124 DQ SU

On pourra considérer alors différents potentiels, par exemple :

1. Une barrière de potentiel de largeur d : V (x) = U0 si |x| < d/2 et V (x) = 0 sinon
(avec bien sûr d inférieur à L).

2. Une barrière de potentiel, moins abrupte, de largeur d : V (x) = U0 cos2 (πx/d) si

|x| < π/2 et V (x) = 0 sinon.

3. Un double puits de potentiel symétrique : V (x) = ax4 − x2 /2 avec un maximum local

en x = 0 et deux minima symétriques de part et d’autre.

4. Un potentiel périodique : V (x) = U0 cos(2πx/d).

5. . . .

2 Évolution temporelle dans le cas d’un puits infini

et d’un potentiel harmonique
On étudie dans cette partie l’évolution temporelle de la fonction d’onde Ψ(x, t) d’un état
décrit par une fonction d’onde initiale Ψ(x, 0). On considèrera d’abord le puits infini et le
potentiel harmonique, pour lesquelles les solutions analytiques exactes sont connues, afin
de vérifier l’évolution temporelle obtenue numériquement.
Cette évolution est obtenue à partir de l’équation 3. Comme dans le cas stationnaire,
on effectuera une discrétisation spatiale. Pour obtenir l’évolution temporelle, on sera égale-
ment amené à effectuer une discrétisation temporelle. Soit Ψn et Ψn+1 les vecteurs de di-
mension N de composantes Ψn,i et Ψn+1,i , représentant respectivement les N composantes
de la fonction Ψ discrétisée spatialement aux temps tn et tn+1 .
Nous utiliserons dans un premier temps l’algorithme d’Euler "explicite" pour obtenir
Ψn+1 en fonction de Ψn . Cet algorithme du premier ordre en δt consiste à faire l’approximation

∂ Ψ(x, tn+1 ) − Ψ(x, tn )

Ψ(x, t)|t=tn ≈ (6)
∂t δt
et conduit à écrire que
Ψn+1 = (1 − iHδt)Ψn (7)
Nous illustrerons sur ces potentiels simples les limites de cet algorithme (précision,
capacité à assurer la conservation de la norme de la fonction d’onde ...). Puis nous mettrons
en œuvre un algorithme de Runge-Kutta d’ordre 4, d’ordre plus élevé en δt et donnant une
estimation plus précise de Ψn+1 que celle donnée par la méthode d’Euler. Pour aller au delà
nous étudierons, dans la mesure du temps disponible, un algorithme permettant d’assurer
une meilleure conservation de la norme de la fonction d’onde.

3/4
LU3PY124 DQ SU

Pour vérifier l’évolution temporelle obtenue, on pourra représenter au cours de temps

la densité de probabilité de présence de la particule ainsi que tracer l’évolution de la norme
de Ψ aux différents temps tn et les valeurs moyennes des opérateurs position, impulsion,
du Hamiltonien...
On pourra partir d’un état initial qui est un état stationnaire et vérifier qu’il n’évolue pas
au cours du temps. On pourra ensuite considérer des superpositions d’états stationnaires
et observer leurs évolutions.
On pourra également considérer un état initial décrit par une fonction d’onde initiale
Ψ(x, 0) qui sera un paquet d’onde gaussien de largeur σ centré à la position x0 . Dans le
cas d’un puits infini, quel phénomène observe-t-on ? Comment varie-t-il en fonction de la
largeur σ ? Dans le cas d’un potentiel harmonique, qu’observe-t-on pour x0 = 0? Et pour
x0 �= 0. Pourquoi peut-on appeler ces états "quasi-classiques" ?

3 Évolution temporelle : oscillations de Bloch

Uen fois la résolution numérique validée, on peut s’intéresser à l’évolution temporelle de
la fonction d’onde pour différents potentiels. A vous d’en imaginer d’intéressants. Vous
pouvez considérer les potentiels décrits dans la partie 1 ou d’autres.
Par exemple, une possibilité est de simuler un phénomène surprenant de la physique
quantique et la physique du solide : les "oscillations de Bloch". Cet effet se produit lorsqu’on
place une particule dans un potentiel périodique et qu’elle est de plus soumis à une force
externe constante F :
V (x) = U0 cos2 (πx/d) − F x (8)
C’est le cas par exemple des électrons dans un solide soumis à un champ électrique statique.
Contrairement à ce que notre intuition "classique" pourrait prédire, la particule ne va pas
"descendre" le long du potentiel mais va au contraire osciller avec une période et une
amplitude inversement proportionnelle à la force appliquée. Cet effet prédit par Félix
Bloch et Clarence Zener au début des années 30 n’a pu être observé qu’une soixantaine
d’années plus tard dans différents systèmes. Il est notamment utiliser pour effectuer des
transferts d’impulsions contrôlés à des atomes neutres.

4/4
 PROJET G

Bille rebondissant sur une plaque en oscillation verticale

November 23, 2021

En 1982, P.J. Holmes a été le premier à formaliser une étude sur le comportement d’une bille rebondissant verticalement
sur une plaque en oscillation sinusoïdale [1]. L’objectif de ce projet est de modéliser et simuler les rebonds d’une telle
bille sur un plateau horizontal en oscillation verticale à une pulsation ω (voir Figure 1).

Figure 1: Une bille rebondissant sur un plateau horizontal rigide oscillant à une pulsation ω.

Ce système unidimensionnel offre une phénoménologie étonnamment riche et complexe. D’une part, lorsque la pulsation
du plateau est basse, il existe des solutions périodiques telles que la bille rebondit n fois en une période du plateau.
D’autre part, en augmentant la pulsation du plateau, on observe une transition vers un régime chaotique, dans lequel le
mouvement de la bille ne présente plus aucune régularité. Cette transition a lieu via une suite de doublement de période,
un phénomène observable numériquement en traçant le diagramme de bifurcation du système. Dans ce régime, le
système présente différentes caractéristiques emblématiques des systèmes chaotiques : sensibilité aux conditions initiales,
convergence vers un attracteur étrange (voir Figure 2), exposant de Lyapunov, hystérésis, etc.

Figure 2: Variation de l’attracteur étrange de Lorenz.

Modélisation du choc La collision entre le plateau et la bille est modélisée par un choc amorti, avec un coefficient
d’amortissement µ. Si la bille possède une vitesse v 0 et le plateau une vitesse vp0 lors du choc, alors l’effet d’amortissement
se traduit par une modification de la vitesse de la bille v 1 après le rebond :

v 1 − vp0 = −µ(v 0 − vp0 )

Régime colle et prochain décollage Il se peut que la bille reste collée sur le plateau pendant un certain temps, si sa
vitesse après le rebond n’est pas supérieure à la vitesse du plateau. Dans ce cas de figure, les chocs ont lieu de manière
très rapprochée. En deçà de 10−2 s entre deux chocs, on considérera que la bille est dans le régime colle. La bille est alors
portée par le plateau pendant une certaine durée. Déterminer une condition sur l’accélération du plateau pour que la
bille en décolle.

1 Pistes d’investigation
1. Faire varier ω pour observer un régime périodique, puis un régime chaotique.
2. Tracer le diagramme de bifurcation. Déterminer la période des régimes périodiques.
3. Tracer les trajectoires de deux billes pour des conditions initiales très proches (v 0 = 0.2, 0.2 + 10−5 ), observer la
divergence des trajectoires, calculer l’exposant de Lyapunov.

4. Modifier progressivement la pulsation du plateau pour observer un phénomène d’hystérésis.

5. Tracer le temps entre deux chocs τn en fonction de la vitesse de la bille vn après le nième choc. Observer l’attracteur
étrange dans le régime chaotique.
6. Modifier le code pour simuler deux billes attachées par un fil sans masse de longueur l, rebondissant sur le même
plateau en oscillation.

References
[1] P.J. Holmes. The dynamics of repeated impacts with a sinusoidally vibrating table. Journal of Sound and Vibration,
84(2):173–189, 1982.