---

Nous avons continué la séance Python de janvier. Nous avons d'une part terminé le travail sur le Julia de exp(z-1.5) : j'ai mis regroupé les correspondantes dans celles de la 2e partie de séance de janvier, vers la fin. Puis nous avons regardé le Julia de sin(z)/2.

---

La fonction f(z)=sin(z)/2 a un point fixe 0 qui est attractif car f'(0)=1/2. Elle son unique valeur asymptotique est l'infini. Les points critiques de f (là où f' s'annule) sont les multiples entiers de π. Ils s'envoient sur 1/2 et -1/2, appelés valeurs critiques. Un théorème de Fatou assure qu'elles dans le bassin d'attraction de 0. Mais ici c'est facile à voir: comme |sin(x)|≤x, il suit que tout l'axe réel est dans le bassin d'attraction de 0.

En fait, c'est pour des valeurs élevées de Im(z) que sin(z) prend des valeurs de grand module (elle devient alors proche de exp(iz)/2i). D'autre part le sinus complexe préserve l'axe imaginaire où il est conjugué au sinus hyperbolique réel.

Premier essai avec une approche élémentaire : en blanc les points de départ dont l'orbite finit par tomber dans le disque de centre 0 et rayon 1/2 (contenu dans le bassin) ; en bleu quand l'orbite finit par passer trop loin de l'axe réel (auquel cas on ne pourrait plus itérer à cause de l'overflow) ; en rouge dans le cas éventuel où l'on atteint max_iter avant l'un des deux cas précédents.

[...]

def N2(z) : return z.real**2 + z.imag**2

def compute(j,i) :
  z = center + (j-width/2+(height/2-i)*1.0j)*xw/width
  k = 0
  reason = NONE
  while True:
    if k>=max_iter :    reason = MAX; break
    if abs(z.imag) > seuil : reason = SEUIL; break
    if N2(z) < 0.5 :    reason = BASIN; break
    z = cmath.sin(z)/2
    k += 1
  if reason==BASIN : return [255,255,255]
  if reason==SEUIL : return [0,0,255]
  return [255,0,0]

[...]

Voici le résultat brut :

Et avec du sur-échantillonnage (j'utilise comme précédemment la réduction bspline avec un logiciel qui, je crois, tient correctement compte du gamma=2.2) :

Note : Aucun pixel rouge, max_iter=50. Les zones rouges doivent être extrêmement petites.

La différence avec exp, c'est que l'ensemble de Julia est de mesure positive et non pas nulle, et il est de plus en plus dense quand Im(z) tend vers ±∞. Il y a un article de McMullen là-dessus : Area and Hausdorff dimension of Julia sets of entire functions dans [Transactions of the American Mathematical Society 300(1), mars 1987].

Cependant, J reste constitué de filaments (voir le travail de Devaney, Rempe, Rottenfusser, Rückert, Schleicher, ...). Comment les faire apparaître ? On peut "repérer" les filaments ainsi : en découpant le plan en morceaux et en les numérotant, on peut coder les orbites des points en fonction des numéro des morceaux visités. Pour un découpage bien choisi, un filament se caractérise comme l'ensemble des points de départ dont l'orbite a un code donné. Pour exp(z-1.5), on peut utiliser les bandes horizontales B_k: Im(z)∈[2kπ-π,2kπ+π].

On voit des grands paquets de filaments, qui se répètent par translation de π. On voit aussi leurs symétriques de l'autre côté de l'axe réel. L'image d'un paquet de filaments est soit l'ensemble des paquets situés en haut, soit ceux en bas. Un paquet du haut sur deux s'envoit sur ceux du haut. Et symétriquement en bas en utilisant f(−z) = −f(z). En fait on peut caractériser les filaments par un codage comme pour exp, avec les demi-bandes verticales

Ci-dessous, on colorie en fonction du signe de z quand il dépasse le seuil (donc la position au-dessus ou en-dessous de la dernière bande visitée).

def N2(z) : return z.real**2 + z.imag**2

def compute(j,i) :
  z = center + (j-width/2+(height/2-i)*1.0j)*xw/width
  k = 0
  reason = NONE
  while True:
    if k>=max_iter :     reason = MAX; break
    if z.imag > seuil :  reason = SEUIL_U; break
    if z.imag < -seuil : reason = SEUIL_D; break
    if N2(z) < 0.5 :     reason = BASIN; break
    z = cmath.sin(z)/2
    k += 1
  if reason==BASIN   : return [255,255,255]
  if reason==SEUIL_U : return [0,0,255]
  if reason==SEUIL_D : return [255,0,255]
  return [255,0,0]

Essayons maintenant d'adoucir les transitions au niveau où le nombre d'itérations diminue de 1, càd aux pré-images des lignes de seuil.

Voici une première méthode, qu'on ne pouvait pas appliquer pour exp(z-1.5).

Explications du code ci-dessous : on a attribué une couleur à un itinéraire fini donné. On se débrouille pour que les 2 couleurs alternées des doigts qui poussent sur un doigt donné aient pour demi-somme la couleur de ce dernier.

def N2(z) : return z.real**2 + z.imag**2

def compute(j,i) :
  z = center + (j-width/2+(height/2-i)*1.0j)*xw/width
  k = 0
  reason = NONE
  a = 0.25
  m = int(math.floor(z.real/pi+0.5))
  if m % 2 == 1 :
    a = 0.75
  b = 0.125
  while True:
    if k>=max_iter :     reason = MAX; break
    if z.imag > seuil :  reason = SEUIL_U; break
    if z.imag < -seuil : reason = SEUIL_D; break
    if N2(z) < 0.5 :     reason = BASIN; break
    z = cmath.sin(z)/mu
    if(abs(z.imag) < seuil) :
      m = int(math.floor(z.real/pi+0.5))
      if m % 2 == 0 :
        a = a-b
      else:
        a = a+b
    b = b/2
    k += 1
  if reason==BASIN   : return [255,255,255]
  if reason==SEUIL_U or reason==SEUIL_D :
    R=a*255.
    return [int(R),0,255]
  return [255,0,0]

On a choisi 1/4 et 3/4 pour la 1e génération, puis une structure en arbre du type:

Cela équivaut à faire la chose suivante: on attribue la valeur 0 aux blocs d'indice pair (qu'ils soient au dessus ou en dessous) et 1 aux blocs d'indice impair. Puis on écrit cette suite de valeurs pour l'itinéraire d'une orbite depuis le point de départ inclus jusqu'au moment où il dépasse le seuil, exclu. Cela donne une suite du type 01100101000101010100. On ajoute "0." devant, "1" derrière, et on considère cela comme un nombre binaire compris entre 0 et 1 : 0.011001010001010101001b. C'est le facteur de couleur (dans un espace linéaire d'intensités).

Défaut : faible contraste, les filaments proches ont une couleur proche.

Idée: inverser l'ordre des valeurs dans l'écriture binaire (plus besoin alors d'introduire un 1 au début/à la fin). Mais on perd la propriété de moyenne. On va pour compenser introduire une interpolation progressive le long d'un filament donné. Cette interpolation est analogue à ce qu'on avait dû faire pour exp(z-1.5).

mu=2.

seuil = 100.0
lseuil=math.log(seuil)
lls=math.log(lseuil)

def N2(z) : return z.real**2 + z.imag**2

def compute(j,i) :
  global mint,maxt
  z = center + (j-width/2+(height/2-i)*1.0j)*xw/width
  k = 0
  reason = NONE
  a = 0
  m = int(math.floor(z.real/pi+0.5))
  if m % 2 == 1 :
    a = 1
  b = a
  while True:
    if k>=max_iter :     reason = MAX; break
    if z.imag > seuil :  reason = SEUIL_U; break
    if z.imag < -seuil : reason = SEUIL_D; break
    if N2(z) < 0.5 :     reason = BASIN; break
    b = a
    m = int(math.floor(z.real/pi+0.5))
    if m % 2 == 0 :
      a = a/2
    else:
      a = (a+1)/2
    z = cmath.sin(z)/mu
    k += 1
  if reason==BASIN   : return [255,255,255]
  if reason==SEUIL_U or reason==SEUIL_D :
    t=(math.log(math.log(math.fabs(z.imag)))-lls)/(lseuil-lls)
    if t>1: t=1
    if t<0: t=0
    a=(1-t)*a+t*b
    R=a*255.
    return [int(R),0,255]
  return [255,0,0]

La couleur d'un itinéraire fini est donc calculée comme suit : On écrit l'itinéraire de l'orbite (toujours en excluant le seuil) par exemple 1111011011000011110. Ensuite on le renverse 0111100001101101111. Puis on ajoute "0." devant et on considère cette écriture comme du binaire: 0.0111100001101101111b. C'est le taux de couleur (en termes d'intensité).

En réalité j'ai légèrement modifié ça par coquetterie, en multipliant par deux la dernière décimale (si c'est un 0 alors ça ne change rien, si c'est un 1 alors il faut reporter la retenue). Dans notre exemple on obtient: 0.011110000110111000b.

Résultat :

Note: comme ici l'interpolation de couleur le long d'un filament est faite en linéaire dans mon code, j'ai donc sauvé l'image en indiquant gamma=1, et fait un rescale bspline qui en tient compte.

---