CAMA : ma21 Surrelaxation pour Gauss-Seidel – Exercice

Cours du 11/05

import numpy as np
import scipy.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){\mathbf{x}}^{k+1}=-U{\mathbf{x}}^k+\mathbf{b}$$

ou

$${\mathbf{x}}^{k+1}=D^{-1}(-L{\mathbf{x}}^{k+1}-U{\mathbf{x}}^k+\mathbf{b})$$

c.a.d.

$$\left[\begin{array}{c}{x}_1^{k+1}\\ x_2^{k+1}\\ ⋮\\ x_n^{k+1}\end{array}\right]=\left[\begin{array}{c}1/a_{11}0\dots 0\\ 01/a_{22}\dots 0\\ ⋮\\ 00\dots 1/a_{nn}\end{array}\right](-\left[\begin{array}{c}00\dots 0\\ a_{21}0\dots 0\\ ⋮\\ a_{n1}{a}_{n2}\dots 0\end{array}\right]\left[\begin{array}{c}{x}_1^{k+1}\\ x_2^{k+1}\\ ⋮\\ x_n^{k+1}\end{array}\right]-\left[\begin{array}{c}0a_{12}\dots a_{1n}\\ 00\dots a_{2n}\\ ⋮\\ 00\dots 0\end{array}\right]\left[\begin{array}{c}{x}_1^k\\ x_2^k\\ ⋮\\ x_n^k\end{array}\right]+\left[\begin{array}{c}{b}_1\\ b_2\\ ⋮\\ b_n\end{array}\right])$$

Notons que je peux mettre

$L{\mathbf{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

${\mathbf{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$ :

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

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

$$(\frac{D}{w}+L){\mathbf{x}}^{k+1}=(\frac{1-w}{w}D-U){\mathbf{x}}^k+\mathbf{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é à un gain d'information lors des opérations matrice vecteur.

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 avec w donné en argument à partir
  de x0.
  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.

Solution

def solve_triangular(L, b):
    x = np.empty(len(b))
    x[0] = b[0] / L[0,0]
    for i in range(1,len(L)):
        x[i] = (b[i] - L[i,:i] @ x[:i]) / L[i,i]
    return x
            
# je teste            
A = np.tril(np.random.randint(10, size=(4,4)))
b = A.sum(axis=1)
solve_triangular(A,b)

array([1., 1., 1., 1.])

def gauss_seidel_r(A, b, x0, w=0.5, n=100):
    D = np.diag(np.diag(A))
    L = np.tril(A, -1)
    U = np.triu(A, 1)
    L = D / w + L
    U = ((1-w) / w) * D - U
    x = x0
    values = []
    for _ in range(n):
        x = solve_triangular(L, U @ x + b)
        values.append(x)
    return np.array(values)

L'algorithme de Wikipedia est difficile alire, lent et il s'agit de Jacobi avec relaxation et non Gauss-Seidel avec relaxation, il est tout ce qu'il ne faut pas faire.

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])

[1. 1. 1. 1.]

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))

plot_convergences(res, np.ones(4))

Resultat

Itération du minimum : 99 99

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

Reponse

gauss_seidel_r(A, b, x0, w=1, n=100)[-1]  # oui

array([1., 1., 1., 1.])

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$).

Solution

seed = 0
while True:
    np.random.seed(seed)
    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)
    res2 = gauss_seidel_r(A, b, x0, w=1, n=100)
    if np.square(res[-1] - np.ones(4)).sum() < 0.01 and np.square(res2[-1] - np.ones(4)).sum() > 1 :
        print(seed)
        break
    seed += 1

87

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

Solution

np.random.seed(87)
A = np.random.randint(10, size=(4,4))
b = A.sum(axis=1)
x0 = np.random.random(4)

plot_convergences(gauss_seidel_r(A, b, x0, w=0.2, n=100), np.ones(4))

Itération du minimum : 95 97

plot_convergences(gauss_seidel_r(A, b, x0, w=1, n=100), np.ones(4))

Itération du minimum : 0 0

É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}$.

Solution

# utilisons une méthode par dichotomie
wmin = 0
wmax = 1
while wmax - wmin > 1E-2:
    w = (wmax + wmin) / 2
    res = gauss_seidel_r(A,b,x0, w, n=1000)[-1]
    if np.square(res - np.ones(4)).sum(axis=-1) < 1E-3:  # converge
        wmin = w
    else:  
        wmax = w
print(f"Intervale de convergence : ]0,{(wmin+wmax)/2}]")

Intervale de convergence : ]0,0.69140625]

plot_convergences(gauss_seidel_r(A, b, x0, w=0.7, n=500), np.ones(4)) # ca converge tres doucement

Itération du minimum : 493 494

plot_convergences(gauss_seidel_r(A, b, x0, w=0.7 + 1E-2, n=500), np.ones(4)) # ca diverge clairement

Il semble que

$\omega =0.7$ soit la limite de la convergence.

# Méthode brutale (Newton serait plus joli mais je ne sais pas si vous connaissez)
# de toute facon c'est très rapide

N = 500
best_w = 0
best_it = N

for i in np.arange(0.01, 0.7, 0.01):
    it_cv = np.argmax(np.square(gauss_seidel_r(A, b, x0, i, 500) - np.ones(4)).sum(axis=-1) < 1E-6) 
    # attention, si la réponse est 0 cela veut dire que cela n'est pas descendu en dessous de 1E-6
    if it_cv > 0 and it_cv < best_it:
        best_w = i
        best_it = it_cv
best_w, best_it

(0.36000000000000004, 154)

plot_convergences(gauss_seidel_r(A, b, x0, w=0.36, n=100), np.ones(4))

Itération du minimum : 99 99

Trouver la valeur optimale de

$\omega$ doit bien sûr pouvoir être fait rapidement. Pour certains problèmes particulier on connait la formule qui donne le

$\omega$ optimal, sinon il faut utiliser des heuristiques sans garanties.