Metodi iterativi per sistemi di equazioni lineari e non lineari

È dato il sistema lineare:

Ax = b

con A ϵ ℜ^{n x n} e x,b ϵ ℜⁿ con det(A) ≠ 0.

Per il teorema di Rouché-Capelli, poiché det(A) ≠ 0, ossia A non è singolare, esiste ed è unica la soluzione v ϵ ℜⁿ data da:

v = A^-1b

Tuttavia, poiché il calcolo della matrice inversa A^-1 è costoso e in generale malcondizionato, si cercano altri metodi per determinare la soluzione, i quali non richiedono il calcolo esplicito della matrice inversa A^-1.

Alla base dei metodi iterativi c’è l’idea do costruire un’opportuna successione di vettori x^(k) tale che:

Metodo di Jacobi

Per semplicità, sia con A ϵ ℜ^{3 x 3} con a_ii ≠ 0 per ogni i = 1,…,3 cioè il sistema lineare:

Tale sistema può essere riscritto come:

Tale riscrittura suggerisce direttamente la formulazione di un metodo iterativo nel modo seguente:

Ove si assume che venga fornito un vettore x⁽⁰⁾, detto vettore d’innesco, a partire dal quale viene generata la successione. Per l’i-esima equazione si può formalizzare in:

Per i = 1, …, n.

Il costo in termini di operazioni per calcolare una nuova iterazione è n². Tale costo dovrà essere a sua volta moltiplicato per il numero di iterazioni effettuate.

Infine, si tenga presente che se l’ipotesi a_ii ≠ 0 per ogni i non è verificata, si possono scambiare fra loro equazioni del sistema così che risulti soddisfatta.

Metodo di Gauss-Seidel

Una variante del metodo di Jacobi può essere semplicemente ottenuta utilizzando le componenti della soluzione già aggiornate. Più precisamente, nel caso n = 3 il metodo Gauss-Seidel è dato da:

Infatti, nella seconda equazione si utilizza il valore x₁^(k+1) appena calcolato al posto del valore x₁^(k) e nella terza equazione si utilizzano i valori x₁^(k+1) e x₂^(k+1) appena calcolati al posto di x₁^(k) e x₂^(k).

Per un sistema di generica dimensione n il metodo di Gauss-Seidel ha la seguente formulazione:

Per i = 1, …, n.

Come nel caso del metodo di Jacobi, il costo in termini di operazioni di tipo moltiplicativo per calcolare una nuova iterazione è n². Tale costo dovrà essere moltiplicato per il numero di iterazioni.

Se la matrice A è diagonalmente dominante in senso stretto, ossia se per ogni i = 1,…,n si ha:

Allora i metodi di Jacobi e Gauss-Seidel convergono.

Esempio:

Quindi:

La soluzione analitica è semplice e fornisce x = [1,2,3]^T.

La matrice è diagonalmente dominante, poiché per ogni elemento della diagonale, la somma degli altri due elementi di riga è minore.

Il seguente codice in C# trova la soluzione con il metodo iterativo di Gauss-Seidel:

using System;
namespace GaussSeidel
{
   class Program
   {
       /// <summary> Solve Ax=b, with the Gauss-Seildel iterative method
       /// static bool GaussSeidel(A, b, x, MaxIter=20, AbsError = 1e-15)</summary>
       /// <param name="double[,] A">A matrix</param>
       /// <param name="double[] x">b vector</param>
       /// <param name="double[] b">b vector</param>
       /// <param name="int MaxIter=20">Maximum number or iterations</param>
       /// <param name="double AbsError = 1e-15">Absolute error for convergence</param>
       /// <returns>bool convergence. True if the solution converges.</returns>
       static bool GaussSeidel(double[,] A, double[] x, double[] b, int MaxIter=20, double AbsError = 1e-15)
       {
           int n = b.Length;
           double[] y = new double[n], last = (double[])x.Clone();
           bool convergence;
           do
           {
               // Gauss-Seidel
               for(int i=0;i<n;++i)
               {
                 x[i] = (b[i] / A[i, i]);
                   for (int j = 0; j < n; ++j)
                   {
                       if (j != i)
                       {
                           x[i] -= ((A[i, j] / A[i, i]) * y[j]);
                           y[i] = x[i];
                       }
                   }
               }
               // check for convergence
               convergence = true;
               for (int i = 0; i < n; ++i)
               {
                   if (Math.Abs(x[i] - last[i]) > AbsError)
                   {
                       convergence = false;
                       break;
                   }
               }
               Array.Copy(x, last, n);
           } while (!convergence && --MaxIter > 0);
           return convergence;
       }
       static double[,] A = new double[,]
       {
           {4,-1,0},
           {-1,4,-1},
           {0,-1,4}
       };
       static double[] b = new double[] { 2, 4, 10 };
       static void Main(string[] args)
       {
           double[] x = { 1, 1, 1 }; // innesco
           bool cnv = GaussSeidel(A, x, b);
           if (cnv)
           {
               string s = string.Join(",", x);
               Console.WriteLine("La soluzione è: {0}", s);
           } else
           {
               Console.WriteLine("La soluzione non converge!");
           }
       }
   }
}

Applicazione ai sistemi non lineari

Facciamo un esempio concreto.

La legge di Darcy-Weisbach stabilisce che la differenza di pressione ai capi di un tubo lungo L metri, di diametro D vale:

Dove λ è il coefficiente d’attrito, ρ e w sono rispettivamente la densità e la velocità del fluido.

In un tubo liscio, il coefficiente d’attrito vale λ = 64 / Re se Re < 2000, mentre vale λ ≈ 0.316∙Re^-0.25 altrimenti (per Re < 10⁵).
Re è il numero di Reynolds:

Dove introduciamo la viscosità dinamica μ.

Se Re è minore o uguale a 2000 il moto è laminare, cioè dominato dalle forze viscose. Se invece è compreso tra 2000 e 4000 c’è una zona di transizione, dove il moto tende a mantenere le proprie caratteristiche, al di sopra è turbolento.

Se conosco la portata volumetrica Q di un fluido noto (ρ e μ noti) che attraversa il tubo, è dunque facile ricavare Δp, perché la velocità media è w = Q/A, quindi calcolo Re. A seconda che sia maggiore di 2000 oppure no, calcolo il coefficiente d’attrito λ appropriato e infine applico la legge di Darcy-Weisbach.

Un po’ meno ovvio è se conosco il Δp e mi chiedono la portata.

Il problema è che questa volta devo partire dalla legge di Darcy-Weisbach, ma non conosco ancora il coefficiente d’attrito λ, perché ancora non so dire se il moto è laminare oppure turbolento.

Anche in questo caso, un algoritmo iterativo ci può venire in soccorso.
Si tratta di ricalcolare iterativamente le tre equazioni, fino a quando si trova una convergenza o al massimo numero di tentativi.

Ed ecco il codice, sempre in C#:

using System;
namespace Darcy
{
    class Program
    {
        static double FlowRateFromdP(double dP, double L, double D,
            double rho, double mu, double lambda, int MaxIter = 100,
            double AbsError = 1e-15)
        {
            if (Math.Abs(dP) < AbsError) return 0.0;
            double sgn;
            if (dP < 0.0)
            {
                dP = -dP;
                sgn = -1.0;
            }
            else sgn = 1.0;
            double w = 0; // flow velocity
            double Re = 0; // Reynolds number
            double wLast = w, ReLast = Re,
                lambdaLast = lambda;
            bool convergence;
            do
            {
                if (lambda > 0.0)
                {
                    w = Math.Sqrt(2 * D * dP / (rho * L * lambda));
                }
                else w = 0.0;
                Re = rho * w * D / mu;
                if (Re<=64) lambda = 1;
                else if (Re <= 2000) lambda = 64 / Re;
                else lambda = 0.3164 / Math.Pow(Re, 0.25);
                convergence = (Math.Abs(w – wLast) < AbsError) &&
                    (Math.Abs(Re – ReLast) < AbsError) &&
                    (Math.Abs(lambda – lambdaLast) < AbsError);
                wLast = w;
                ReLast = Re;
                lambdaLast = lambda;
            } while (!convergence && –MaxIter > 0);
            return sgn*w * (0.25 * Math.PI * D * D);
        }
        static void Main(string[] args)
        {
            double dP = 7.85e+5;
            double L = 1.0;
            double D = 1.0e-2;
            double rho = 850.0;
            double mu = 0.051;
            double lambda = 0;
            double Q = FlowRateFromdP(dP, L, D, rho,
                mu, lambda, AbsError: 1e-6);
            Console.WriteLine(“Q = {0} m3/s”,Q);
            Console.ReadKey();
        }
    }
}