Successive Over Relaxation

In analisi numerica, il metodo Successive over-relaxation, comunemente detto SOR o metodo del rilassamento successivo, è una variante del metodo di Gauss-Seidel per la risoluzione di sistemi di equazioni lineari, caratterizzato da una convergenza più rapida. Fu introdotto dal matematico statunitense David M. Young Jr. nel 1950.

Formulazione

Si consideri un sistema di equazioni lineari di incognita :

dove

La matrice può essere decomposta come dove è la componente diagonale, e sono le componenti triangolari superiori e inferiori strette. Introducento un parametro di rilassamento nel metodo di Gauss-Seidel, si ottiene:

per La formulazione vettoriale è:

Convergenza

La scelta del parametro di rilassamento non è un problema semplice e dipende dalle proprietà della matrice dei coefficienti. Nel 1947, Ostrowski ha dimostrato che, se è simmetrica e definita positiva, il metodo SOR converge se e solo se

Raggio spettrale della matrice di iterazione del metodo SOR. Il grafico mostra la dipendenza dal raggio spettrale della matrice di iterazione del metodo di Jacobi .

Inoltre, nel caso in cui la matrice dei coefficienti sia anche tridiagonale, è possibile determinare un parametro ottimale in corrispondenza del quale la velocità di convergenza sia la più elevata possibile:

dove è la matrice di iterazione del metodo di Jacobi e il suo raggio spettrale. Inoltre, in corrispondenza di tale valore il fattore asintotico di convergenza è pari a

Algoritmo

Esempio con Python e NumPy

Di seguito viene riportato un esempio di implementazione in Python del metodo SOR:

import numpy as np

def successive_over_relaxation(A, b, x0, kmax, tol, omega):
    n = A.shape[0]
    x_old = x0.astype(float).copy()
    x = x_old.copy()

    k = 0
    err = tol + 1

    while err > tol:

        if k >= kmax:
            raise RuntimeError("Numero massimo di iterazioni raggiunto.")

        for i in range(n):

            if A[i, i] == 0:
                raise ValueError(f"Elemento diagonale nullo nella riga {i}")

            s = b[i]

            for j in range(i):
                s -= A[i, j] * x[j]

            for j in range(i + 1, n):
                s -= A[i, j] * x_old[j]

            x[i] = omega * (s / A[i, i]) + (1 - omega) * x_old[i]

        err = np.linalg.norm(x - x_old, ord=np.inf)
        x_old = x.copy()
        k += 1

    return x, k

Voci correlate

Collegamenti esterni