Le code dans cet article a été écrit avec Code::Blocks et la SDL 2.
Vous pouvez trouver
ici un guide pour installer ces logiciels.
Bien qu'il soit basé sur la SDL, je n'utilise pas ses fonctions directement. J'ai écrit une petite bibliothèque
avec quelques fonctions basiques pour faciliter la compréhension et la portabilité dans un autre langage.
Vous pouvez en apprendre plus sur cette bibliothèque
ici.
Pour cet article j'ai mis à jour la classe des nombres complexes
La
méthode de Newton est un algorithme pour trouver les zéros d'une fonction en appliquant à plusieurs reprises
cette formule:
xn+1 = xn - f(xn) / f'(xn)
Où f'(x) est la dérivée de la fonction f(x).
Pour transformer ça en fractale on doit l'appliquer à des nombres complexes.
Comme pour l'ensemble de Mandelbrot l'écran va représenter le plan complexe et chaque pixel sera un nombre complexe.
Et on va appliquer la formule jusqu'à ce que l'on atteigne un nombre d'itérations maximum.
Mais comme pour l'ensemble de Mandelbrot on peut sortir de cette boucle d'itération si une condition est remplie.
Pour l'ensemble de Mandelbrot c'était quand la valeur de z devenait trop grande.
Mais pour une fractale de Newton, la plupart des pixels vont tendre vers un zéro de la fonction alors on va tester
quand ils sont assez proche d'un de ces zéros.
Typiquement pour une fractale de Newton de base, on utilise la fonction:
f(z) = z3 - 1
Mais on verra plus tard qu'on peut essayer avec d'autres fonctions.
On aura aussi besoin de sa dérivée:
f'(z) = 3z2
Et on doit connaitre ses racines:
z1 = 1
z2 = -0.5 + i * √3 / 2
z3 = -0.5 - i * √3 / 2
Sur la page anglaise de la
fractale de Newton il y a un morceau de code pour la dessiner.
Et on va transformer ça en un programme qui marche.
Comme pour l'ensemble de Mandelbrot on va avoir besoin d'un nombre maximum d'itération et d'une valeur limite pour
sortir de la boucle.
#define MAX_ITERATIONS 100
#define TOLERANCE 1e-6
Ensuite, on définit notre fonction et sa dérivée.
// the function
Complex f(Complex z)
{
// z^3 - 1
static Complex one(1, 0);
return z * z * z - one;
}
// derivative of the function
Complex df(Complex z)
{
// 3 * z^2
return (z * z) * 3;
}
Dans la fonction main on va définir quelques valeurs pour dire où la fenêtre sera centrée et sa taille dans le plan
complexe.
static double zoom = 2.0 / (double)SCREEN_HEIGHT;
static double centerX = (double)SCREEN_WIDTH / 2.0;
static double centerY = (double)SCREEN_HEIGHT / 2.0;
On va avoir besoin des racines de la fonction:
// roots of the function
Complex roots[3];
roots[0] = Complex(1, 0);
roots[1] = Complex(-0.5, sqrt(3) / 2);
roots[2] = Complex(-0.5, -sqrt(3) / 2);
Chaque pixel prendra une couleur en fonction de la racine qu'on a atteinte dans la boucle d'itération.
Alors on va définir une palette. Mais on va y ajoute une valeur à la fin pour les pixels qui n'ont atteint aucune
racine après le nombre max d'itérations.
// palette
static Color palette[4];
palette[0] = Color(255, 0, 0);
palette[1] = Color(0, 255, 0);
palette[2] = Color(0, 0, 255);
palette[3] = Color(0, 0, 0);
Maintenant on va dessiner la fractale. Alors on commence une boucle sur tous les pixels de l'écran:
// draw the fractal
for (int y = 0; y < SCREEN_HEIGHT; y++)
{
for (int x = 0; x < SCREEN_WIDTH; x++)
{
On convertit les coordonnées du pixel en nombre complexe:
// convert the pixel pos to a complex number
Complex z;
z.re = (x - centerX) * zoom - 0.75;
z.im = (y - centerY) * zoom;
On va utiliser une variable pour stocker la racine qu'on a atteint. Par défaut elle est sur la dernière couleur de
la palette.
// by default the pixel is black if the calculations fail
int color = 3;
On lance la boucle d'itération
// iterate
int iter;
for (iter = 0; iter < MAX_ITERATIONS; iter++)
{
Voilà la formule de Newton
z = z - f(z) / df(z);
On boucle sur chaque racine pour tester si on en a atteint une
// are we near a root ?
for (int i = 0; i < 3; i++)
{
Complex delta = z - roots[i];
double distance = delta.squareMag();
Si on a atteint une racine, on stocke son numéro et on sort de la boucle de test des racines.
if (distance < TOLERANCE)
{
color = i;
break;
}
}
Et on sort aussi de la boucle d'itération:
if (color != 3)
break;
}
Et finalement on dessine le pixel.
// draw the pixel
gfx.setPixel(x, y, palette[color]);
Comme le dessin peut prendre du temps, à la fin de chaque ligne, on va l'afficher et gérer les évènements pour
tester si l'utilisateur a fermé la fenêtre.
On fait ça pour éviter que le programme n'ait pas l'air de répondre pendant le dessin.
}
gfx.render();
sys.processEvents();
if (sys.isQuitRequested() == true)
break;
}
A la fin, on rajoute une boucle pour attendre que l'utilisateur ferme la fenêtre quand le dessin est terminé.
// wait until we quit
while (sys.isQuitRequested() == false)
{
sys.processEvents();
sys.wait(16);
}
gfx.quit();
return EXIT_SUCCESS;
Télécharger le code source
Télécharger l'exécutable pour Windows
Et voila le résultat:
Cette image a l'air un peu plate.
Alors on va changer la couleur des pixels en fonction du nombre d'itérations qu'on a faites.
Remarque: j'ai renommé la variable "color" en "palIdx" parce que c'est plus parlant.
Maintenant, avant qu'on affiche le pixel, on va mettre ce code:
// draw the pixel
Color color = palette[palIdx];
if (palIdx != 3)
{
float coef = 1 - log(iter + 1) / 5.0;
if (coef < 0)
coef = 0;
color.r *= coef;
color.g *= coef;
color.b *= coef;
}
gfx.setPixel(x, y, color);
Télécharger le code source
Télécharger l'exécutable pour Windows
Je calcule un coefficient à partir du logarithme du nombre d'itérations.
Ce coefficient devrait rester entre 0 et 1.
Ensuite je multiplie les valeurs R, V, B de la couleur par ce coefficient.
Le résultat est de que la couleur va devenir plus sombre quand le nombre d'itérations va augmenter.
Et voila l'image qu'on obtient:
Maintenant on peut voir les racines de la fonction.
Ce sont les points clairs qui sont dans les grandes zones à l'extérieur de la fractale.
Maintenant je vais essayer de reproduire la plupart des images qu'on peut voir sur la
page wikipedia de la fractale
de Newton.
On va commencer avec l'image z
3 - 2z + 2.
Si vous voulez changer la fonction dans ce code vous aurez aussi besoin de sa dérivée et de ses racines.
On verra d'autres méthodes qui évitent ça, mais pour l'instant une rapide visite sur
Wolfram Alpha devrait vous
donner tout ce dont vous avez besoin.
Alors voilà ce qu'on doit changer dans le code:
// the function
Complex f(Complex z)
{
// z^3 - 2z + 2
return z * z * z - z * 2 + Complex(2, 0);
}
// derivative of the function
Complex df(Complex z)
{
// 3 * z^2 - 2
return (z * z) * 3 - Complex(2, 0);
}
[...]
// roots of the function
Complex roots[3];
roots[0] = Complex(-1.7693, 0);
roots[1] = Complex(0.88465, 0.58974);
roots[2] = Complex(0.88465, -0.58974);
Télécharger le code source
Télécharger l'exécutable pour Windows
J'ai aussi changé les couleurs et le coefficient d'ombrage pour être plus proche de l'image d'origine et voila ce
qu'on obtient:
La particularité avec cette fonction c'est que pour certaines valeurs de départ on n'arrive pas à une racine.
Ces zones sont en rouge sur l'image.
Celle-ci a plus de racines, ça implique quelques changements quand on teste la variable palIdx.
Mais à part ça il n'y a pas beaucoup de changements.
// the function
Complex f(Complex z)
{
// z^8 + 15z^4 - 16
Complex z4 = z*z*z*z;
return z4 * z4 + z4 * 15 - Complex(16, 0);
}
// derivative of the function
Complex df(Complex z)
{
// 8z^7 + 60z^3
Complex z3 = z*z*z;
Complex z4 = z3 * z;
return (z4 * z3) * 8 + z3 * 60;
}
[...]
// roots of the function
Complex roots[8];
roots[0] = Complex(-1, 0);
roots[1] = Complex(1, 0);
roots[2] = Complex(0, -1);
roots[3] = Complex(0, 1);
roots[4] = Complex( 1.4142, 1.4142);
roots[5] = Complex(-1.4142, 1.4142);
roots[6] = Complex( 1.4142, -1.4142);
roots[7] = Complex(-1.4142, -1.4142);
Télécharger le code source
Télécharger l'exécutable pour Windows
f(z) = z5 - 3iz3 - (5+2i)z2 + 3z + 1
// the function
Complex f(Complex z)
{
// z^5 - 3iz^3 - (5+2i)z^2 + 3z + 1
return z * z * z * z * z - Complex(0, 3) * z * z * z - Complex(5, 2) * z * z + z * 3 + Complex(1, 0);
}
// derivative of the function
Complex df(Complex z)
{
// 5z^4 - 9iz^2 - (10+4i)z + 3
return z * z * z * z * 5 - Complex(0, 9) * z * z - Complex(10, 4) * z + Complex(3, 0);
}
[...]
// roots of the function
Complex roots[5];
roots[0] = Complex(-1.28992, -1.87357);
roots[1] = Complex(-0.824853, 1.17353);
roots[2] = Complex(-0.23744, 0.0134729);
roots[3] = Complex(0.573868, -0.276869);
roots[4] = Complex(1.77834, 0.963437);
Télécharger le code source
Télécharger l'exécutable pour Windows
Celle là a l'air à l'envers par rapport à l'image de wikipedia, mais la forme est la même.
// the function
Complex f(Complex z)
{
// z^6 + z^3 - 1
Complex z3 = z * z * z;
return z3 * z3 + z3 - Complex(1, 0);
}
// derivative of the function
Complex df(Complex z)
{
// 6z^5 + 3z^2
return z * z * z * z * z * 6 + z * z * 2;
}
[...]
// roots of the function
const int nbRoots = 6;
Complex roots[nbRoots];
roots[0] = Complex(-1.17398499670533, 0);
roots[1] = Complex(0.851799642079243, 0);
roots[2] = Complex(0.586992498352664, 1.016700830808605);
roots[3] = Complex(0.586992498352664, -1.016700830808605);
roots[4] = Complex(-0.425899821039621, -0.737680128975117);
roots[5] = Complex(-0.425899821039621, 0.737680128975117);
Télécharger le code source
Télécharger l'exécutable pour Windows
Celle là n'a pas l'air aussi symétrique que celle de wikipedia. Je ne sais pas pourquoi.
Ce n'est peut-être pas la même partie de la fractale.
La page wikipedia anglaise montre une variante de la formule où un coefficient est ajouté:
zn+1 = zn - a * f(zn) / f'(zn)
Avec cette formule les racines de la fonction ne correspondent plus.
Alors on ne va pas utiliser les tableaux roots et palette dans ce code.
La première image qui utilise cette formule utilise une méthode de coloration différente.
By Asympt (talk) - self-made, Public Domain, Link
Elle utilise la fonction z
3 - 1 avec a = -0.5 et effectue 40 itérations.
Puis on prend l'argument de z pour trouver la couleur.
J'ai essayé de la reproduire exactement et ça m'a pris un peu de temps pour trouver les bonnes formules pour la
couleur:
// draw the pixel
Color color;
color = Color(0, 0, 0);
float angle = z.arg();
if (angle <= M_PI / 3)
color.r = 127 + 128 * sin(-angle);
if (angle >= -M_PI * 2/3 && angle <= M_PI * 2/3)
color.g = 255 * sin(angle * 3/4 + M_PI / 2);
if (angle >= -M_PI / 3)
color.b = 127 + 128 * sin(angle);
gfx.setPixel(x, y, color);
Télécharger le code source
Télécharger l'exécutable pour Windows
Et voila ce que j'obtiens:
Maintenant l'image originale présente des bandes parce qu'elle utilise une palette de 64 couleurs.
On peut reproduire ça en "arrondissant" l'angle:
float angle = z.arg();
angle = angle + M_PI;
angle = (int)(angle * 32 / M_PI);
angle = angle * M_PI / 32;
angle = angle - M_PI;
Télécharger le code source
Télécharger l'exécutable pour Windows
Ensuite j'ai essayé un schéma de couleurs plus simple
// draw the pixel
Color color;
float angle = z.arg();
color.r = 127 * (1 + sin(angle));
color.g = 127 * (1 + sin(angle + 2 * M_PI / 3));
color.b = 127 * (1 + sin(angle + 4 * M_PI / 3));
gfx.setPixel(x, y, color);
Télécharger le code source
Télécharger l'exécutable pour Windows
Et voilà ce que j'ai obtenu:
Ca démontre que l'aspect d'une fractale peut changer énormément quand vous avez une palette différente.
f(z) = z2 - 1 ; a = 1 + i
Voici les changements pour la 2ième image qui utilise cette variante.
// the function
Complex f(Complex z)
{
// z^2 - 1
return (z * z) - Complex(1, 0);
}
// derivative of the function
Complex df(Complex z)
{
// 2z
return z * 2;
}
[...]
Complex a(1, 1);
Télécharger le code source
Télécharger l'exécutable pour Windows
L'étrange cas de f(z) = sin(z)
Maintenant revenons à la formule standard de Newton.
L'image montrant f(z) = sin(z) est un peu plus ardue à obtenir.
Ecrire les formules est facile:
// the function
Complex f(Complex z)
{
// sin(z)
return Complex::Sin(z);
}
// derivative of the function
Complex df(Complex z)
{
// cos(z)
return Complex::Cos(z);
}
Mais cette fonction a un nombre de racines infini de la forme n * π où n est un entier.
Alors cette fois encore on ne va pas utiliser les tableaux roots ni palette.
Mais on va utiliser une astuce pour trouver la racine qu'on a atteinte.
D'abord on divise z par π:
// are we near a root ?
Complex delta = z / Complex(M_PI, 0);
Ensuite on compare ce résultat avec un complexe créé en prenant sa partie entière:
Complex delta2 = delta - Complex((int)delta.re, 0);
double distance = delta2.squareMag();
if (distance < TOLERANCE)
{
Finalement on peut choisir une couleur qui dépend de cet entier, qui est n:
found = true;
float n = (int)delta.re * 0.072345 * M_PI;
color.r = 127 * (1 + sin(n));
color.g = 127 * (1 + sin(n + 2 * M_PI / 3));
color.b = 127 * (1 + sin(n + 4 * M_PI / 3));
Télécharger le code source
Télécharger l'exécutable pour Windows
Et après avoir ajouté notre coefficient d'ombrage on obtient cette belle image:
Remarquez que les 2 cercles noirs ici n'apparaissent qu'à cause du coefficient d'ombrage, mais la fractale a
vraiment 2 cercles noirs qu'on peut voir en zoomant un peu plus.
La ligne noire verticale devrait probablement disparaitre si vous décalez un peu les coordonnées complexes de la
fenêtre.
La méthode de Newton n'est pas le seul algorithme que les mathématiciens ont essayé pour trouver les racines d'une
fonction.
La version française de la
page wikipedia parle d'autres variantes.
La première est la méthode de la sécante où l'on remplace la formule de Newton par celle-ci:
zn+1 = zn - f(zn) * (zn - zn-1) / (f(zn) - f(zn-1))
Ici on n'a pas besoin de la dérivée de la fonction mais on va devoir utiliser des astuces parce que la formule
utilise z
n-1.
Ce programme et les suivants vont être des versions modifiées du 2ième programme de cet article.
Maintenant au lieu d'un simple "z" on aura un tableau de 3 valeurs:
// convert the pixel pos to a complex number
Complex z[3];
for (int i = 0; i < 3; i++)
{
z[i].re = (x - centerX) * zoom - 0.75 + (2-i) * TOLERANCE;
z[i].im = (y - centerY) * zoom;
}
Remarquez qu'on ajoute une petite valeur à la partie réelle, car la méthode ne marche pas si z
n-1 est égal à z
n.
Alors la formule d'itération va ressembler à ça:
z[2] = z[1] - f(z[1]) * (z[1] - z[0]) / (f(z[1]) - f(z[0]));
z[0] = z[1];
z[1] = z[2];
Télécharger le code source
Télécharger l'exécutable pour Windows
Et voici l'image qu'on obtient:
Si vous observez les nombreuses zones noires, vous pouvez en déduire que cette méthode n'est pas très efficace pour
trouver une racine de z
3 - 1
La 2ième variante décrite dans la page française utilise la méthode de Halley:
zn+1 = zn - 2 * f(zn) * f'(zn) / (2 * [f'(zn)]2 - f(zn) * f''(zn))
C'est une méthode cubique qui utilise la dérivée seconde de la fonction.
// the function
Complex f(Complex z)
{
// z^3 - 1
static Complex one(1, 0);
return z * z * z - one;
}
// derivative of the function
Complex df(Complex z)
{
// 3 * z^2
return (z * z) * 3;
}
// second derivative of the function
Complex ddf(Complex z)
{
// 6 * z
return z * 6;
}
Et la formule ressemble à ça:
Complex fz = f(z);
Complex dz = df(z);
Complex ddz = ddf(z);
z = z - fz * dz * 2 / (dz * dz * 2 - fz * ddz);
Télécharger le code source
Télécharger l'exécutable pour Windows
Le résultat est un peu moins passionnant que la fractale de Newton classique:
La méthode de Householder
La 3ième variante présentée dans la page française est la méthode de Householder.
C'est une combinaison des méthodes de Newton et de Halley:
zn+1 = zn - (1 + h) * f(zn) / f'(zn)
Où
h = f(zn) * f''(zn) / (2 * [f'(zn)]2)
La formule est implémentée comme ceci:
Complex fz = f(z);
Complex dz = df(z);
Complex ddz = ddf(z);
Complex h = (fz * ddz) / (dz * dz * 2);
z = z - (Complex(1, 0) + h) * fz / dz;
Télécharger le code source
Télécharger l'exécutable pour Windows
Et le résultat est un peu plus intéressant:
La page anglaise parle d'une fractale Nova sans en donner d'exemple.
La formule est une extension de la première variante qu'on a vue.
On n'introduit plus seulement un coefficient "a" mais on ajoute aussi une constante:
zn+1 = zn - a * f(zn) / f'(zn) + c
L'itération ressemble maintenant à ça:
// iterate
int iter;
Complex a(1, 0);
Complex c(0.3, 0.3);
for (iter = 0; iter < MAX_ITERATIONS; iter++)
{
Complex oldZ = z;
z = z - a * f(z) / df(z) + c;
Et comme avant on ne peut pas compter sur les racines de la fonction d'origine, alors on va utiliser la même
technique de coloration en fonction de l'argument.
Mais la formule converge vers une racine quand même.
Alors au lieu d'exécuter un nombre fixe d'itérations, on va sortir de la boucle de quand la différence entre
z
n+1 et z
n sera assez petite:
Complex delta = z - oldZ;
double distance = delta.squareMag();
if (distance < TOLERANCE)
break;
}
Télécharger le code source
Télécharger l'exécutable pour Windows
Et voici un exemple avec a = 1 et c = 0.3 + 0.3i
Comme la page wikipedia anglaise l'explique, il y a 2 façons de calculer une fractale Nova.
C'est la même différence qu'il y a entre l'ensemble de Mandelbrot et une fractale de Julia.
Pour la première version de la fractale Nova on a initialisé "z" avec les coordonnées du pixel et "c" est resté
constant.
Maintenant pour cette version "c" sera la position du pixel et on va initialiser z avec une valeur constante.
Voici les petites modifications apportées au code:
// iterate
int iter;
Complex a(1, 0);
Complex c = z;
z = Complex(-1, 0);
for (iter = 0; iter < MAX_ITERATIONS; iter++)
{
Complex oldZ = z;
z = z - a * f(z) / df(z) + c;
Télécharger le code source
Télécharger l'exécutable pour Windows
Et voici ce qu'on obtient avec a = 1 et z
0 = -1
Ca ressemble à des images que j'ai trouvées
ici.
