FEUILLE DE TP 7

---

Méthode de Newton (complexe...)

---

Manipulation d'images en Python

On commence par importer le module matplotlib.image, puis on importe l’image image0.png et on l’affiche

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.image as mpimg

img0 = mpimg.imread('image0.png')
img1 = mpimg.imread('image1.png')

fig = plt.figure(figsize=(16, 8))
for k, img in enumerate([img0, img1]):
    ax = fig.add_subplot(1, 2, k+1)
    ax.set_xticks([])
    ax.set_yticks([])
    ax.grid(False)
    ax.imshow(img)

Question 1

Quel est le type des objets img0 et img1 ? Vous pourrez utiliser la fonction type.

print(type(img0))
print(type(img1))

<class 'numpy.ndarray'>
<class 'numpy.ndarray'>

Question 2

Les objet img0 et img1 sont des numpy arrays de taille $(N_x, N_y, 4)$ où $N_x$ est le nombre de pixels (de points) dans la direction horizontale et $N_y$ est le nombre de pixels dans la direction verticale.

La couleur de chaque pixel est codée par un vecteur [R, G, B, A] de taille 4 dit RGBA (ou RVBA en français rouge-vert-bleu-alpha). Chaque coordonnée varie entre 0 et 1 (ou 0 et 255 selon la convention). [1, 0, 0, 1] code le rouge, [0, 1, 0, 1] le vert, [0, 0, 1, 1] le bleu, [1, 1, 1, 1] le blanc et [0, 0, 0, 1] le noir. Le dernier nombre code la transparence (0 complètement transparent et 1 opaque). Par exemple:

print(img0.shape)
print('La couleur du pixel de img0 de coordonnées (1,100) est donnée par [R G B A]=', img0[1,100,:])
print('La couleur du pixel de img1 de coordonnées (50,100) est donnée par [R G B A]=', img1[50,100,:])

(480, 640, 4)
La couleur du pixel de img0 de coordonnées (1,100) est donnée par [R G B A]= [0.9372549 0.9372549 0.9372549 1.       ]
La couleur du pixel de img1 de coordonnées (50,100) est donnée par [R G B A]= [0.  0.4 1.  1. ]

• Créez une image de taille $N \times N$ pixels entièrement noire, puis blanche, puis rouge.
• Créer une autre image de taille $N \times N$ dont la moitié droite est rouge et la moitié gauche est bleue.
• Essayez de jouer avec la transparence...

N = 16
liste_alpha = [1, 0.75, 0.5, 0.25, 0]

liste_couleur = [
    [0, 0, 0],  # noir
    [1, 1, 1],  # blanc
    [1, 0, 0],  # rouge
    [0, 1, 0],  # vert
    [0, 0, 1]   # bleu
]

for color in liste_couleur:
    fig = plt.figure(figsize=(len(liste_alpha), 1))

    for k, alpha in enumerate(liste_alpha):
        v = color + [alpha]
        img = np.empty((N, N, 4))
        img[..., :] = v
        ax = fig.add_subplot(1, len(liste_alpha), k+1)
        ax.set_xticks([])
        ax.set_yticks([])
        ax.grid(False)
        ax.imshow(img)

color1, color2 = [1, 0, 0, 1], [0, 0, 1, 1]
n = 5
N = 2**n
fig = plt.figure(figsize=(n, n))
for i in range(n):
    for j in range(n):
        ax = fig.add_subplot(n, n, n**2-n*i-j)
        img = np.empty((N, N, 4))
        for k in range(N):
            for l in range(N):
                if ((k >> (i+1)) % 2) ^ ((l >> (j+1)) % 2) :
                    img[k, l, :] = color1
                else:
                    img[k, l, :] = color2
        ax.set_xticks([])
        ax.set_yticks([])
        ax.grid(False)
        ax.imshow(img)

Fractale de Newton

On considère la fonction holomorphe $f:\mathbb{C} \rightarrow \mathbb{C}$ telle que $$ f(z) = z^3-1. $$

L'équation $f(z)= 0$ admets trois racines dinstinctes: $$ {r}_0 = -\frac{1}{2}+ \frac{\sqrt{3}}{2} i,\quad {r}_1 = -\frac{1}{2}- \frac{\sqrt{3}}{2} i\,,\quad r_2 = 1 . $$

La méthode de Newton peut-être appliquée sans modifications aussi dans ce cas, en utilisant la dérivée complexe $f'(z) = 3z^2$.

Ici nous voulons étudier les bassins d'attractions des différentes racines de $f$.

Question 3

Afin de définir le nombre complexe $z = x + i y$ en python, on définit z = complex(x, y) par exemple

z = complex(4,3)
print(f'Soit z={z}, alors z^2={z**2}')

Soit z=(4+3j), alors z^2=(7+24j)

Définissez la fonction $f(z) = z^3-1$ ainsi que sa dérivée $f'$ et un ndarray appelé roots de taille 3 contenant les trois racines distinctes de f.

def f(z): 
    return z**3 - 1

def df(z): 
    return 3*z**2
    
roots = np.array([
    complex(-1, np.sqrt(3))/2,
    -complex(1, np.sqrt(3))/2,
    1
])

Question 4

Implémentez une méthode de Newton modifiée Newton_mod à partir de la fonction Newton implémentée au TP précédent. On demande les modifications suivantes : on va utiliser le fait qu’on a une expression algébrique des racines de $f$ afin d’étudier plus simplement la dépendance de la méthode de Newton en fonction de la condition initiale et on remplace le critère d’arrêt par le critère d’arrêt: $$ |x_k − r_0| < \epsilon \quad \text{ OR }\quad |x_k − r_1| < \epsilon\quad\text{ OR }\quad |x_k − r_2| < \epsilon \quad\text{ OR }\quad k > N_{max}. $$ La fonction Newton_mod doit avoir en sortie un argument additionnel par rapport à la fonction Newton : un entier $b$ pouvant prendre 4 valeurs ($b=-1$ si la la méthode n’a pas convergé (c'est-à-dire si on dépasse le nombre d'iterations maximal), ou $b=j$ si la méthode a convergé vers la racine $r_j$ avec $j=0, 1, 2$).

def Newton_mod(f, df, x0, roots, epsilon, Nmax):
    """
    Algorithme de Newton pour résoudre f(x) = 0 
    
    Parameters
    ----------
    f: function
        fonction f
    df: function
        dérivé de f
    x0: float 
        approximation initiale
    roots: ndarray
        les zéros de la fonction f
    epsilon: float
        tolerance dans le critere d'arret
    Nmax: int
        nombre d'iterations maximal

    Returns
    -------
    niter: float
        nombre d'iteration
    xstar: float
        solution approchée
    xL: ndarray
        la liste des itérés de la méthode
    b: int
        l'indice du zéro s'il est trouvé (-1 sinon)
    """
    x = x0
    fx, dfx = f(x), df(x)
    niter = 0
    xL = [x0]
    while min(abs(x-roots)) >= epsilon and niter < Nmax:
        if abs(dfx) < 1.e-17:
            return -1, 0, x0
        x -= fx / dfx
        fx, dfx = f(x), df(x)
        niter += 1
        xL.append(x)
    b = -1
    for k, l in enumerate(roots):
        if abs(x - l) < epsilon:
            b = k
    return niter, x, np.asanyarray(xL), b

Question 5

• Définissez un vecteur $p$ de taille $N$, obtenu en divisant l'interval $[-1,1]$ en $N-1$ sous-intervalles égaux.
• Définissez une images noire de taille $N \times N$ comme précédemment et appelez la img0.

Chaque pixel d'img0 correspond à un point dans le carré $[-1,1]\times [-1,1]$.

• Pour chaque pixel correspondant aux indices $(i_x,i_y)$, appellez la méthode Newton_mod en prenant pour condition initiale le nombre complexe $$ z_0 = a+ i b , \text{ où } a= p[i_x] \text{ et } b= p[i_y] . $$
• Coloriez le pixel $(i_x,i_y)$ de img0 en noir si la méthode n’a pas convergé et sinon en rouge, vert ou bleu selon la racine vers laquelle elle a convergé.
• Essayez avec $N = 32$ puis $N = 512$.

N = 512
Nmax = 100
p = np.linspace(-1, 1, N)
img = np.zeros((N, N, 4))
img[..., -1] = 1

for i, pi in enumerate(p):
    for j, pj in enumerate(p):
        b = Newton_mod(f, df, complex(pi, pj), roots, 1.e-14, Nmax)[-1]
        if b >= 0:
            img[i, j, b] = 1

fig = plt.figure(figsize=(9, 9))
ax = fig.add_subplot(1, 1, 1)
ax.grid(False)
ax.set_xticks([])
ax.set_yticks([])
ax.imshow(img);

Question 6

Sur une nouvelle image, coloriez de manière similaire chaque pixel sur une échelle de gris proportionnellement au nombre d’itérations $n$ nécessaires à la convergence de la méthode, $n/N_{max} \in [0, 1]$; c'est-à-dire en saisissant le code [R G B] donné par $[n/N_{max},n/N_{max},n/N_{max}]$.

N = 512
Nmax = 100
p = np.linspace(-1, 1, N)
img = np.zeros((N, N, 3))

for i, pi in enumerate(p):
    for j, pj in enumerate(p):
        n, _, _, b = Newton_mod(f, df, complex(pi, pj), roots, 1.e-14, Nmax)
        if b >= 0:
            img[i, j, :] = n

nmax = img.max()
img /= nmax

fig = plt.figure(figsize=(9, 9))
ax = fig.add_subplot(1, 1, 1)
ax.grid(False)
ax.set_xticks([])
ax.set_yticks([])
ax.imshow(img);