Exercice ma21¶

In [ ]:
import numpy as np
import numpy.linalg as lin
import matplotlib.pylab as plt

%matplotlib inline

np.set_printoptions(precision=3, linewidth=150, suppress=True)

On va augmenter le rayon de convergence la méthode Jacobi amméliorée faite en TD à savoir la méthode de Gauss-Seidel.

On étudiera sa convergence dans différents cas.

Gauss-Seidel¶

Lorsqu'on calcul le x suivant avec Jacobi on ne profite pas du fait que N est triangulaire et donc qu'on connait la nouvelle valeur de $x_n$ lorsqu'on calcule $x_{n-1}$. Avec Gauss-Seidel on utilise toujours la dernière valeur calculée ce qui accélère la convergence.

Pour résumer Gauss-Seidel d'un point de vu matriciel on a :

  • D = la matrice diagonale extraite de A : D = np.diag(np.diag(A))
  • L = la matrice stritecement triangulaire inférieure de A : L = np.tril(A, -1)
  • U = la matrice stritecement triangulaire supérieure de A : U = np.triu(A, 1)

et une itération est donnée par la formule suivante :

$$ (D + L)\, {\bf x}^{k+1} = -U\; {\bf x}^k + {\bf b} $$ ou $$ {\bf x}^{k+1} = D^{-1} \, ( -L\, {\bf x}^{k+1} - U\; {\bf x}^k + {\bf b}) $$ c.a.d. $$ \begin{bmatrix} x_{1}^{k+1} \\ x_{2}^{k+1} \\ \vdots \\ x_{n}^{k+1} \\ \end{bmatrix} = \begin{bmatrix} 1/a_{11} \quad 0 \quad \ldots \quad 0 \\ 0 \quad 1/a_{22} \quad \ldots \quad 0 \\ \vdots \\ 0 \quad 0 \quad \ldots \quad 1/a_{nn} \\ \end{bmatrix} \; \left( \; - \begin{bmatrix} 0 \quad 0 \quad \ldots \quad 0 \\ a_{21} \; 0 \quad \ldots \quad 0 \\ \vdots \\ a_{n1} \, a_{n2} \; \ldots \quad 0 \\ \end{bmatrix} \; \begin{bmatrix} x_{1}^{k+1} \\ x_{2}^{k+1} \\ \vdots \\ x_{n}^{k+1} \\ \end{bmatrix} - \begin{bmatrix} 0 \; a_{12} \; \ldots \; a_{1n} \\ 0 \quad 0 \; \ldots \; a_{2n} \\ \vdots \\ 0 \quad 0 \; \ldots \; 0 \\ \end{bmatrix} \; \begin{bmatrix} x_{1}^k \\ x_{2}^k \\ \vdots \\ x_{n}^k \\ \end{bmatrix} + \begin{bmatrix} b_{1} \\ b_{2} \\ \vdots \\ b_{n} \\ \end{bmatrix} \; \right) $$

Notons que je peux mettre $L\, {\bf x}^{k+1}$ à droite du signe égal si je résoud mon système ligne par ligne en commencant par le haut puisque dans ce cas les ${\bf x}^{k+1}$ utilisés sont connus. C'est ce qu'on a fait lors du dernier TP.

Surrelaxation de Gauss-Seidel¶

Comme on a fait avec Jacobi, on introduit de l'inertie avec $w$ :

$$ {\bf x}^{k+1} = w \, D^{-1} \, ( -L\, {\bf x}^{k+1} - U\; {\bf x}^k + {\bf b}) + (1-w) \; {\bf x}^k $$

Vérifiez que l'on arrive à l'écriture matricielle suivante :

$$ \left(\frac{D}{w} + L\right)\, {\bf x}^{k+1} = \left(\frac{1-w}{w} \, D - U\right)\; {\bf x}^k + {\bf b} $$

Écrit ainsi on voit que cette méthode consiste à avoir les éléments de la diagonale des 2 cotés de l'égalité. On peut interpréter cela comme un avantage lié à une meilleure répartition de l'information contenue dans la matrice (à tester pour savoir).

Programmons Gauss-Seidel surrelaxé¶

On écrira deux fonctions :

  • solve_triangular(L,b) qui retourne la solution de L x = b lorsque L est triangulaire inférieure
  • gauss_seidel_r(A, b, x0, w, n) qui fait n iteration de Gauss-Seidel en démarrant à x0 avec w le coefficient de relaxation donné en argument. Cette fonction retourne un tableau des valeurs de x calculées (donc tableau en 2D).

Comme toujours, attention à limiter les for et à faire le plus possible d'opérations vectorielles et matricielles.

In [ ]:
In [ ]:
In [ ]:
np.random.seed(123)
A = np.random.randint(10, size=(4,4))
b = A.sum(axis=1)
x0 = np.random.random(4)

res = gauss_seidel_r(A, b, x0, w=0.2, n=100)
print(res[-1])
In [ ]:
def plot_convergences(values, result):
    error = np.square(values - result).sum(axis = -1) / np.square(result).sum(axis=-1)
    error2 = np.square(np.diff(values)).sum(axis = -1) / np.square(values).sum(axis=-1)
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14,4))
    ax1.plot(range(len(error)), error)
    ax1.set_title('Erreur absolue normalisée')
    ax1.semilogy();
    ax2.plot(range(len(error2)), error2)
    ax2.set_title('Erreur relative normalisée')
    ax2.semilogy()
    print("Itération du minimum :",np.argmin(error), np.argmin(error2))
In [ ]:
plot_convergences(res, np.ones(4))

Est-ce que la méthode de Gauss-Seidel non relaxée converge dans ce cas ?

In [ ]:

Le bon cas¶

Trouver un seed qui permet de générer un cas qui ne converge pas avec Gauss-Seidel de base mais qui converge avec la relaxation ($w=0.2$).

In [ ]:

Tracer les courbes de convergence pour le cas retenu avec et sans relaxation.

In [ ]:
In [ ]:
In [ ]:

Étude de $w$¶

Toujours dans notre cas retenu, indiquer quel est l'intervale de valeurs de $w$ qui garantit la convergence pour notre système matriciel A x = b avec toujours le même x0 et un nombre d'itérations à déterminer.

Trouver la valeur optimiale de $w$ pour converger le plus rapidement pour ce cas.

La précision demandée pour l'intervale et la valeur optimale est de $10^{-2}$.

In [ ]:
In [ ]: