Frattali
"Caos: quando il presente determina il futuro, ma un'approssimazione del presente non determina approssimativamente il futuro."
Immaginiamo di far rotolare una biglia lungo un piano inclinato. Indichiamo con $x$ la posizione da cui la facciamo partire lungo il bordo alla sommità del piano e indichiamo con $y$ la posizione che raggiunge alla base del piano. Se facciamo variare la posizione di partenza della biglia lungo il bordo superiore del piano di una quantità $\Delta x$ la posizione che essa raggiungerà alla base del piano varierà a sua volta di una quantità $\Delta y = \Delta x$.
Se sostituiamo il piano inclinato con una superficie curva la variazione $\Delta y$ dipenderà dalla variazione $\Delta x$ in modo più complicato. Ma se la curvatura della superficie è abbastanza regolare, saremo sempre in grado dato un $\varepsilon$ arbitrariamente piccolo di trovare un $\delta$ sufficientemente piccolo tale per cui se $\Delta x < \delta$ allora $\Delta y < \varepsilon$.
Comportamenti caotici
Se facciamo rotolare una biglia lungo una superficie estremamente complicata, immaginiamoci ad esempio un pendio solcato dall'erosione, questa caratteristica verrà meno. Una piccolissima variazione della posizione iniziale potrebbe essere sufficiente per cambiare del tutto la sorte della biglia. Ad esempio facendola partire da $x$ potrebbe giungere a valle all'estrema sinistra della rampa, mentre facendola partire da $x + \Delta x$ potrebbe raggiungere la base nel lato opposto.
In questo caso dato un $\varepsilon$ sufficientemente piccolo, vale a dire richiedendo che la biglia raggiunga una posizione $y + \Delta y$ molto vicina a $y$ variando la posizione iniziale, non saremo in grado di trovare un $\delta$ tale per cui $\Delta y < \varepsilon$ se $\Delta x < \delta$. In questo caso parleremo di comportamento caotico.
PS: Lo studio dei sistemi caotici è importante tra le altre cose nel cosiddetto problema degli n corpi (vedi Newton.exe e Lagrange).
Il moto della biglia è un sistema continuo, possiamo rappresentarlo come un'infinità di trasformazioni infinitesime regolate da un'equazione differenziale. Nei sistemi discreti invece, si applica una certa trasformazione finita ripetute volte, iterandone l'azione. Per continuare con il nostro esempio possiamo immaginare una trasformazione $f$ che corrisponda al moto della biglia per un piccolo intervallo di tempo. $f(x)$ sarà la posizione della biglia dopo la prima iterazione, $f \big( f(x) \big)$ dopo la seconda, eccetera.
PS: Non c'entra nulla, ma la biglia che rotola è una metafora appopriata anche per l'algoritmo di apprendimento delle reti neurali.
Iterazione di funzioni complesse
Molti frattali, incluso il famoso insieme di Mandelbrot, sono sistemi discreti che mostrano comportamenti caotici generati iterando funzioni a variabile complessa. Se la rappresentazione geometrica dei numeri complessi non ti è familiare potresti dare un'occhiata al mio programma Argand.exe prima di continuare. Se ti interessa la definizione generale di frattale, non necessariamente legata a funzioni a variabile complessa, potresti invece dare un'occhiata alla pagina di R-Paint, in cui tra le altre cose descrivo le differenti definizioni di dimensione (algebrica, topologica e frattale). Infine se ti interessano le funzioni a variabile complessa in generale, e non solo nel caso specifico riguardante i frattali, potrebbe interessarti il mio programma ImMap.
Equicontinuità
La definizione di insieme di Julia per l'iterazione di una certa funzione $f(z)$ riguarda i punti per i quali piccoli cambiamenti dell'argomento $z$ causano grandi cambiamenti nei risultati. Il complementare dell'insieme di Julia, l'insieme di Fatou, è costituito dai punti per i quali piccoli cambiamenti di $z$ causano piccoli cambiamenti nei risultati. Per quest'ultimo è possibile dare una definizione rigorosa che ricalca la definizione di funzione continua, ma che in questo caso si applica ad una famiglia di funzioni (le funzioni ottenute iterando la funzione di partenza).
"Sirena"
Questa definizione rigorosa è data dal concetto di equicontinuità. Per facilitarne la comprensione formale ricordiamo la classica definizione "epsilon-delta" di continuità di funzione reale a variabile reale che tutti abbiamo studiato (o avremmo dovuto studiare) a scuola e che fa uso della definizione di limite. Una funzione reale avente dominio $X \subseteq \mathbb{R}$ si dice continua in $x_0 \in X$ se $\forall \varepsilon > 0 , \exists \delta > 0$ tale che $\forall x \in X$ $$ \left| x_0 - x \right| < \delta \Rightarrow \Big| f \left( x_0 \right) - f \left( x \right) \Big| < \varepsilon $$
A questo punto la definizione di equicontinuità per una famiglia di funzioni complesse dovrebbe essere abbastanza banale, essendo del tutto analoga. Una famiglia $\mathcal{F}$ di funzioni da $X \subseteq \mathbb{C}$ a $\mathbb{C}$ si dice equicontinua in $x_0 \in X$ se $\forall \varepsilon \in \mathbb{R}^+ , \exists \delta \in \mathbb{R}^+$ tale che $\forall f \in \mathcal{F} , \forall x \in X$ $$ d \left( x_0 , x \right) < \delta \Rightarrow d \Big( f \left( x_0 \right) , f \left( x \right) \Big) < \varepsilon $$ dove $d$ è l'usuale metrica indotta dalla norma data dal modulo dei numeri complessi. Nel nostro caso $\mathcal{F} = \left\{ f_n : n \in \mathbb{N} \right\}$ dove $f_n$ è l'$n$-esima iterazione della funzione, quindi $f_1(z) = f(z)$, $f_2(z) = f \Big( f(z) \Big)$, eccetera.
Così abbiamo una definizione rigorosa di "regolare" e definiamo il "caos" come il suo complementare.
"Tim Burton"
Alcuni esempi
Applicando questa idea si possono rappresentare gli insiemi di Julia per molte interessanti funzioni. Alcune di quelle che ho trovato particolarmente interessanti finora sono: $$ f(z) = \frac{\left( z + c \right)^2}{z^3 + c} $$ $$ f(z) = \frac{z^3 + c}{z - c} $$ $$ f(z) = z + \frac{z + c}{z^5 + c} $$ che possono essere viste come versioni più elaborate della classica $f(z) = z^2 + c$.
Ci sono poi tutta una serie di funzioni in cui si introduce l'operazione di potenza per un numero complesso. $$ f(z) = c^z $$ $$ f(z) = \sqrt[z]{c} $$ $$ f(z) = z^z + c $$ $$ f(z) = \sin(z) + c $$ Tutte queste funzioni mostrano una particolare forma "a cactus" (non sono riuscito a scoprire se questa particolare forma ha un nome più serio, che di certo si meriterebbe). Ho incluso anche la funzione seno nell'elenco in quanto anche essa mostra la medesima figura, infatti le funzioni trigonometriche sono strettamente legate all'esponenziale complessa, in particolare si ha infatti $$ \sin(z) = \frac{e^{i \, z} - e^{-i \, z}}{2 \, i} $$ (per alcune utili note riguardo ai calcoli con i numeri complessi vedi Appendice A).
L'algoritmo
Il metodo che ho scelto per individuare i valori di $z$ a cui corrispondono comportamenti regolari (e di conseguenza individuare i valori che portano a comportamenti caotici), consiste nel fissare un $\delta$ (non importa il valore preciso, ma soltanto l'ordine di grandezza) e verificare il comportamento della somma dei moduli delle differenze tra le immagini della funzione in $z$ e $z + \delta$. Avere differenze che crescono molto iterando la funzione non significa necessariamente essere di fronte ad una discontinuità, quindi a rigore il metodo non individua effettivamente i valori non equicontinui, ma è un buon compromesso se si vuole ottenere immagini in tempi brevi (magari per ottenere animazioni in realtime, visto che spesso è molto interessante vedere come l'immagine cambia variando $c$).
In pratica si calcola
$$ \varepsilon = \sum_n \Big| f_n(z + \delta) - f_n(z) \Big| $$
sommando un numero di termini sufficientemente grande.
La scelta del valore di $\delta$ dipende da quanto grande è l'area che si vuole visualizzare nel piano complesso. Un valore tipico per un frattale con $z$ che prende i valori da $-1 - i$ a $1 + i$ potrebbe essere $\delta = 0.001$. Il numero di iterazioni $N$ è invece nell'ordine delle centinaia, da un valore minimo di $N = 100$ fino anche a $500$ (ovviamente aumentando il numero di iterazioni i tempi di elaborazione si allungano notevolmente, cosa che potrebbe rendere difficile ottenere animazioni in realtime). Il valore ottimale per entrambi può dipendere dalla specifica $f(z)$.
Questa è la funzione principale del codice che ho scritto per generare i frattali mostrati in questa pagina:
// epsilon = SUM_n |f_n(z + delta) - f_n(z)|
vec3 fractal(vec2 z, vec2 c)
{
vec2 near = z + vec2(DELTA);
float epsilon = 0.0;
for (int n = 0; n < N; ++n)
{
z = f(z, c);
near = f(near, c);
epsilon += distance(z, near);
}
float hue = 0.5 + phase(z) / TWICE_PI;
return hsl2rgb(hue, S, L(epsilon));
}
(ancora, per i dettagli riguardo al codice, vedi Appendice A, mentre per il senso delle ultime due righe che calcolano il colore vedi il paragrafo successivo).
Il colore
Una volta ottenuto per ogni pixel (e quindi per ogni valore di $z$) il valore di $\varepsilon$ per un certo $\delta$ si può passare a calcolare il colore vero e proprio. La luminosità è una funzione di $\varepsilon$. Bisogna scegliere una funzione che dato un $\varepsilon \in \left[ 0 , \infty \right)$ restituisca un valore per la luminosità $L \in \left[ 0 , 1 \right]$. Ci sono molte differenti possibilità, la più semplice che ho usato finora è $$ L = 1 - \frac{1}{1 + \varepsilon} $$ L'importante, da un punto di vista puramente estetico, è cercare di avere una funzione che non restituisca semplicemente valori prossimi a $0$ per le regioni regolari e prossimi a $1$ per le regioni caotiche, ma che restituisca una buona gamma di valori intermedi nel confine tra queste regioni. Questo dipende molto anche dalla specifica $f(z)$, quindi anche la scelta della funzione $L(\varepsilon)$ può cambiare in base allo specifico frattale.
La tonalità indicata con $H$ (dall'inglese hue) è la fase di $f_n(z)$ con $n = N$, quindi $$ H = \arg \left( f_N(z) \right) $$
Infine la saturazione viene scelta a seconda dell'immagine specifica, molto spesso semplicemente 50%, quindi $S = 0.5$. Passare dalle coordinate HSL (hue, saturation, lightness) alle coordinate RGB necessarie per impostare fisicamente il colore del pixel è molto semplice (vedi Appendice B).
Conclusioni
Ricapitolando i parametri che determinano i frattali in questa pagina, in ordine di importanza per l'influenza che hanno sul risultato finale, sono:
- la scelta della funzione $f(z)$
- la scelta del parametro $c$ contenuto nella funzione
- la regione da disegnare (tipicamente intorno all'origine del piano complesso)
- il valore di $\delta$
- il numero di iterazioni $N$
- la funzione $L(\varepsilon)$ per la luminosità
- il valore $S$ per la saturazione
Risulta particolarmente interessante vedere come il frattale cambia variando il parametro $c$. In questa playlist di YouTube puoi trovare alcuni esempi (il primo video è una carrellata di immagine statiche ottenute da differenti funzioni, i successivi sono tutti video registrati in realtime).
Appendice A: operazioni con i numeri complessi in GLSL
La generazione dei frattali sembra fatta apposta per gli shader delle schede video, che eseguono in parallelo i calcoli che altrimenti dovrebbero essere eseguiti pixel per pixel. Per realizzare i miei shader ho scelto il linguaggio GLSL di OpenGL. Il codice riportato potrebbe essere ulteriormente ottimizzato, ma ho preferito scegliere un compromesso tra efficienza e leggibilità.
Iniziamo con delle semplici funzioni che non fanno altro che wrappare funzioni builtin in GLSL, in modo da esplicitarne il significato quando i nostri vettori bidimensionali rappresentano numeri complessi.
// |z| = (x^2 + y^2)^(1 / 2)
float modulus(vec2 z)
{
return length(z);
}
// |z|^2 = x^2 + y^2
float squareModulus(vec2 z)
{
return dot(z, z);
}
// arg(z) = arctan(y / x)
// [-PI, PI]
float phase(vec2 z)
{
return atan(z.y, z.x);
}
Definiamo poi una funzione per passare dalla rappresentazione polare alla rappresentazione cartesiana di un numero complesso.
// rho * (cos(theta) + i * sin(theta))
vec2 fromPolar(float rho, float theta)
{
return rho * vec2(cos(theta), sin(theta));
}
Continuiamo definendo alcune funzioni per le operazioni di base con i numeri complessi: la coniugazione, l'inverso, la moltiplicazione e la divisione (se fosse necessario qualche chiarimento rimando alla pagina di Argand.exe). Per la somma non è necessario fare nulla, visto che la somma di numeri complessi corrisponde esattamente alla somme di vettori che è disponibile nativamente in GLSL.
// z* = x - i * y
vec2 conjugate(vec2 z)
{
return vec2(z.x, -z.y);
}
// 1 / z = z* / (z z*) = z* / |z|^2
vec2 inverse(vec2 z)
{
return conjugate(z) / squareModulus(z);
}
// z * w
vec2 multiply(vec2 z, vec2 w)
{
return vec2(z.x * w.x - z.y * w.y,
z.x * w.y + z.y * w.x);
}
// i * z
vec2 multiplyByI(vec2 z)
{
return vec2(-z.y, z.x);
}
// z / w = z * (1 / z)
vec2 divide(vec2 z, vec2 w)
{
return multiply(z, inverse(w));
}
Passiamo alle operazioni che coinvolgono l'elevamento a potenza e il logaritmo, che poi sono quelle con cui si possono ottenere risultati davvero interessanti. Per la funzione esponenziale, utilizzando semplicemente la proprietà del prodotto di potenze con la stessa base, si ha $$ e^z = e^{\text{Re}(z) + i \, \text{Im}(z)} = e^{\text{Re}(z)} e^{i \, \text{Im}(z)} $$ quindi il modulo risulta $|e^z| = e^{\text{Re}(z)}$ e la fase semplicemente $\arg(e^z) = \text{Im}(z)$.
// e^z = e^(x + i y) = e^x * e^(i * y)
vec2 exponential(vec2 z)
{
return fromPolar(exp(z.x), z.y);
}
Per il logaritmo naturale (logaritmo in base $e$) si ha per definizione $$ w = \ln(z) \iff e^w = z $$ e potendo scrivere nuovamente $$ z = e^w = e^{\text{Re}(w) + i \, \text{Im}(w)} = e^{\text{Re}(w)} e^{i \, \text{Im}(w)} $$ si vede che $e^{\text{Re}(w)} = |z| \Rightarrow \text{Re}(w) = \ln|z|$ e che $\text{Im}(w) = \arg(z)$, quindi $$ \ln(z) = \ln|z| + i \, \arg(z) $$ Notiamo una cosa molto particolare: l'esponenziale e il logaritmo complessi in un certo senso invertono i ruoli di parte reale e immaginaria con modulo e fase.
// ln(z) = ln|z| + i * arg(z)
// (principal value)
vec2 naturalLogarithm(vec2 z)
{
return vec2(log(modulus(z)), phase(z));
}
Il logaritmo per una base arbitraria può facilmente essere ricondotto al logaritmo naturale. Sempre per la definizione di logaritmo si ha $$ w^{\log_w(z)} = z $$ e prendendo il logaritmo naturale di ambo i membri si ottiene $$ \ln(w^{\log_w(z)}) = \ln(z) $$ da cui, essendo che $e^{\ln(w)} = w \Rightarrow \left( e^{\ln(w)} \right)^u = w^u \Rightarrow e^{u \ln(w)} = w^u \Rightarrow u \ln(w) = \ln(w^u)$, si ha $$ \log_w(z) \ln(w) = \ln(z) $$ e quindi $$ \log_w(z) = \frac{\ln(z)}{\ln(w)} $$ (non abbiamo usato nessuna proprietà di $e$, quindi la cosa è valida per una base qualsiasi).
// log_w(z) = ln(z) / ln(w)
vec2 logarithm(vec2 z, vec2 w)
{
return divide(naturalLogarithm(z), naturalLogarithm(w));
}
Infine la potenza con base e esponente complessi può essere facilmente ricondotta all'esponenziale e al logaritmo naturale. Essendo per definizione $z = e^{\ln(z)}$, si ha $$ z^w = \left( e^{\ln(z)} \right)^w = e^{\ln(z) \, w} $$
// z^w = (e^ln(z))^w = e^(ln(z) * w)
vec2 power(vec2 z, vec2 w)
{
return exponential(multiply(naturalLogarithm(z), w));
}
La sezione "matematica" dello shader si conclude con la funzione che calcola il seno complesso $$ \sin(z) = \frac{e^{i \, z} - e^{-i \, z}}{2 \, i} $$
// sin(z) = (e^(i * z) - e^(-i * z)) / (2 * i)
vec2 sine(vec2 z)
{
vec2 iz = multiplyByI(z);
return divide(exponential(iz) - exponential(-iz), vec2(0.0, 2.0));
}
Appendice B: conversione da HSL a RGB
La rappresentazione dei colori tramite le coordinate RGB (red, green, blue) è la più immediata da un punto di vista tecnico, ma non è molto intuitiva dal punto di vista della percezione umana. Dato un colore non è così facile indovinare quale siano approssimativamente le tre componenti RGB. Esiste un differente sistema di coordinate molto più vicino all'intuito, il sistema HSL (hue, saturation, lightness).
Diamo prima una descrizione della trasformazione da RGB a HSL, che a mio parere permette di capire poi più rapidamente l'algoritmo di conversione da HSL a RGB. In un sistema di riferimento cartesiano immaginiamo un cubo di lato $1$ avente un vertice centrato nell'origine e i $3$ spigoli adiacenti a tale vertice in corrispondenza dei $3$ assi. Facciamo corrispondere ad ognuno dei punti del cubo un colore le cui coordinate RGB sono individuate dalle coordinate cartesiane del punto. Avremo il vertice $(0, 0, 0)$ che corrisponde al nero, il vertice $(1, 0, 0)$ al rosso puro, eccetera. Consideriamo la diagonale che congiunge il vertice nero $(0, 0, 0)$ con il vertice bianco $(1, 1, 1)$. I punti che giacciono in questa diagonale corrispondono alla scala di grigi. Questi punti hanno saturazione $0$ e un valore di luminosità che va da $0$ (per il nero) all'$1$ (per il bianco). La tonalità è indifferente, essendo nulla la saturazione. Allontanandoci dalla diagonale aumentiamo la saturazione e la tonalità corrisponderà all'angolo rispetto alla diagonale. Saturazione e tonalità corriponderebbero quindi al raggio $\rho$ e all'angolo $\phi$ di un sistema di coordinate cilindriche il cui asse $z$ corrisponde alla diagonale considerata, con la differenza che in realtà abbiamo a che fare con un prisma esagonale (la proiezione del cubo sul piano perpendicolare alla diagonale è un esagono) e non un cilindro. Per convenzione si fa corrispondere l'angolo $0$ con il rosso.
[...]
// HSL => RGB
vec3 hsl2rgb(float hue, float saturation, float lightness)
{
float hexhue = 6.0 * hue;
float frahue = fract(hexhue);
vec3 color;
if (hexhue < 1.0) color = vec3( 0.5, -0.5 + frahue, -0.5);
else if (hexhue < 2.0) color = vec3( 0.5 - frahue, 0.5, -0.5);
else if (hexhue < 3.0) color = vec3(-0.5, 0.5, -0.5 + frahue);
else if (hexhue < 4.0) color = vec3(-0.5, 0.5 - frahue, 0.5);
else if (hexhue < 5.0) color = vec3(-0.5 + frahue, -0.5, 0.5);
else color = vec3( 0.5, -0.5, 0.5 - frahue);
float chroma = saturation * (1.0 - abs(2.0 * lightness - 1.0));
return color * chroma + lightness;
}
