Dans son article sur la méthode de Newton, Tan Lei a expliqué comment cette méthode engendre des dessins fractals.

Voici à gauche un dessin fractal typique, généré par l’application de la méthode de Newton avec $P(z)=z^3-1$, donc la recherche des trois racines de ce polynôme en nombres complexes. Tous les points des zones d’une certaine couleur, lorsqu’ils sont pris comme point de départ pour l’itération, mènent à la même racine : par exemple, tous les points dans les zones rouges mènent à la racine $z=1$.

Il y a un nombre infini de zones rouges, vertes et jaunes, et les frontières entre eux sont fractales.

Si nous pouvons imposer un certain relief aux différentes zones, on aura sans doute des images plus intéressantes. On peut ensuite emprunter des techniques de lancer de rayon (raytracing) pour dessiner des paysages fractals « newtoniens ».

Le lancer de rayon sur Newton

On va s’arranger pour que les bords fractals entre les zones restent dans un plan, mais que tout le reste soit plus bas, selon une certaine fonction dont on reparlera plus tard.
On va donc créer un paysage avec un nombre infini de vallées et de crêtes.

On va se mettre sur un point de vue au dessus des crêtes, et tirer des rayons vers le paysage pour trouver les intersections de ces rayons avec le « sol ». Il est impossible de trouver ces points d’intersection analytiquement : le ’sol’ est une surface fractale, qu’on ne peut pas exprimer avec une formule $w=f(x,y)$. [1] Pour chaque point sur cette surface, on doit effectuer une itération pour en trouver les coordonnées. On a donc besoin d’une méthode algorithmique pour trouver l’intersection avec un rayon : faire avancer un point sur ce rayon à petits pas. Quand on aura finalement trouvé ce point, il ne reste alors qu’à afficher à l’écran la couleur du « sol » en ce point. Cette couleur sera déterminée par la zone (toutes les zones qui mènent à la même racine ont la même couleur), la position d’une source lumineuse et éventuellement des ombres jetées par des crêtes...

Où se trouve alors le « sol » ? Pour commencer, définissons le sol comme une simple fonction du nombre d’itérations $u$ qui fait converger la méthode de Newton pour le point en question. Ce point a des coordonnées $(x,y)$ dans le plan et une coordonnée verticale $w$. On fait l’itération de Newton $z_{n+1}=z_{n}-\frac{P(z)}{P'(z)}$ avec le point complexe $z=x+iy$, on s’arrête quand $|z_{n+1}-z_{n}|<\epsilon$, où $\epsilon$ est une petite valeur, typiquement $10^{-4}$ et on trouve un nombre d’itérations $u$. Pour créer nos vallées, définissons alors la profondeur de la vallée en ce point comme, par exemple, $w=-k/u^3$, avec $k$ une constante de notre choix.

Pour déterminer les rayons à tirer, on prend un plan qui sera le plan de projection, donc l’écran de l’ordinateur. C’est un plan perpendiculaire par rapport à un vecteur (en orange dans la figure) du point de vue ($PDV$) vers un certain point à son choix dans le plan horizontal (Dans la figure on a pris l’origine du plan horizontal comme point de mire). On tire un rayon vers chaque point (pixel) de l’écran, et on obtiendra donc une projection de $ABCD$ vers $A'B'C'D'$. On va faire avancer un point sur chaque rayon et s’arrêter quand le point rencontre le « sol ».

Comme le relief du paysage se trouve sous le plan horizontal, on peut prendre comme point de départ l’intersection du rayon avec ce plan.

On avance notre point sur le rayon à petits pas qui sont tous de la même longueur $d$ et, à chaque pas on fait l’itération : on connaît les coordonnées $(x_p, y_p, w_p)$ du point et on fait l’itération avec $z=x_p+iy_p$. Si la position verticale $w_p$ du point sur le rayon est plus haute que $w=-k/u^3$, c’est qu’on n’a pas encore rencontré le sol et on peut continuer à avancer. Après un certain nombre de pas, la position verticale sera en effet inférieure à $w$ et on saura alors qu’on est allés trop loin ! Le point est enterré sous le sol.

Sur la figure à gauche, le point $B$ est la position sur le rayon qui a franchi la surface de la vallée (en bleu) et le point précédent, c’était $A$. Pour trouver l’intersection exacte avec la fonction qui décrit la vallée, on démarre maintenant une recherche par dichotomie.

Le point $M_1$ est le point à mi-chemin entre $A$ et $B$. Si $M_1$ est au-dessus de la vallée, comme c’est le cas sur la figure, on oublie le point $A$ et on calcule le point $M_2$, à mi-chemin entre $B$ et $M_1$. (Si par contre $M_1$ est en dessous, on oublie le point $B$ et on calcule le point à mi-chemin entre $A$ et $M_1$.) Si $M_2$ est en dessous, on calcule le point à mi-chemin entre $M_1$ et $M_2$, etc. En continuant selon ce principe, on trouve assez rapidement le point exact sur le rayon qui est l’intersection avec la fonction de la vallée, avec une précision de notre choix : il suffit d’imposer la fin de l’itération quand la distance entre des points consécutifs devient inférieure à une limite de notre choix (par exemple $10^{-8}$).

Mais il y un problème. Il est clair sur la figure à droite que notre méthode ne marche plus quand les rayons coupent les crêtes. Les points consécutifs $A$ et $B$ sont tous les deux au-dessus de la courbe bleue et si on ne fait rien, le point suivant sera encore plus loin sur le rayon que $B$. On aura manqué la crête !

On pourrait avancer à très petits pas (diminuer considérablement la distance $d$ entre les points $A$ et $B$) pour avoir plus de chance de ne rien manquer, mais cela « coûte » en temps de calcul, et on n’est jamais sûr de n’avoir rien manqué.

Pour améliorer la méthode à pas de longueur constante, il nous faut une estimation de la distance horizontale entre le point sur le rayon et les crêtes, parce que cela permet d’adapter les pas : des grands pas quand on est loin, de tous petits pas quand on est proche d’une crête. La longueur des pas $d$ doit donc devenir variable.

Pour des ensembles fractals comme celui de Mandelbrot, il y a des formules qui donnent une estimation de la distance entre un point et le bord de l’ensemble, mais pour les fractales à la Newton ces formules sont un peu plus difficiles à trouver.

Heureusement il y a bien une formule empirique [2] qui marche très bien.

Sur la figure à gauche on voit comment l’estimation de distance ($ED$) permet de s’approcher très prudemment de la crête. La longueur des pas diminue quand on s’approche, et on ne va pas passer les pics.

Avec cette formule, on a tous les outils nécessaires. Dès qu’on connaît la position d’un point sur l’intersection, il nous faut aussi la normale à la surface. On peut obtenir celle-ci en tirant quelques autres rayons qui donnent d’autres points d’intersection, et alors calculer une approximation pour la normale par un produit vectoriel. L’angle de la normale avec un vecteur vers la source lumineuse va déterminer l’intensité de la couleur à ce point. (voir ombrage Phong )

Si on veut montrer des ombres, on tire un rayon dès la source lumineuse vers le point d’intersection. Si ce rayon arrive à un autre point d’intersection, alors notre point et dans l’ombre et on adapte la couleur.

Tout cela veut dire qu’on a un très grand nombre de calculs à faire pour la réalisation des images. Toutes les images dans cet article viennent du logiciel Ultrafractal. Une image de petite taille prend moins de 30 secondes, mais pour des images de grande taille ( par exemple $1600*1000$ pixels), avec anticrénelage il faut compter environs deux heures de calcul.

Paysages newtoniens

Reprenons $P(z)=z^3-1$ et regardons le résultat :

Ceci montre bien les positions des racines comme des puits profonds. Pour des points qui sont assez prêts des racines, le nombre d’itérations $u$ n’est pas très élevé, donc $w=-k/u^3$ est assez grand en valeur absolu. On remarque aussi qu’il y a des puits (moins profonds) dans toutes les zones. Cela indique des points qui convergent très rapidement vers les racines de l’équation. Dans ce cas, les racines sont $1, e^{i2\pi/3}, e^{-i2\pi/3}$. Ces valeurs comme points de départ de l’itération convergent après une itération. Si on commence l’itération avec $-\frac{1}{2}, \frac{1}{2}e^{i\pi/3}$ ou $\frac{1}{2}e^{-i\pi/3}$, on trouve que l’itération converge vers une des racines après juste deux itérations. On trouve aussi des puits dont le nombre d’itération est trois, quatre, cinq, etc.

Cette façon de dessiner ce paysage est donc intéressante d’un point de vue didactique mais elle l’est moins d’un point de vue artistique. Les escaliers ne s’accordent pas très bien avec l’aspect plus organique des bords fractals entre les zones d’attraction.

On a défini le relief en fonction du nombre d’itérations qui est par définition entier. On essaie alors de définir un nombre d’itérations « continu », comme ceci :
\[u_{L}=u+1+ \frac{\ln\left( \frac {\ln(\epsilon)} {\ln(|z_{n+1}-z_{n}|) } \right) } {\ln(2)},\]
ce qui est lié à la convergence quadratique de l’itération.

On a maintenant un relief lisse dans les différentes zones :

Avec les images fractales, il est toujours intéressant de zoomer sur des petits détails pour découvrir la répétition à l’infini des structures. Malheureusement, notre surface à base du nombre d’itérations ne se prête pas très bien aux zooms. Cela donne des problèmes près des crêtes.

Une fonction pour les reliefs qui marche mieux consiste à prendre la profondeur des vallées comme une fonction de l’estimation de distance. On peut simplement prendre $w=-k.ED$, mais pour éviter que la profondeur croisse trop « loin » des crêtes, on peut aussi prendre (avec $k$ et $m$ des constantes à choisir) :
\[w=-k.\frac{(1-e^{-m.ED})}{(1+e^{-m.ED})}\]

Toujours pour $P(z)=z^3-1$, cela donne :

On voit tout d’abord que les racines sont devenues des pointes au lieu des puits ! C’est à cause de la formule pour $ED$ qui s’annule pour les racines.

Nous pouvons maintenant zoomer comme nous le souhaitons, la limite étant imposée par la précision des calculs. (Le logiciel de l’auteur travaille normalement avec une précision de 15 décimales. On peut augmenter la précision à plus de 20 décimales, mais alors les temps de calculs sont considérablement plus longs...).

Voici quelques autres paysages, d’abord avec $P(z)=z^5-1$.

Un zoom sur $P(z)=z^8-180z^6+378z^4-420z^2+315$ :

Et finalement un zoom sur $P(z)=z^8+15z^4-16$ :

Le lapin de Douady

Dans l’article déjà cité, Tan Lei a montré comment on peut construire un fractal newtonien où apparaît un ensemble de Julia qu’Adrien Douady a nommé le « lapin » (Figure à droite). Cet ensemble de Julia dans le fractal de Newton est une zone où l’itération ne converge vers aucune des racines, donc une zone où la méthode échoue. On peut construire ceci en 3D aussi.

Prenons l’équation \[P(z)=z^3+(c-1)z-c\] Il est clair que $z=1$ est une racine. Les deux autres racines dépendent de la valeur de $c$. On soumet l’équation à la méthode de Newton mais cette fois-ci, on prend comme valeur initiale $z=0$ et on met $c$ égal à la valeur complexe des coordonnées du pixel. Cela veut dire qu’on a une équation différente pour chaque point du plan.

On obtient ceci [3] :

On fait un zoom sur le cadre :

On voit que la figure contient beaucoup d’ensembles de Mandelbrot... et on zoome encore une fois sur l’un d’eux :

On sait que l’ensemble de Mandelbrot est comme une carte. Chaque point dans l’ensemble engendre un ensemble de Julia. Pour un lapin de Douady on peut prendre le point qui est marqué en rouge sur l’image ci-dessus. C’est un endroit qui va générer un ensemble de Julia type « lapin » dans l’ensemble de Mandelbrot classique. Les coordonnées de ce point rouge sont $(0.292855027662,-1.60534922008)$.

Si maintenant on prend l’équation $P(z)=z^3+(c-1)z-c$ avec $c=0.292855027662-1.60534922008.i$, et si on effectue l’itération normale de la méthode de Newton, on obtient un nouveau paysage décoré du lapin !

On peut ainsi construire les ensembles de Julia plus ou moins comme on veut. En voici un autre :

Avec d’autres équations, on peut construire d’autres ensembles de Julia. L’image ci-dessous est générée par $P(z)=z^6+z+c$ :

Pour une plus grande collection d’images newtoniennes, on peut visiter cette page.