Programmation Fonctionnelle, TP1

Version & licenses
Creative Commons License

Programmation Fonctionnelle, TP1 : Ensemble de Mandelbrot.

Guyslain Naves

L'objectif de ce travail est de réaliser des dessins de fractales, en utilisant le même moteur de calcul sur plusieurs polynômes complexes. Ceci nous amènera à définir des fonctions simples, à apprendre à faire des sortes de boucles à partir de récursion, et aussi à écrire quelques fonctions d'ordre supérieur.

Image without legend

Nous rappelons qu'un nombre complexe est de la forme $c = a + i b$, où $a$ est la partie réelle de $c$ et $b$ sa partie imaginaire. Nous pouvons aussi le représenter comme un vecteur du plan 2-dimensionel $(a,b)$. Nous coderons les nombres complexes comme des couples de flottants, ainsi, $c = 0.123 + 4.56 i$ sera encodé par la valeur (notez comment sont construits les couples)

  1. let c = (0.123, 4.56)

Vous utiliserez les fonctions du module Pervasive de la librairie standard pour manipuler les nombres flottants et accéder aux fonctions mathématiques standards. L'API des modules est disponible sur cette page : la libraire standard (ajoutez l'adresse à vos favoris).

Mise en place

Consulter le sommaire du cours pour mettre en place votre environnement de travail. Vous trouverez aussi une aide contenant des exemples de syntaxe.

Avant de vous lancer dans le TP, prenez le temps de découvrir l'éditeur que vous utilisez et l'interprèteur Ocaml. Vous pouvez reprendre les exercices du TD et les compiler par exemple pour vous familiariser avec ce nouvel environnement.

Convergence d'une suite récursive complexe

Question : Algèbre des complexes

L'addition et la multiplication complexes sont donnés par $(a + i b) + (c + i d) = (a+c) + i (b + d)$ et $(a + i b) \times (c + i d) = ac - bd + i (bc + ad)$. Utilisez cela pour coder l'addition, la multiplication et la puissance par un entier positif de nombres compexes en Ocaml.

Pour la puissance, utilisez une définition récursive, en utilisant le mot-clé rec.

Utilisez la boucle d'interaction pour tester vos fonctions.

Ensuite, vous pouvez définir des opérateurs infixes pour utiliser ces fonctions :

  1. let ( +! ) = add
  2. let ( *! ) = multiply
  3. let ( ^! ) = power
Question : Module d'un complexe

Le module d'un complexe est sa distance à l'origine dans le plan complexe. Si $c = a + i b$, alors le module de $c$ est $\left|c\right| = \sqrt{a^2 + b^2}$. Coder la fonction module d'un nombre complexe.

Pour la racine carrée, la fonction idoine est définie dans la librairie standard, dont vous pouvez voir la documentation ici.

Calcul d'une suite récursive

Une suite récursive est une séquence infinie de valeurs définies par une relation de récurrence, par exemple par $u_{n+1} = f(u_{n}), u_0 = c$. Étant donné la fonction $f$ et la constante $c$, nous voulons calculer les valeurs de cette suite.

Calculer la valeur suivante

Nous représentons la $n$e valeur de la suite par une paire $(u_n,n)$. Écrire une fonction next qui étant donné $f$ et la $n$e valeur, calcule la $(n+1)$e valeur.

  1. val next : ('value -> 'value) -> ('value * int) -> ('value * int)

  2. let res = next (fun x -> 2*x) (16,4)
  3. val res : int * int = (32,5)
Question : Répéter un calcul $n$ fois

Nous allons écrire une version fonctionnelle de la boucle for. Le principe d'un for est de faire $n$ fois la même opération. Étant données une fonction $f$ et une valeur $v$, nous voulons donc calculer $f (f (f \ldots (f v)\ldots ))$, avec $n$ appels à $f$.

Définissez une fonction récursive repeat pour effectuer ce calcul.

  1. val repeat : int -> f:('value -> 'value) -> 'value -> 'value

Ainsi, pour calculer le $n$e terme, il suffit de faire :

  1. let u_0 = 1
  2. let (u_100,_) = repeat 100 ~f:(next (fun x -> x+3)) (u_0,0)

  3. val u_100 : int = 301

Utilisez repeat pour afficher "Hello world!" dix fois.

Question : Répéter un calcul jusqu'à une condition

Passons à l'équivalent fonctionnel du while. En s'inspirant de repeat, écrivez une fonction qui applique $f$ jusqu'à ce que le résultat soit supérieur ou égal à 100 :

  1. val repeat_until_at_least_100 : f:(int -> int) -> int -> int

Ensuite, remplacez la condition au moins 100 par un prédicat sur la valeur $v$ que la fonction prend en argument :

  1. val repeat_until :
  2. predicat:('value -> bool) ->
  3. f:('value -> 'value) ->
  4. 'value -> 'value

On peut alors redéfinir :

  1. let repeat_until_at_least_100 = repeat_until ~predicat:(fun n -> n >= 100)
Question : Exemple d'évaluation d'une suite

Soit la suite d'entiers définie par la fonction collatz :

  1. let collatz i =
  2. if i mod 2 = 0 then i/2
  3. else 3*i + 1
  4. let collatz_0 = 13

Calculez à l'aide de repeat_until le plus petit indice $n$ tel que le $n$e terme de la suite de Collatz soit 1.

Indice de divergence

On considère, pour un complexe donné $c$, le polynôme complexe $P_c(z) = z^2 + c$. L'ensemble de Mandelbrot est défini comme l'ensemble des complexes $c$ tels que la suite $u_n$ générée par $P_c$ avec $u_0 = 0$ converge. On le voit en noir sur l'image suivante :

Image without legend

Savoir si $u_n$ converge ou diverge est un problème difficile, nous nous contentons d'utiliser une heuristique simple. Nous fixons un rayon $r \geq 2$ réel (par exemple, $r = 10^{10}$), et une borne entière $B$ (par exemple, $B=100$). Si $\left| P^{(n)}_c(0) \right| \geq r$ pour un $n \leq B$, alors $u_n$ diverge certainement. Sinon, on supposera (peut-être à tort) que $u_n$ converge. Autrement dit, nous calculons au plus $B$ itérations de la suite, et dès que nous obtenons un module supérieur à $r$, nous savons que la suite diverge.

$r = 2$ suffit pour être sûr que la suite diverge. En prenant une grande valeur de $r$, nous obtenons une meilleure précision sur la vitesse à laquelle la suite diverge (au prix d'un temps de calcul plus long).

Écrivez une fonction utilisant repeat_until, qui étant donné $r$, $B$ et $P_c$, détermine le plus petit $n$ tel que $\left| P^{(n)}_c(0) \right| \geq r$ ou $n > B$. On traduit ensuite cet entier en flottant (via la fonction float_of_int), qui est négatif si $n > B$.

  1. val diverging_term :
  2. radius:float ->
  3. bound:int ->
  4. f:(float * float -> float * float) ->
  5. float * float -> (float * float) * int

  6. val divergence :
  7. radius:float ->
  8. bound:int ->
  9. f:(float * float -> float * float)
  10. -> float * float -> float

La valeur du nombre d'itérations effectuées sera plus tard traduite en couleur pour donner la couleur du point de coordonnée $c$.

Rendu graphique

On utilise la librairie Graphics, que l'on charge avec le code suivant :

  1. #load "graphics.cma";;
  2. open Graphics
  3. let _ = open_graph ""
  4. let _ = resize_window 400 400

Si tout s'est passé correctement, une fenêtre blanche devrait s'être ouverte, dans laquelle nous allons dessiner.

Question : Itération sur un intervalle.

Écrivez une fonction récursive iter_interval, telle que iter_interval ~f ~min ~max applique $f$ sur tous les entiers de l'intervalle $[\text{min},\text{max}]$.

  1. val iter_interval : f:(int -> unit) -> min:int -> max:int -> 'a

Nous disposons donc maintenant d'une autre version de la boucle for, qui permet d'exécuter un calcul sur tous les éléments d'un intervalle. Par exemple,

  1. let () = iter_interval ~f:(fun n -> Printf.printf "%d\n" n) ~min:1 ~max:100

provoque à l'exécution l'impression en sortie standard des 100 premiers entiers strictement positifs.

Question : Itération sur un rectangle d'entier.

En déduire une fonction iter_rectangle telle que iter_rectangle ~f ~xmin ~xmax ~ymin ~ymax applique f x y pour tout x $\in$ [xmin,xmax] et tout y $\in$ [ymin, ymax].

  1. val iter_rectangle :
  2. f:(int -> int -> unit) ->
  3. xmin:int -> xmax:int ->
  4. ymin:int -> ymax:int -> 'b =
  5. <fun>

Ceci nous donne donc une sorte de boucle for en deux dimensions. Utilisez-la pour calculer les tables de multiplications jusqu'à 9.

Question : Affichage d'une fonction sur le plan.

À un nombre complexe $a + i b$, on associe le point du plan $(a,b)$.

Complétez la fonction suivante, qui affiche une image selon les arguments :

  • color_of_point : (float *. float) -> color donne la couleur associée à un nombre complexe (donc un point du plan),
  • center est le point au centre de l'image affichée,
  • width et height sont les dimensions du rectangle du plan à afficher.

Nous avons déjà écrit la fonction qui transforme un pixel en une couleur, et celle qui remplit un pixel. Il suffit donc d'appeller draw_pixel pour tous les pixels de la fenêtre graphique.

  1. let display ~color_of_point ~center ~width ~height =
  2. let xmin = 0 and ymin = 0 in
  3. let xmax = Graphics.size_x () and ymax = Graphics.size_y () in

  4. let to_complex (x,y) =
  5. (2. *. float x /. float xmax -. 1.,
  6. 2. *. float y /. float ymax -. 1.
  7. )
  8. in
  9. let rescale (re,im) = (re *. width /. 2., im *. height /. 2.) in
  10. let translate complex = complex +! center in
  11. let set_color_of_pixel pixel =
  12. pixel |> to_complex |> rescale |> translate |> color_of_point |> set_color
  13. in

  14. let draw_pixel x y =
  15. set_color_of_pixel (x,y);
  16. Graphics.plot x y
  17. in
  18. ... (* TODO *)

Testez votre code avec cet appel. Le résultat obtenu est-il celui auquel vous vous attendiez ?

  1. let _ =
  2. display
  3. ~color_of_point:(fun (x,y) -> rgb (int_of_float x) (int_of_float y) 0)
  4. ~center:(127.5,127.5)
  5. ~width:256. ~height:256.
Question : calcul de la couleur d'un point.

Étant donné $P_c$, la couleur attribuée au point $c$ est le noir si l'indice de divergence est supérieur à la borne sur le nombre d'itération. Sinon, la couleur attribuée est fonction du nombre d'itération. Les choix possibles de couleurs sont illimités. Une idée générale est d'attribuer deux couleurs d'autant plus proches que les indices de divergence sont eux-mêmes proches, afin d'obtenir des gradients de couleurs. Nous allons coder des gradients de couleurs. On commence par faire un dessin en noir et blanc.

Question : une première image en noir et blanc

Définissez le polynome complexe $P_c(z) = z^2 + c$ comme une fonction de $c$ et $z$. Puis codez une fonction qui étant donné $c$, calcule l'indice de divergence de $P_c$ en appelan la fonction divergence et retourne une couleur, blanc ou noir, selon que le nombre d'itérations dépasse ou non la borne. Enfin, affichez le résultat avec la fonction display :

  1. val color_of_point_black_and_white : (float * float) -> color
  2. let _ =
  3. display
  4. ~color_of_point:color_of_point_black_and_white
  5. ~center:(-1.,0.)
  6. ~width:5. ~height:5.

Amélioration de l'image

Question : mixer des couleurs

On représente les couleurs par des triplets de flottants compris entre 0 et 1 en notation RGB. On souhaite faire un mélange de couleurs avec différentes proportions. Par exemple 2 unités de bleu avec 3 unités de jaune. Il s'agit donc de calculer une moyenne pondérée de deux couleurs, un barycentre.

Nous utilisons la fonction élémentaire suivante pour calculer le barycentre de deux flottants :

  1. let barycenter coef1 coef2 x1 x2 =
  2. (x1 *. coef1 +. x2 *. coef2) /. (coef1 +. coef2)

Nous pouvons ensuite l'utiliser pour écrire une fonction pour les triplets de flottants. De façon générale, il est toujours possible de prendre un opérateur sur les flottants et d'en faire un opérateur sur les triplets de flottants. C'est ce que fait la fonction suivante :

  1. let lift_to_triple ~op (l1,m1,r1) (l2,m2,r2) =
  2. (op l1 l2, op m1 m2, op r1 r2)

Utiliser ces deux fonctions pour coder le mélange de deux couleurs :

  1. val mix :
  2. (float * float * float) -> float ->
  3. (float * float * float) -> float ->
  4. (float * float * float) = <fun>
Question : Créer un gradient

Nous devons transformer un flottant en couleur. Nous fixons deux flottants min et max et deux couleurs associées, ce sont les limites du gradient. Si la valeur value dont on cherche la couleur est entre min et max, on calcule le barycentre des deux couleurs avec les coefficients max -. value et value -. min, et cela crée un gradient linéaire.

Si value n'est pas dans les bornes, on utilise une autre fonction appelée otherwise pour calculer la couleur associée.

Implémentez cette idée en complétant la fonction suivante :

  1. let gradient ~min ~max ~color_min ~color_max ~otherwise value =
  2. ... (* TODO *)

Créez un gradient simple et ajoutez des couleurs à l'image précédemment en noir et blanc.

Question : Continuiser la valeur de divergence
Image without legend

La valeur de divergence étant entière, les images font apparaitre des zones uniformes de couleurs qui nuisent à leur aspect esthétique. Cette partie propose un amélioration pour lisser les couleurs. Pour cela, il faut raffiner le calcul de cette valeur. Notons $d$ le degré du polynome $P_c$ utilisé (pour l'instant, $d = 2$). Utilisez la formule suivante :

$$ n - \frac{\log \frac{\log |z_n|}{\log r}}{\log d}$$

avec :

  • $(z_n,n)$ le $n$e terme de la suite,
  • $r$ le rayon choisi,
  • $d = 2$.
Image without legend
Question : Combiner des gradients

En utilisant un seul gradient linéaire entre $0$ et $100$, nous obtenons une image assez décevante parce que beaucoup plus de points sont proches de $0$ que de $100$, et donc la linéarité du gradient est contre-effective.

De plus, une seul gradient ne permet pas d'avoir d'aussi jolies couleurs qu'avec trois ou quatre gradients qui se partagent l'ensemble des flottants possibles.

Créez des schémas de couleurs utilisant plusieurs gradients se partageant l'ensemble des flottants compris entre $0$ et $100$.

Quelques points intéressants sur lesquels zoomer :

  • $-0.80 + 0 \cdot i$ avec largeur et hauteur de $4$,
  • $-0.82 - 0.19 \cdot i$ avec largeur et hauteur de $0.038$,
  • $0.001643721971153 + 0.822467633298876 \cdot i$ avec largeur et hauteur $0.0000003$,
  • $-0.743643887037151 + 0.13182590420533 \cdot i$ avec des largeurs de plus en plus faibles.
Image without legend
Image without legend
Image without legend

À vous de jouer...

Essayer de tracer d'autres fractales, en utilisant des polynomes différents, par exemple $z \rightarrow z^3 + c$, en testant des fonctions de colorations plus subtiles et en zoomant sur des points intéressants du plan.