David Madore's WebLog: Jouons à analyser la forme des continents

Je cherchais à me faire une idée intuitive un peu plus claire de la notion mathématique de décomposition en harmoniques sphériques (voir ici pour une explication très sommaire) : or la meilleure façon de comprendre une notion mathématique est probablement de s'amuser avec — je me suis dit que pour avoir une fonction raisonnablement « parlante » sur la sphère avec laquelle faire joujou, un candidat assez naturel est la forme des continents. J'ai donc analysé cette fonction en harmoniques sphériques ; plus exactement, j'ai pris la fonction qui vaut −1 sur la terre et +1 sur la mer, histoire d'être mieux centré vers 0, mais c'est peu important (ça va juste introduire des facteurs ½ pénibles un peu partout dans la suite), et en faisant semblant que la Terre est une sphère. Ce calcul n'a, bien sûr, rien d'original, même si le genre de fonction qu'on analyse pour des applications plus sérieuses seraient plutôt l'altitude, le champ de gravité ou quelque chose de ce goût. Je tire mes données géographiques de cette page (Earth Specular Map 8K). J'ai utilisé la bibliothèque SHTns pour faire les calculs (après une tentative pitoyable pour les faire moi-même, cf. ci-dessous).

L'image à gauche de ce texte montre les sommes partielles de cette décomposition en harmoniques sphériques : en haut, le niveau ℓ=0, en-dessous la somme des niveaux ℓ=0 et ℓ=1, puis la somme des niveaux ℓ≤2, et ainsi de suite (à chaque fois, toutes les valeurs de m, c'est-à-dire −ℓ≤m≤ℓ, sont mises pour chaque ℓ, donc si on veut, la première ligne montre 1 terme, le suivant la somme de 4 termes, puis la somme de 9 et ainsi de suite). La Terre est vue en double projection orthographique, c'est-à-dire comme si elle était vue de l'infini : hémisphère nord à gauche, hémisphère sud à droite, le pôle correspondant au centre de chaque disque, le méridien de Greenwich comme le segment horizontal reliant les pôles — tout ceci devrait être assez clair sur les dernières images où on commence vraiment à voir la forme des continents ; mais bien sûr, cette façon de projeter n'a vraiment rien à voir avec le calcul lui-même, qui est porte sur la sphère. L'image de droite montre chaque niveau d'harmoniques séparément (si on veut, chaque ligne de l'image de droite est donc la différence entre la ligne correspondante de l'image de gauche et la précédente : elle montre donc ce qui a changé ; de nouveau, à chaque fois, toutes les valeurs de m, c'est-à-dire −ℓ≤m≤ℓ, sont sommées pour le ℓ correspondant). On peut cliquer sur chacune des lignes de l'image pour la voir en plus gros. Sur l'image de gauche (sommes partielles), même si j'ai tronqué la fonction à −1 et +1, on voit assez nettement les artefacts classiques qui résultent d'une troncature de la transformée de Fourier (ici sphérique mais peu importe).

L'intérêt de cette décomposition en harmoniques sphériques est qu'elle est naturelle pour la sphère : ce que je veux dire, c'est qu'elle ne dépend pas du choix des coordonnées — de la position des pôles. Pour dire les choses autrement, si on fait tourner la sphère n'importe comment, chacun des niveaux ℓ de la décomposition (et, a fortiori, la somme des niveaux ≤ℓ) tourne de la même façon. (Il est essentiel ici de sommer tous les m : si on ne prenait que les termes avec m=0, par exemple, on obtiendrait une moyenne selon les cercles de latitude, et ça, ça dépend du choix des pôles.) Pour dire les choses encore autrement, et de façon un peu plus savante, quand on applique une rotation de la sphère, chaque harmonique sphérique Y[ℓ,m] est transformé en une combinaison linéaire des Y[ℓ,m′] pour le même ℓ (mais pour l'ensemble des −ℓ≤m′≤ℓ) : l'espace vectoriel engendré par les Y de niveau (exactement) ℓ est stable par rotations (c'est une représentation de SO(3), et c'est même, pour ceux qui savent ce que ça veut dire, la représentation irréductible de plus haut poids ℓ).

En fait, pour un algébriste, la meilleure façon de présenter les choses est certainement la suivante : l'espace vectoriel engendré par les Y de niveau ≤ℓ est tout simplement l'espace vectoriel des polynômes sur la sphère de degré ≤ℓ. (Attention cependant, comme x²+y²+z²=1 sur la sphère, le degré d'un polynôme y est mal défini ; je parle ici de l'espace, qui est de dimension (ℓ+1)², des restrictions à la sphère de l'espace — lui-même de dimension (ℓ+1)(ℓ+2)(ℓ+3)/6 — des polynômes de degré ≤ℓ en x,y,z. On peut aussi préférer utiliser les polynômes harmoniques, c'est-à-dire dont le laplacien 3D est nul : pour ceux-là, la restriction à la sphère est une bijection, le degré est bien défini et coïncide avec la graduation par ℓ.) On peut même dire mieux : si on introduit le produit scalaire défini par l'intégration sur la sphère (normalisée pour avoir surface 1), alors la composante en harmoniques de niveau ≤ℓ d'une fonction f est la projection orthogonale, pour ce produit scalaire, de f sur l'espace vectoriel des polynômes sur la sphère de degré ≤ℓ. Quant aux harmoniques sphériques réelles Y elles-mêmes, si je ne m'abuse, on peut dire que Y[0,0], Y[1,0], Y[1,1], Y[1,−1], Y[2,0], Y[2,1], Y[2,2], Y[2,−1], Y[2,−2], Y[3,0], etc. (ordonnées par ℓ puis par m en mettant les valeurs négatives après les positives), s'obtiennent par orthonormalisation de Gram-Schmidt à partir des polynômes 1, z, x, y, z², xz, x², yz, xy, z³, xz², x²z, x³, yz², xyz, x²y, etc. (ordonnés par degré total, puis par degré ≤1 en y, puis par degré en x). On obtient ainsi : Y[0,0] = 1 ; Y[1,0] = √3·z ; Y[1,1] = √3·x ; Y[1,−1] = √3·y ; Y[2,0] = √5·(z²−½x²−½y²) ; Y[2,1] = √15·xz ; Y[2,2] = √15·(½x²−½y²) ; Y[2,−1] = √15·yz ; Y[2,−2] = √15·xy ; Y[3,0] = √7·(z³−(3/2)x²z−(3/2)y²z) ; Y[3,1] = √42·(xz²−¼x³−¼xy²) ; etc.

Encore une autre façon de voir le niveau ℓ de la décomposition en harmoniques sphériques d'une fonction f est, peut-être à une constante près dont je ne suis pas très sûr, comme la convolée de cette fonction avec Y[ℓ,0] (j'insiste : convoler avec Y[ℓ,0] donne la projection sur tous les Y[ℓ,m] de ce niveau) : en général, la convolution de deux fonctions sur la sphère n'a pas de sens (on ne peut pas ajouter deux points sur la sphère), mais elle en a quand l'une des fonctions convolées est zonale, c'est-à-dire qu'elle ne dépend que de la latitude. En l'occurrence, Y[ℓ,0] vaut, à un coefficient de normalisation près, P[ℓ](cos(θ)) où P[ℓ] est un polynôme de Legendre et θ désigne la colatitude (=π/2 moins la latitude).

Du coup, les niveaux de la décomposition en harmoniques sphériques ont donc une vraie signification par rapport à la fonction sommée.

Le terme ℓ=0, ou ce que les physiciens appellent le terme monopôle, est simplement la moyenne de la fonction : dans l'exemple que j'ai pris, il nous renseigne donc sur la proportion de terre et de mer. Je trouve une moyenne de 0.4283, ce qui, compte tenu du fait que j'ai mis la terre à −1 et la mer à +1, signifie qu'il y aurait (1+0.4283)/2 soit 71.41% de mer, et 28.59% de terre ferme, sur la Terre. Je suppose que les mesures peuvent varier selon ce qu'on compte exactement comme terre et mer, notamment dans les régions polaires — je donne ici simplement ce qui résulte de l'image dont je suis partie, et je ne sais pas vraiment quelle est sa source — et peut-être quand on tient compte de l'aplatissement de la Terre, mais cette valeur est au moins réaliste. Pour dire les choses autrement, si on imagine que les terres émergées ont une densité surfacique constante égale à 1 sur la surface de la sphère (et que la mer a une densité nulle), ce qu'on mesure ici est la masse totale (c'est une façon bizarre de formuler les choses, mais la comparaison à la masse va être utile pour comprendre les deux termes suivants comme un terme de barycentre et un terme de moment d'inertie).

Le terme ℓ=1, ou terme dipôle, calcule la somme (ou la moyenne) des coordonnées x, y et z contre la fonction, donc donne aussi une information sur la Terre qui a un sens intuitif assez clair : sa direction correspond au barycentre des terres émergées, ce qui se rapporte au genre de problème dont je parlais ici. Mon calcul place ce barycentre à 44.4° de latitude (nord) et 29.0° de longitude (est), du côté de Constanța en Roumanie. Ceci colle au moins grossièrement avec ce qu'on trouve sur Wikipédia, mais celle-ci a l'air surtout de citer des crackpots qui veulent plus ou moins que ce centre ait un rapport avec la Grande Pyramide, et je ne vois pas de raison de penser que mon calcul serait moins bon que le leur (de nouveau, ça dépend sans doute surtout de ce qu'on compte comme terres émergées dans les régions arctiques).

Maintenant, il faut souligner ceci : ce dont je parle ci-dessus est la notion bien définie (en général) de barycentre sphérique, qui est tout simplement la projection sur la sphère (depuis son centre) du barycentre calculé en 3D (j'ai déjà dû citer le joli article de Galperin, A concept of the mass center of a system of material points in the constant curvature spaces, Comm. Math. Phys. 154 (1993) 63–84) ; mais dans le terme dipôle, il a bien trois composantes réelles (puisqu'il y a trois harmoniques sphériques au niveau 1, Y[1,0], Y[1,1] et Y[1,−1]), i.e., ce terme dipôle a une amplitude et pas juste une direction. Il donne donc aussi la profondeur du barycentre 3D. Mon calcul donne un moment dipolaire de la terre émergée de norme 0.0996, c'est-à-dire 34.83% du moment monopolaire (0.2859, la proportion de terre émergée, cf. ci-dessus), c'est-à-dire qu'il place le barycentre des terres émergées à 34.83% du rayon de la Terre à partir de son centre (soit à (x,y,z)=(0.2176,0.1205,0.2439) si z est orienté du centre vers le pôle nord, et x du centre vers le point de longitude 0 sur l'équateur).

(J'espère ne pas avoir mal placé un √3 ou ½ quelque part dans ce calcul : les harmoniques sphériques de niveau 1 avec la convention de normalisation que j'utilise sont Y[1,0]=√3·z, Y[1,1]=√3·x et Y[1,−1]=√3·y, du coup il y a des √3 qui se promènent ; il y a aussi un −2 à cause de ma convention sur les valeurs de la fonction, et il faut encore diviser par la valeur 0.2859 du terme monopôle si on veut obtenir la position du barycentre 3D.)

Les choses deviennent vraiment intéressantes avec le terme ℓ=2, ou terme quadrupôle, parce que, pour parler de façon un peu imagée, c'est le premier niveau de la décomposition qui a une forme : alors que Y[0,0] est juste une constante, et que toute combinaison linéaire de Y[1,m] pour m∈{0,1,−1} peut se déduire de toute autre par une rotation de la sphère et multiplication par une constante, ce n'est plus le cas pour les Y[2,m].

Comment interpréter ce terme quadrupôle ? Il calcule la somme (ou la moyenne) des monômes de degré 2 dans coordonnées x, y et z (c'est-à-dire x², xy, y², xz, etc., sachant que x²+y²+z²=1 est déjà contenu dans le terme monopôle) contre la fonction considérée. Si on reprend l'idée que les terres émergées sont une masse de densité surfacique constante placée sur la sphère, de même que le terme monopôle donne la masse totale et le terme dipôle (divisé par le terme monopôle) donne l'emplacement du centre d'inertie, le terme quadrupôle, pour sa part, donne essentiellement les moments d'inertie de cette distribution de masse, ou un peu plus précisément, la manière dont ces moments d'inertie sont orientés et diffèrent les uns des autres. Les détails sont un petit peu pénibles à calculer exactement, parce que d'une part les Y[2,m] sont des combinaisons non totalement évidentes des monômes de degré 2 (il y a notamment des √5 ou √15 qui se promènent) ; d'autre part il faut décider si on veut considérer l'inertie autour du centre de la sphère ou autour du centre d'inertie (ce qu'on fait normalement en physique), ce dernier cas revenant en quelque sorte à « retrancher » le terme dipôle ; et enfin, ce qu'on appelle la matrice d'inertie I n'est pas la même chose que la matrice des moments quadrupolaires Q (il y a une sorte de dualité : I = trace(Q) − Q, voir sur Wikipédia).

Si je suis la convention de considérer les axes par le centre de la sphère (ce qui est plus naturel pour la décomposition en harmoniques sphériques), alors je trouve les trois axes quadrupolaires suivants pour la forme des continents (i.e., terres émergées) : l'un par le point de latitude 56.0° (nord) et de longitude 69.6° (est) (par ici au nord du Kazakhstan, ou, bien sûr, son antipode) qui a un moment quadrupolaire de 0.0275 (plus 0.0953 de moment monopolaire dans chaque direction, ce qui fait 0.1228 au total), un autre par le point, perpendiculaire au précédent, de latitude 21.0° (nord) et de longitude −165.8° (ouest) (quelque part au large de Hawaï ou son antipode) pour un moment quadrupolaire de −0.0214 (plus 0.0953 égale 0.0739), et le troisième par le point, perpendiculaire aux deux précédents, de latitude 25.5° (nord) et de longitude −65.3° (ouest) (vers les Bermudes ou antipode) pour un moment quadrupolaire de −0.0061 (plus 0.0953 égale 0.0892). Peut-être des crackpots géomanciens peuvent-ils nous expliquer la signification mystique de ces points sur la surface de la Terre (les sommets de l'octaèdre quadrupolaire des continents) et nous dire s'il faut construire des pyramides à ces endroits-là pour améliorer la calibration feng shui de la planète.

(Je répète que les axes que j'ai donnés ci-dessus sont les axes d'inertie de la forme des continents par le centre de la sphère ; si on veut prendre les axes d'inertie par le centre d'inertie calculé ci-dessus, alors il y en a un qui part dans la direction 5.8° de latitude et 11.7° de longitude, un autre vers 42.1° de latitude et 107.0° de longitude, et un troisième vers 47.3° de latitude et −84.6° de longitude — mais ils ne passent pas par ces points, puisqu'ils partent du centre d'inertie. Ce calcul est probablement moins intéressant ou naturel que celui des axes quadrupolaires à partir du centre, donc je n'insiste pas : j'ai simplement mentionné les moments d'inertie comme une façon de se représenter le moment quadrupolaire.)

À partir du moment octopôle (ℓ=3), les choses deviennent vraiment difficiles à représenter de façon intuitive. Pour bien se figurer les choses, il faudrait trouver une description permettant de comprendre intuitivement les représentations correspondantes de SO(3) (ℓ=0 correspond à la représentation triviale, ℓ=1 à la représentation standard de SO(3) sur ℝ³, ℓ=2 à celle sur les matrices symétriques — ou formes quadratiques — de trace nulle, mais ensuite je n'en ai pas vraiment de présentation claire à part qu'elles sont la plus grosse composante de la puissance symétrique 𝒮^ℓ(ℝ³) et j'imagine décrite par des annulations de traces partielles) ; et éventuellement décrire les invariants de ces représentations : je ne l'ai pas fait.

Ce que je peux faire, en revanche, c'est calculer le spectre en puissance, c'est-à-dire la norme 2 de chaque niveau d'harmoniques séparément [note : il est peut-être plus courant de prendre la norme 2 au carré ; ici je prends bien la norme 2 elle-même] :

(pour la fonction qui vaut +1 sur la mer et −1 sur la terre). Les premières valeurs sont 0.4283 (la moyenne de la fonction), 0.3450 (l'amplitude du dipôle, ici égale à √3 fois celle considérée ci-dessus à cause de la normalisation, et encore ×2 à cause du choix de la fonction), puis 0.1938 (l'amplitude du quadrupôle), 0.3151 (ℓ=3), 0.3062 (ℓ=4), 0.3227 (ℓ=5), 0.1370 (ℓ=6), et 0.1798 (ℓ=7). Ces amplitudes se voient aussi assez bien sur l'image ci-dessus, à droite. Les amplitudes de ℓ=5 et dans une moindre mesure ℓ=7, ont l'air assez importantes essentiellement à cause de l'antarctique. La décroissance de la fonction ci-dessus est environ en ℓ^−0.84 (à vue de nez : je n'ai pas fait un fit précis) : si je ne me plante pas, ça veut dire que la fonction appartient aux espaces de Sobolev H^s sur la sphère pour s<0.34 environ ; et je subodore que ceci signifie que, en un certain sens, la dimension fractale des côtes des continents vaut environ 2 − 2×0.34 = 1.32 (le chiffre est, au moins, vaguement plausible), mais je n'y connais rien. S'il y a des gens qui peuvent m'éclairer, qu'ils n'hésitent pas à s'exprimer.

Parlant de ne rien y connaître, une des leçons que j'ai apprises est qu'il ne faut pas s'improviser numéricien : avant de décider d'utiliser la bibliothèque SHTns pour faire les calculs, j'ai pensé les faire moi-même. Mais d'une part j'avais la flemme de programmer une transformée de Fourier rapide (en longitude, l'analyse en harmoniques sphériques est une simple transformée de Fourier), donc mon algorithme était en complexité O(N⁴), avec N le nombre de niveaux, au lieu de O(N³). À la limite ce n'est pas très important : je peux laisser le programme tourner toute la nuit ; mais le problème était surtout que je faisais une intégrale sur la sphère (typiquement, de la fonction f multipliée par une harmonique sphérique Y dont on veut calculer l'amplitude dans la fonction) en prenant une grille régulière en latitude+longitude (celle donnée par l'image de la forme des continents) et en sommant sin(θ) fois la valeur à intégrer en chaque point de la grille (puis je divise par la normalisation, mais ce n'est pas ce qui est important). Autrement dit, je faisais une intégration par la méthode des rectangles. Le problème est que pour les harmoniques sphériques avec ℓ grand et m nul (ou en tout cas petit), Y[ℓ,m] est très concentré au pôle, et le facteur sin(θ) s'y annule, ce qui donne une fonction intégrée prenant des valeurs importantes sur un petit anneau autour du pôle, et variant rapidement dans cette région : ceci conduit à une erreur systématique très important avec la méthode des rectangles, et j'obtenais des coefficients complètement faux pour ces harmoniques (et ça se voyait nettement à la reconstitution). Je me doutais quand même que quelque chose du genre risquait de se produire, et j'avais cru contourner le problème en divisant par l'intégrale de Y² (censée valoir 1) calculée selon la même méthode, mais le problème est que les erreurs d'intégrations font que les harmoniques sphériques ne sont même plus orthogonales (il ne suffit pas de les normaliser). Moralité : il vaut mieux utiliser des bibliothèques que d'autres gens ont écrites.

Si certains veulent jouer avec les données elles-mêmes (et je vous y encourage : les mathématiques sont faites pour qu'on s'amuse avec), je les ai mises ici (trois colonnes séparées par des tabulations : ℓ, m et la valeur du coefficient devant Y[ℓ,m] ; il faut préciser qu'il s'agit des harmoniques sphériques réelles, normalisées en convenant que la mesure totale de la sphère vaut 1, et sans la « phase de Condon-Shortley »).