Matteo Basei

Una collezione di piccoli programmi realizzati a scopo didattico.

Il meraviglioso mondo dei frattali

"Caos: quando il presente determina il futuro, ma un'approssimazione del presente non determina approssimativamente il futuro."

Edward Norton Lorenz

Immaginiamo di far rotolare una biglia lungo un piano inclinato. Indichiamo con $x$ la posizione da cui la facciamo partire 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$.

Frattale ottenuto iterando la funzione f(z) = ((z + c)^2) / (z^3 + c).
Frattale ottenuto iterando la funzione $f(z) = (z + c)^2 / (z^3 + c)$.

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.

Un altro frattale ottenuto iterando la funzione f(z) = ((z + c)^2) / (z^3 + c) per un differente valore di c.
Un altro frattale ottenuto iterando la funzione $f(z) = (z + c)^2 / (z^3 + c)$ per un differente valore di $c$.

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. Se invece 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. Infine se ti interessano le funzioni a variabile complessa in generale potrebbe interessarti il mio programma ImMap.

Frattale ottenuto iterando la funzione f(z) = (z^3 + c) / (z - c)
Frattale ottenuto iterando la funzione $f(z) = (z^3 + c) / (z - c)$.

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).

Frattale ottenuto iterando la funzione f(z) = z + (z + c) / (z^5 + c)
Frattale ottenuto iterando la funzione $f(z) = z + (z + c) / (z^5 + c)$.

Questa definizione rigorosa è data dal concetto di equicontinuità: 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.

Frattale.
$f(z) = c^z$

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).

Frattale.
$f(z) = \sqrt[z]{c}$

L'algoritmo

[...]

$$ \varepsilon = \sum_n \Big| f_n(z + \delta) - f_n(z) \Big| $$

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)$.

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 = clamp(0.5 + phase(z) / 6.2831853071790.01.0);
    
    return hsl2rgb(hue, S, L(epsilon));
}
Frattale.
$f(z) = z^z + c$

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 + \epsilon} $$ 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(\epsilon)$ può cambiare in base allo specifico frattale.

La tonalità $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).

Frattale.
$f(z) = \sin(z) + c$

Conclusioni

Ricapitolando i parametri che determinano i frattali in questa pagina, in ordine di importanza per l'influenza che hanno sul risultato finale, sono:

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).

Frattale.
$f(z) = z + \sqrt[10]{\left( -1 \right)^z}$

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.

float modulus(vec2 z)
{
    return length(z);
}

float squareModulus(vec2 z)
{
    return dot(z, z);
}

float phase(vec2 z)
{
    return atan(z.y, z.x);
}

[...]

vec2 fromPolar(float rho, float theta)
{
    return rho * vec2(cos(theta), sin(theta));
}

Definiamo poi 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.

vec2 conjugate(vec2 z)
{
    return vec2(z.x, -z.y);
}

vec2 inverse(vec2 z)
{
    return conjugate(z) / squareModulus(z);
}

vec2 multiply(vec2 z, vec2 w)
{
    return vec2(z.x * w.x - z.y * w.y,
                z.x * w.y + z.y * w.x);
}

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)$.

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.

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).

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} $$

vec2 power(vec2 z, vec2 w)
{
    return exponential(multiply(naturalLogarithm(z), w));
}

Appendice B: conversione da HSL a RGB

[...]