F. Fasso`.    Universita` di Padova.    Laurea di I livello in Matematica.    Minicorso a.a. 2007-2008.

Julia sets

Lez. 3. Julia sets

Avvertenza: dimensione files

Attenzione: 
I grafici creati possono essere di grandi dimensioni:
- salvare spesso il notebook (meglio se con un nome diverso, per evitare che venga sovrascritto; volendo, usare il pacchetto salva.m, ma fare attenzione alla dimensione dei files: cancellare le celle di output prima di salvare)
- cancellare spesso i grafici che non servono piu`

Teoria: Iterazione di mappe polinomiali complesse: Julia sets e insieme di Mandelbrot

Si studiano le orbite 
    z,  (Subscript[P, c])(z),  ((Subscript[P, c])^2)(z) , ...., ((Subscript[P, c])^n)(z) ,  .........
della mappa complessa  (cioe` da  C a  C)  
      z -> (Subscript[P, c])(z) : = (z^2)+c
che dipende dal parametro complesso  c.  
(Si intende ((Subscript[P, c])^2)(z) = (Subscript[P, c])((Subscript[P, c])(z))  etc.)

Interessa il comportamento asintotico di queste orbite. In particolare: restano limitate o vanno all'infinito?

Per ogni  c, si definisce l'insieme dei punti che "non scappano" come
    (Subscript[K, c])={z ∈ C : (lim_(n -> ∞)) |((Subscript[P, c])^n)(z)| < ∞}
mentre l'insieme di quelli che "scappano" e` il suo complementare
    (Subscript[E, c]) =  C - (Subscript[K, c])

Esempio:  (Subscript[K, 0]) e` il disco unitario. 
A parte questo caso ed il caso  c = –2, la struttura dell'insieme  (Subscript[K, c]) e` estremamente complessa - ed interessante

Proposizione
Qualunque sia  c ∈ C,  gli insiemi    (Subscript[K, c])  ed   (Subscript[E, c])  sono non vuoti.

Dim: i punti fissi   z=(Subscript[P, c])(z)  non scappano.  I punti esterni ad un disco di raggio abbastanza grande scappano; quest'ultimo punto si prova con semplici stime, che portano al seguente:

Lemma 
Qualunque sia  c ∈ C,   (Subscript[K, c])  e` contenuto nella chiusura del disco di raggio 
        (Subscript[R, c]) =\(TraditionalForm\`\(\(\ \)\( .5 + \@\(\(\(.25\)\(+\)\) | c | \)\)\)\)  

NB: (Subscript[R, c]) ≤ max(|c|,2)  Dunque ogni punto che sta a distanza ≥ max(|c|,2) dall'origine "scappa". (Questa e` la stima fatta a lezione)

Definizione:  la frontiera  (Subscript[J, c]) di  (Subscript[K, c])  si chiama  Julia set del polinomio  (Subscript[P, c])

Teorema (Julia, Fatou)
1)  (Subscript[K, c]) e` compatto
2)  (Subscript[K, c]) e` connesso se e solo se  l'origine ∈  (Subscript[K, c])
3)  Se (Subscript[K, c]) non e` connesso, allora e` totalmente disconnesso e non numerabile ("Cantor dust").

Domanda: come sono fatti i Julia set?
Ci proponiamo di disegnarli, per diversi valori di c, e di studiarne la struttura (per es, auto-similarita`)

Algoritmo

Strategia

Il procedimento per la costruzione dei Julia set e` in linea di principio molto semplice:

1) Per stabilire (approssimativamente) se un punto z del piano complesso appartenga o no a (Subscript[K, c]) si itera la mappa un certo numero N di volte e si controlla se la sua orbita resti all'interno del cerchio di raggio (Subscript[R, c]):
- se esce da esso, allora certamente z non appartiene a (Subscript[K, c])
- se vi resta, si "stabilisce" che appartiene a (Subscript[K, c])
Se scappa, si colora il punto di un colore scelto, se non scappa di un altro.

Naturalmente, bisogna che N sia abbastanza grande. In pratica si puo` prendere N da alcune decine ad alcune centinaia, a seconda della accuratezza desiderata; per avere un'indicazione empirica del fatto che l'N scelto sia abbastanza grande si puo` provare ad aumentarlo (per es, raddoppiarlo) e controllare che i risultati restino (quasi) invariati.

2) Si ripete questa operazione su una griglia nel piano complesso C - cioe` in (R^2)

3) Si produce infine una figura, di tipo raster

Immagini raster

Vediamo innanzittutto come creare un' immagine raster, opportunamente colorata.

Cosa e` una immagine raster?
From Wikipedia, the free encyclopedia:
A raster graphics image, digital image,or bitmap, is a data file or structure representing a generally rectangular grid of pixels, or points of colour, on a computer monitor, paper, or other display medium

In Mathematica, l' immagine raster di una matrice rettangolare si puo` creare nel modo piu` semplice con il comando RasterArray.
Se c11,....,cnm sono n × m  direttive grafiche che specificano un colore, allora 
         RasterArray[ { {c11,c12,....,c1m}, ....,{cn1,....cnm}}]
produce una primitiva grafica che consiste di un reticolo di m x n quadrati (ciacuno di lato unitario) colorati secondo la direttiva grafica.  
Questa primitiva grafica deve naturalmente essere trasformata in un oggetto grafico e poi mostrata:
        Show[ Graphics[     RasterArray[ { {c11,c12,....,c1m}, ....,{cn1,....cnm}}]   ]]

Esempio:

RasterArray[   Table[
        Hue[ ((i + j)/6)], {i, 2}, {j, 2}]]

% // Graphics // Show;

Beh, non sembrano quadrati. Ma sappiamo come rimediare:

Show[Graphics[
          RasterArray[   Table[Hue[ ((i + j)/6)], {i, 2}, {j, 2}]]],AspectRatio -> Automatic]

Dimensioni del RasterArray:  se tab e` una matrice m×n allora 
            RasterArray[ { {c11,c12,....,c1m}, ....,{cn1,....cnm}]
produce un rettangolo di lati m × n. Il rettangolo di vertici  (xmin,ymin)  ed  (xmax,ymax)  si ottiene con
            RasterArray[ { {c11,c12,....,c1m}, ....,{cn1,....cnm},  { {xmin,xmax},  {ymin,ymax}} ]

Iterazione del polinomio su un punto, e colorazione

Scriviamo ora una funzione
        julia[z,numit,c]
che determina le iterazioni successive di (Subscript[P, c]) a partire da un punto z , controlla se ciascuna di essa esce dal disco di raggio (Subscript[R, c]), nel qual caso si arresta e produce un colore (= primitiva grafica di colore). Se l'iterazione supera  numit  volte senza uscire da tale disco, la procedura si arresta e produce un altro colore.

Scelta dei colori: per esempio, Hue[0] (rosso) e Hue[.7] (azzurro)

L'iterazione della mappa, ed il controllo di uscita dal cerchio di raggio (Subscript[R, c]), si puo` fare in varii modi. Uno possibile e` con una While loop. In Mathematica, la sintassi di While e` la seguente: 
               While[test, body] evaluates test, then body, repetitively, until test first fails to give True.
Da notare che
       Unless Return[ ] or Throw[ ] are generated, the final value returned by While is Null.
Pero`,  si puo` far produrre al ciclo while il valore di un iteratore.

Possiamo allora fare qualcosa del tipo:
(Inizializzazione)
      k = 1 ; w=z   
(Ciclo while)
       While[   Abs[w] < (Subscript[R, c]) , w = (Subscript[P, c])(w)  , k++  ]
(Uscita)
Quando la procedura si arresta, il numero di iterazioni e` registrato nella variabile k. Basta allora produrre
       Hue[ .7 Floor[k/numit] ]
La procedura puo` essere convenientemente assemblata in un Module.

julia[z_,numit_,c_]:=Module[{k = 1,w=z,r=.5+(Sqrt[.25+Abs[c]])},
                                                                     While[Abs[w]<r   &&   k<numit, w = (w^2) + c ; k++];  
                                                                      Hue[.7 Floor[(k/numit)]]]

Prova, con c=0 (sappiamo quale deve essere il risultato):

{julia[.9,10,0], julia[1.1,10,0]}

Valutazione di "julia" su griglia

Dobbiamo ora calcolare julia[z ,numit,c]  per tutti i punti  z  di una  griglia che riempe una certo rettangolo nel piano complesso.

Questo si fa con una Table. Si puo` rappresentare il punto z nella forma  x + I y e fare una Table bidimensionale su x ed y:
               Table[  julia[x+I y ,numit,c],  {y,ymin,ymax,dy}, {x,xmax,xmin,dx} ]
ove xmin, xmax etc determinano la regione rettangolare e  dx e dy sono i passi della griglia.
Questo produce la desiderata matrice quadrata ("Array") di primitive grafiche di colore, che si trasformera` poi nell'immagine Raster..

Scriviamo allora una funzione
             julgrid[{xmin_,xmax_,ymin_,ymax_},ndiv_,numit_,c_]
che esegue questa operazione. Conviene tenere variabili
-  gli estremi   (xmin,xmax),  (ymin,ymax)   degli intervalli che formano il rettangolo (perche` poi vorremo poter "zoomare" nella figura)
-  i passi  della griglia, o equivalentemente i numeri di righe e colonne.
Per semplicita` prenderemo uguali i numeri di righe e colonne (ndiv), cosicche` i passi sono  ((xmax - xmin)/ndiv)  ed  ((ymax - ymin)/ndiv).

julgrid[{xmin_, xmax_, ymin_, ymax_}, ndiv_, numit_, c_]:=
Table[ julia[x + I y, numit, c] ,  {y, ymin, ymax, ((ymax - ymin)/ndiv)}, {x, xmin, xmax, ((xmax - xmin)/ndiv)}]

(NB: Attenzione all'ordine degli iteratori nella Table qui sopra: prima y, poi x (se no la figura viene ruotata oppure, per evitarlo, bisogna trasporre la matrice). Per capirlo, provare con qualche saso semplice)

Funzione che crea il filled Julia set

Non resta ora che trasformare la array cosi` prodotta in una RasterArray, trasformarla in un oggetto grafico, e mostrarla. Per maggior leggibilita`, scriviamo due funzioni separate che creano l'oggetto grafico e la figura (cosi` potremo aggiungere opzioni a Show senza riscrivere tutto).

GRA[{xmin_, xmax_, ymin_, ymax_}, ndiv_, numit_, c_]:=
    Graphics[ 
     RasterArray[ julgrid [{xmin, xmax, ymin, ymax}, ndiv, numit, c] , {{xmin, ymin}, {xmax, ymax}} ]
        ]

(Attenzione al diverso ordine delle coordinate dei vertici)

JP[{xmin_, xmax_, ymin_, ymax_}, ndiv_, numit_, c_]:=
Show[
        GRA [{xmin, xmax, ymin, ymax}, ndiv, numit, 
          c], 
         AspectRatio -> 
          Automatic]

Prova

Vediamo se JP funziona bene (l'output grafico e` sufficientemente chiaro?) e se e` sufficientemente veloce.

Se c=0, (Subscript[K, c]) e` il disco unitario.
Attenzione: provare con piccoli valori di numit e, soprattutto, ndiv, facendoli crescere a poco a poco.

JP[{-1. , 1. , -1. , 1. }, 50, 30, 0. ]

Attenzione: se a c si da il valore intero 0 si ottengono dei messaggi di underflow.

Attenzione: e` molto importante che tutti i valori delle variabili reali siano numeri floating point: confrontate la velocita` di esecuzione di
            juliaplot[{-1.,1.,-1.,1.},20,10,0.]
e di 
            juliaplot[{-1,1,-1,1},20,10,0].

Miglioramento dell'immagine

JP[{-1. , 1. , -1. , 1. }, 100, 50, -1. ]; // 
      Timing

Ci sono un paio di miglioramenti da fare all'immagine creata da juliaplot:
  1. Aggiungere degli assi (o una cornice) con i valori delle coordinate
  2.  Vorremmo poter scegliere la dimensione dell'immagine

1. E` piu` semplice aggiungere una cornice, con l'opzione   Frame-> True, che gli assi:

JP[{xmin_, xmax_, ymin_, ymax_}, ndiv_, numit_, c_]:=
Show[   GRA[{xmin, xmax, ymin, ymax}, ndiv, numit, c],
         AspectRatio -> Automatic,Frame -> True]

JP[{-1. , 1. , -1. , 1. }, 100, 50, -1. ]; // 
      Timing

2. Le dimensioni si fissano con l'opzione     ImageSize. In questo caso, per avere un facile controllo dei grafici che produrremo piu` avanti, conviene specificare la lunghezza, chiamiamola    ImSize , del solo lato maggiore. Come si puo` fare?

ImageSize->ImSize  Min[1,(sidex/sidey)]

JP[{xmin_, xmax_, ymin_, ymax_},ndiv_,numit_,c_,ImSize_]:=
Show[ GRA[{xmin, xmax, ymin, ymax}, ndiv, numit, c],
           AspectRatio -> Automatic,Frame -> True,
           ImageSize -> 
              ImSize   Min[1, ((xmax - xmin)/(ymax - ymin))] ]

JP[{-2. , 2. , -1. , 1. }, 50, 30, -1, 200];

JP[{-1. , 1. , -2. , 2. }, 50, 30, -1. , 
        200];

JP[{-2. , 2. , -1. , 1. }, 300, 40, -1.  +  .3  I, 
          300]; // Timing

Provando, si osserva facilmente che in certi casi la figura viene tutta rossa: tutti i punti scappano. Il motivo e` che se il Julia set e` totalmente disconnesso, ed il numero di iterazioni e` alto, tutti i punti scappano. Un rimedio e` il seguente:

Prime prove

JP[{-2. , 2. , -1. , 1. }, 300, 40, -1.  +  .3  I, 
          300]; // Timing

Provando, si osserva facilmente che in certi casi la figura viene tutta rossa: tutti i punti scappano. Il motivo e` che se il Julia set e` totalmente disconnesso, ed il numero di iterazioni e` alto, tutti i punti scappano. Un rimedio e` il seguente:

Visualizzazione della velocita` di fuga

Per visualizzare i Julia set anche in casi nei quali essi sono "Cantor dust"  basta colorare i punti del piano con un colore che dipende dal  numero di iterazioni al quale la loro orbita scappa: in questo modo i punti di (Subscript[K, c]), anche se isolati, sono circondati da aloni colorati che li mettono in evidenza. 
Inoltre, cosi` si ottengono figure molto belle.

Funzioni

Per farlo, basta modificare l'uscita di   juliac  :  invece di    Floor[(k/numit)]   serve    (k/numit)

Non c'e` invece bisogno di modificare le funzioni che producono il rasterarray

Esempi

Riguardiamo adesso gli esempi precedenti

JP[{-2. , 2. , -1. , 1. },300,40,-1.  +  .3  I,300];//Timing

JP[{-2. , 2. , -1. , 1. },500,300,-1.  +  .3  I,300];//Timing

Zoomiamo dentro quest'ultimo

JP[{- .1,  .1, - .1,  .1},500,300,-1.  +  .3  I,300];//Timing

JP[{0.02,  .06, - .02,  .02},500,300,-1.  +  .3  I,300];//Timing

E poi

JP[{-2. , 2. , -1. , 1. },500,50,-1.  + 1.1  I,300];//Timing

Altri esempi

Esempi di Julia sets  connessi:
    c = -1 ,   c = .3+.5I
Esempi di Julia sets  totalmente disconnessi (Cantor dust):
    c = .3,    c = I
Si riesce a disegnarli?

JP[{-2. , 2. , -1. , 1. }, 400, 100, -1. , 400]; // 
      Timing

JP[{-1.2, 1.2, -1.5, 1.5}, 400, 150,  .3 +  .5 I, 400]; // 
      Timing

JP[{-1, 1, -1.5, 1.5}, 900, 35,  .3, 400]; // 
      Timing

Esercizio

A cosa e` dovuto il fatto che in alcune figure si vedono come degli aloni,  come nel caso seguente?

JP[{-2. , 2. , -1. , 1. },300,40,-1.  +  .3  I,300];