À l'occasion de la publication de mon article dans la revue "Physical Review E", j'entame aujourd'hui une série de billets pour retracer l'historique de ce projet et pour aborder toutes sortes de problèmes mathématiques intéressants (tout du moins, pour moi).
Tout commence un soir solitaire de Novembre 2021. Alors que je nettoie ma vaisselle distraitement, je suis absorbé par les bulles qui se forment petit à petit dans mon évier. Des sphères de tailles variables, tangentes deux à deux qui forment une structure complexe et aléatoire. J'entreprends le soir même (ou plutôt la nuit même) de concevoir un modèle simplifé pour comprendre ce phénomène et de le simuler numériquement. Il me faudra quelques jours pour sortir la tête de mes simulations et découvrir que ce phénomène a déjà un nom : "Random Apollonian Packing" (RAP) et appartient à la grande famille des problèmes d'empilement de sphères.
Toutefois, ne brûlons pas les étapes et commençons par aborder une construction qui partage beaucoup de similarités avec le modèle RAP.
Cercle d'Apollonius et théorèmes de Descartes
Le Cercle d'Apollonius désigne une famille de fractales en deux dimensions, construites sur la base de quelques règles simples :
- Soit un sous-ensemble de trois cercles tangents deux à deux
- Il existe exactement deux autres cercles, chacun tangent aux trois premiers, c'est le Théorème de Descartes
- On ajoute ces deux nouveaux cercles à l'ensemble des cercles tracés
- On obtient de nouveaux sous-ensembles de trois cercles tangents deux à deux sur lesquels on reproduit cette procédure
En choisissant trois cercles tangents deux à deux complètement arbitraires, on construit de cette façon le cercle d'Apollonius de manière récursive. Commençons par énoncer ce fameux théorème de Descartes. Celui-ci s'énonce en utilisant la courbure des cercles. La courbure d'un cercle est définie par $k = \pm 1/r$ , où $r$ est son rayon. Plus le cercle est grand, plus sa courbure est petite, et vice versa. Le signe plus dans la courbure s'utilise pour un cercle qui est tangent extérieurement aux autres cercles. Dans le cas d'un cercle tangent intérieurement on utilise le signe moins.
Théorème de Descartes
Si quatre cercles tangents entre eux ont pour courbure $k_{i}$ (pour $i = 1..4$), alors $$(k_{1}+k_{2}+k_{3}+k_{4})^{2}=2\,(k_{1}^{2}+k_{2}^{2}+k_{3}^{2}+k_{4}^{2}) $$
Si on fixe les trois courbures $k_1, k_2$ et $k_3$, cette équation est un polynôme du second degré en $k_4$ qui possède en général deux racines. On peut implémenter cette solution facilement en Python et obtenir les courbures des deux nouveaux cercles :
import numpy as np
def descartes_theorem(k1, k2, k3):
""""Theoreme de Descartes pour trois cercles tangents deux a deux
Parameters
==========
- k1: float
Courbure du premier cercle
- k2: float
Courbure du second cercle
- k3: float
Courbure du troisième cercle
Returns
=======
tuple of floats
Courbures des deux nouveaux cercles
"""
return (
k1 + k2 + k3 - 2 * np.sqrt(k1 * k2 + k2 * k3 + k3 * k1),
k1 + k2 + k3 + 2 * np.sqrt(k1 * k2 + k2 * k3 + k3 * k1),
)
Pour obtenir les centres de ces cercles, on peut interpréter les coordonnées $(x, y)$ de ces centres en termes de variables complexes $z = x + iy$ et utiliser le théorème complexe de Descartes :
Théorème complexe de Descartes
Soient quatre cercles de courbure $k_{i}$ et de centre $z_i$ (pour $i = 1..4$), alors $$(k_{1}z_{1}+k_{2}z_{2}+k_{3}z_{3}+k_{4}z_{4})^{2}=2\,(k_{1}^{2}z_{1}^{2}+k_{2}^{2}z_{2}^{2}+k_{3}^{2}z_{3}^{2}+k_{4}^{2}z_{4}^{2})$$
Si on fixe les trois cercles $(z_1, k_1), (z_2, k_2)$ et $(z_3, k_3)$, il s'agit de résoudre un polynôme complexe du second degré en $z_4 k_4$. En Python, la solution est similaire :
def complex_descartes_theorem(z1, z2, z3, k1, k2, k3):
""""Theoreme complexe de Descartes pour trois cercles tangents deux a deux
Parameters
==========
- z1: complex
Centre du premier cercle
- z2: complex
Centre du second cercle
- z3: complex
Centre du troisième cercle
- k1: float
Courbure du premier cercle
- k2: float
Courbure du second cercle
- k3: float
Courbure du troisième cercle
Returns
=======
tuple of floats
Solutions de z4*k4
"""
return (
k1 * z1
+ k2 * z2
+ k3 * z3
- 2 * np.sqrt(k1 * k2 * z1 * z2 + k2 * k3 * z2 * z3 + k3 * k1 * z3 * z1)
), (
k1 * z1
+ k2 * z2
+ k3 * z3
+ 2 * np.sqrt(k1 * k2 * z1 * z2 + k2 * k3 * z2 * z3 + k3 * k1 * z3 * z1)
)
Construction du cercle d'Apollonius
En utilisant les deux théorèmes de Descartes, il est possible de trouver deux solutions pour la courbures $k_4$ et deux solutions pour $z_4 k_4$. Ce qui fait en principe quatre combinaisons pour le couple $(z_4, k_4)$ dont deux seulement sont valables et doivent être retenues. Pour cela, on vérifie si les relations entre les coordonnées et les courbures sont compatibles $$|z_i - z_4| = \frac{1}{k_i} + \frac{1}{k_4} \quad \text{si } k_i, k_4 > 0$$
def find_next_cercle(z1, z2, z3, k1, k2, k3):
""""Trouve le cercle suivant dans la construction du cercle d'Apollonius.
Parameters
==========
- z1: complex
Centre du premier cercle
- z2: complex
Centre du second cercle
- z3: complex
Centre du troisième cercle
- k1: float
Courbure du premier cercle
- k2: float
Courbure du second cercle
- k3: float
Courbure du troisième cercle
Returns
=======
tuple of floats
Centre et courbure du prochain cercle
"""
# Deux possibilites pour k4
# On ne garde que la plus grande racine, la plus petite correspondant a un cercle
# deja construit
_, k4 = descartes_theorem(k1, k2, k3)
# Deux possibilites pour z4 * k4
z4k4m, z4k4p = complex_descartes_theorem(z1, z2, z3, k1, k2, k3)
# Donc deux possibilites pour z4
z4m, z4p = z4k4m / k4, z4k4p / k4
# score et score_2 calculent la diference entre les rayons calcules comme abs(zi - zj)
# et les rayons calcules comme 1 / ki + 1 / kj
# score_1 teste z4m
score_1 = (
np.abs(z2 - z4m) + np.abs(z3 - z4m)
)
score_1 -= 1 / k2 + 1 / k3
score_1 -= 2 / k4
# score_2 teste z4p
score_2 = (
np.abs(z2 - z4p) + np.abs(z3 - z4p)
)
score_2 -= 1 / k2 + 1 / k3
score_2 -= 2 / k4
# On retient les combinaisons valables
if abs(score_1) < abs(score_2):
return z4m, k4
else:
return z4p, k4
Parmis les deux nouveaux cercles que l'on va trouver, en général l'un d'eux sera déjà présent dans notre ensemble et l'autre sera nouveau avec une courbure plus grande. C'est pourquoi on ne garde que le cercle associé à $k_{4+}$ dans la fonction ci-dessus.
Il est maintenant facile de créer notre cercle d'Apollonius, il nous suffit d'ajouter itérativement les différents cercles et de garder à jour une liste des triplets de cercles tangents entre eux.
Pour obtenir une figure plus homogène, on aimerait dessiner les cercles dans l'ordre des rayons décroissants. Pour cela nous allons utiliser une PriorityQueue, soit une liste de priorité. Cette liste va prendre en entrée des tuples et va les classer en fonction du premier argument du tuple. Dans notre cas, ce premier argument est la courbure du prochain cercle, et les autres arguments renseignent les indices du tiplets de cercles tangents.
from queue import PriorityQueue
def add_circles(centers, curvatures, todo, cutoff):
""""Forme le cercle d'Apollonius avec tous les cercles de courbure < cutoff.
Parameters
==========
- centers: list
Liste des centres de cercles
- curvatures: list
Liste des courbures de cercles
- todo: PriorityQueue
File d'attente des triplets de cercles
- cutoff: float
Plus grande courbure consideree
Returns
=======
tuple of lists
Centres, courbures et file d'attente
"""
new_index = len(centers)
k4 = 0
while k4 < cutoff:
k4, i, j, k = todo.get()
k1, k2, k3 = curvatures[i], curvatures[j], curvatures[k]
z1, z2, z3 = centers[i], centers[j], centers[k]
z4, k4 = find_next_cercle(z1, z2, z3, k1, k2, k3)
centers.append(z4)
curvatures.append(k4)
todo.put((descartes_theorem(k1, k2, k4)[1], i, j, new_index))
todo.put((descartes_theorem(k2, k3, k4)[1], i, k, new_index))
todo.put((descartes_theorem(k3, k1, k4)[1], j, k, new_index))
new_index += 1
return centers, curvatures, todo
On initialise notre construction du cercle d'Apollonius avec trois cercles. On l'aide un peu en ajoutant le cercle circonscrit à la main.
# Conditions initiales
k1 = -1 / (1 + np.sqrt(3) / 2)
k2 = 2 / np.sqrt(3)
centers, curvatures = [], []
centers.append(0.0)
centers.append(1 + 0j)
centers.append(np.exp(2j * np.pi / 3))
centers.append(np.exp(-2j * np.pi / 3))
curvatures.append(k1)
curvatures.append(k2)
curvatures.append(k2)
curvatures.append(k2)
todo = PriorityQueue()
todo.put((descartes_theorem(k1, k2, k2)[1], 0, 1, 2))
todo.put((descartes_theorem(k1, k2, k2)[1], 0, 2, 3))
todo.put((descartes_theorem(k1, k2, k2)[1], 0, 1, 3))
todo.put((descartes_theorem(k2, k2, k2)[1], 1, 2, 3))
Et on visualise tout ça avec Matplotlib, en faisant varier la valeur de la courbure maximale !
import matplotlib.pyplot as plt
import numpy as np
plt.rc("text", usetex=True)
plt.rc("font", family="serif")
def plot_configuration(ax, centers, curvatures):
for center, curvature in zip(centers, curvatures):
circle = plt.Circle(
(center.real, center.imag),
1 / abs(curvature),
fill=False,
edgecolor="blue",
linewidth=1e-1,
)
ax.add_patch(circle)
ax.set_xlim(-2, 2)
ax.set_ylim(-2, 2)
fig = plt.figure(figsize=(10, 10))
plot_configuration(fig.add_subplot(2, 2, 1), centers, curvatures)
centers, curvatures, todo = add_circles(centers, curvatures, todo, 1e1)
plot_configuration(fig.add_subplot(2, 2, 2), centers, curvatures)
centers, curvatures, todo = add_circles(centers, curvatures, todo, 1e2)
plot_configuration(fig.add_subplot(2, 2, 3), centers, curvatures)
centers, curvatures, todo = add_circles(centers, curvatures, todo, 1e3)
plot_configuration(fig.add_subplot(2, 2, 4), centers, curvatures)
plt.show()
Propriétés statistiques
Le cercle d'Apollonius est une fractale, et on peut aisément sonder ses propriétés statistiques. On va s'intéresser en particulier à la distribution des courbures. Pour cela, on va tracer la fonction de répartition $\mathbb{P}(K \leq k)$
fig, ax = plt.subplots()
centers, curvatures, todo = add_circles(centers, curvatures, todo, 1e5)
cumulant = np.array(range(1, len(curvatures) + 1)) / len(curvatures)
k_range = np.array(sorted(curvatures))
ax.loglog(k_range, cumulant)
ax.set_xlabel(r"$k$")
ax.set_ylabel(r"$P(K \leq k)$")
plt.show()
On obtient une belle loi de puissance à large courbure, qui est de plus en plus lisse à mesure qu'on augmente le nombre de cercles dans notre fractale.
Une régression linéaire permet de trouver l'exposant de cette loi de puissance. Pour cela utilisons la fonction polyfit de numpy
lower, upper = 1e4, 6e4
m, c = np.polyfit(
np.log(k_range[(k_range > lower) & (k_range < upper)]),
np.log(cumulant[(k_range > lower) & (k_range < upper)]),
1,
)
print(m)
1.3055483712995455
Ce nombre est une approximation de la dimension de Hausdorff du cercle d'Apollonius, dont les déterminations les plus récentes donnent approximativement $1.3057$ (voir la suite OEIS correspondante). On n'est pas loin !
Conclusion et la suite
On a vu dans ce billet comment produire la fractale du cercle d'Apollonius en utilisant quelques théorèmes simples de géométrie, et comment optimiser la construction de cette fractale avec une liste de priorités. On a aussi vu que la distribution des rayons dans cette fractale prend la forme d'une loi de puissance avec un certain exposant lié à la dimension fractale.
Dans un prochain billet, on va découvrir ensemble la version aléatoire de cette figure, le Random Apollonian Packing et apprendre à la construire de manière relativement optimisée.
Sur ce, je vous laisse avec une belle image vectorisée du cercle d'Apollonius !