Pré-requis

Plutôt que de réinventer la roue, on utilise le module de calcul formel sympy qui nous permettra d'avoir accès aux nombres complexes.

Dans un terminal, tapez ce qui suit :

  wget http://sympy.googlecode.com/files/sympy-0.6.5.tar.gz
  tar zxvf sympy-0.6.5.tar.gz
  cd sympy-0.6.5

Ensuite, nous vous conseillons de créer un répertoire dans votre home, au nom de python par exemple, qui contiendra l'ensemble des modules, non-installés à l'IUT, que vous utiliserez :

  mkdir ~/python

Reste alors à y copier la bibliothèque sympy :

  cp -R sympy ~/python

Reste à tester le tout...

  cd ~/python
  python
  >>> from sympy import *
  >>> I
  I
  >>> I**2
  -1

On a donc bien notre imaginaire pur !

Chez vous, vous pouvez installer sympy, et tout autre module, à sa place normale dans votre système, en tapant la commande (dans le répertoire sympy-0.6.5) :

  sudo python setup.py install

Liens sympy :

• Site.
• Documentation.

Premier contact avec les complexes

Dans sympy, l'imaginaire pur (dont le carré est égal à -1) est noté I.

Opérations élémentaires

Addition, soustraction

Définissons deux complexes a et b, et montrons comment réaliser les opérations élémentaires avec :

  >>> a=2+3*I
  >>> b=5-2*I

  >>> a+b
  7 + I

  >>> a-b
  -3 + 5*I

Multiplication, puissance

  >>> a*b
  (2 + 3*I)*(5 - 2*I)

Demandons à sympy de développer l'expression :

  >>> expand(a*b)
  16 + 11*I

La puissance ne pose aucun problème :

  >>> a
  2 + 3*I
  >>> a**5
  (2 + 3*I)**5
  >>> expand(a**5)
  122 - 597*I

Quotient

Passons au quotient :

  >>> a/b
  (2 + 3*I)/(5 - 2*I)

On peut en demander la simplification :

  >>> simplify(a/b)
  4/29 + 19*I/29

Conjugaison

On peut enfin obtenir la partie conjuguée...

  >>> a
  2 + 3*I
  >>> a.conjugate()
  2 - 3*I

Module

On peut calculer le module $\sqrt{x^2+y^2}$ du complexe $x+iy$, avec la fonction abs :

  >>> b=2+3*I
  >>> abs(b)
  13**(1/2)

Ce qui donne bien la valeur souhaitée. On rappelle que la puissance 1/2 est la racine carrée ; pour s'en convaincre :

  >>> pprint(abs(b))
    ⎽⎽⎽⎽
  ╲╱ 13 

Equation du second degré

L'une des premières utilisations des complexes est la résolution des équations du second degré, comme $x^2+2x+3$ :

  >>> x = Symbol('x')
  >>> solve(x**2+2*x+3,x)
  [-1 + I*2**(1/2), -1 - I*2**(1/2)]

On rappelle en effet que si une équation du second degré $ax^2+bx+c=0$ a un discriminant $\Delta = b^2-4ac$...

• strictement positif, alors les racines sont réelles : $\frac{-b\pm\sqrt{\Delta}}{2a}$
• nul, alors on a une racine double : $\frac{-b}{2a}$
• strictement négatif, alors les racines sont complexes conjuguées : $\frac{-b\pm i\sqrt{|\Delta|}}{2a}$

Travaux pratiques

Créer votre propre fonction, qui reproduit le comportement de la fonction solve :

1. Demandez à l'utilisateur de saisir son équation du second degré.
2. Calculez le déterminant.
3. Retournez la, ou les racines, éventuellement complexes.

Application : M de Mandelbrot

Le but de ce qui suit est de représenter la fractale de Mandelbrot.

On rappelle qu'au point du plan cartésien de coordonnées $(x,y)$, on peut associer le complexe $x+iy$.

Attach:dd.png Δ

Présentation de la fractale

Pour chaque point $c = x+iy$ du plan complexe, on regarde le nombre d'itérations nécessaires pour que la suite récurrente

  $z_0     = 0$
  $z_{n+1} = z_n^2+c$

dépasse, en module, une certaine limite fixée. La couleur du pixel en position $c$ dépendra de ce nombre d'itérations.

Travaux Pratiques : détermination du nombre borné d'itérations

• Créer la fonction $nombreIt1(c,borne)$ telle que
  • le premier paramètre $c$ est un nombre complexe
  • le second paramètre $borne$ est un réel positif
  • la fonction retourne le plus petit nombre $n$ tel que le module de $z_{n}$ est supérieur à la $borne$. Pour cela, on utilisera la méthode $evalf()$ qui, appliquée à une expression sympy formelle, retourne sa valeur flottante comme dans l'exemple suivant:

  >>> from sympy import *

  >>> abs(a)
  34**(1/2)
  >>> abs(a).evalf()
  5.83095189484530

• Vérifier la fonction avec $nombreIt1(0.5,2)$, $nombreIt1(0.8*I,2)$ et $nombreIt1(0.5*I,2)$
• Pour éviter le dernier problème, construire la fonction $nombreIt(c,borne,itermax)$ telle que
  • les mêmes autres contraintes que $nombreIt1(c,borne)$ mais
  • itermax est le nombre que la fonction retournera si tous les $z_i$ sont de module plus petit que $c$ avec $1 \le i \le n$

Plus $itermax$ sera grand, plus le découpage du M de Mandelbrot sera précis, mais plus long sera le temps de calcul.

Tracé de l'image

On utilise le module Image. On commence par créer une nouvelle image M, en niveaux de gris (mode 'L'), de taille $32\times 32$ ;

  >>> From PIL import Image 
  >>> monImage = Image.new('L',(32,32))
  >>> monImage.show(command = 'eog')

On parcours alors chaque pixel $(k,l)$ de l'image, en regardant à chaque fois le nombre d'itérations nécessaire à la suite $z_0 = 0$, $z_{n+1} = z_n^2+c$ ($c=k+il$) pour dépasser, en module, 2.

Ce nombre d'itération, tronqué à la valeur 20, est le niveau de gris du pixel $(k,l)$...

  >>> for k in range(32):
  ...     for l in range(32):
  ...         couleur = nombreIt(k/32.+I*l/32.,2,20)
  ...         monImage.putpixel((k,l),couleur)
  ...     print k,
  ... 

On a donc réalisé le M de Mandelbrot. Seulement, chaque couleur de pixel va de 0 à 20, avec 256 niveaux de gris possibles. On va donc faire une homothétie sur ces niveaux de gris : les multiplier par $\frac{256}{20}$, et afficher le résultat...

Travaux pratiques : affichage complet

• Construire le Mandelbrot de 64 sur 64 en 20 niveaux de gris, 40 niveaux de gris...

Problème du temps de calcul

Sympy est pour l'instant trop lent pour exécuter une telle quantité de calculs. Dans ce qui suit, on ne fait plus des calculs complexes, mais on utilise les coordonnées des pixels:

  def nombreIt(c,borne,itermax):
      n,z=0,[0,0]
      while z[0]**2+z[1]**2<borne**2 and n<itermax:
          z=[z[0]**2-z[1]**2+c[0],2*z[0]*z[1]+c[1]]
          n+=1
      return n
  taille=512
  monImage = Image.new('RGB',(taille,taille))

  for k in range(taille):
      for l in range(taille):
          couleur = nombreIt([k/float(taille),l/float(taille)],2,50)/50.*256
          monImage.putpixel((k,l),(couleur,0,0))
  monImage.show(command = 'eog')

On en a profité pour changer la couleur et la taille de l'image, et pour augmenter le nombre d'itérations (50) pour améliorer l'approximation :

Travaux Pratiques : recentrer l'image

Il reste un point à améliorer : on ne voit qu'une partie du M. Modifier votre code pour obtenir :