[C++] Algoritmo di eliminazione di Gauss

[oleg.fresi]

Stavo studiando l'algoritmo di eliminazione di Gauss, solo che non avevo la più pallida idea di come implementarlo. Ho trovato in rete questo codice:

int main()
{
	int i, j, k, n; 

	float a[10][10], b, x[10];
	cout << "Inserisci la dimensione della matrice" << endl;
	cin >> n;
	cout << "inserisci gli elementi della matrice" << endl;
	for (i = 1; i <= n; i++) {
		for (j = 1; j <= (n + 1); j++) {
			cout << "A[" << i << ", " << j << " ]=";
			cin >> a[i][j];
		}
	}
	
	for (j = 1; j <= n; j++) {
		for (i = 1; i <= n; i++) {
			if (i != j) {
				b = a[i][j] / a[j][j];
				for (k = 1; k <= n + 1; k++) {
					a[i][k] = a[i][k] - b * a[j][k];
				}
			}
		}
	}
	cout << "\nLa soluzione è:\n";
	for (i = 1; i <= n; i++) {
		x[i] = a[i][n + 1] / a[i][i];
		cout << "x" << i << "=" << x[i] << " ";
	}

Mi sono chiari il primo e l'ultimo pezzo, ma non quello centrale.

for (j = 1; j <= n; j++) {
		for (i = 1; i <= n; i++) {
			if (i != j) {
				b = a[i][j] / a[j][j];
				for (k = 1; k <= n + 1; k++) {
					a[i][k] = a[i][k] - b * a[j][k];
				}
			}
		}
	}

Qualcuno potrebbe per favre spiegarmi cosa effettua quel pezzo di codice?

[Super Squirrel]

Ora non ho tempo, ma ne approfitto per consigliarti il seguente link:

https://it.wikipedia.org/wiki/Metodo_di ... e_di_Gauss

[oleg.fresi]

Ti ringrazio, ma l'algoritmo in se l'ho capito, cioè saprei farlo manualmente su carta, ciò in cui trovo difficoltà è nella comprensine di quel pezzo di codice in cui viene formalizzato il procedimento.

[vict85]

Dove viene definito in quel codice il vettore \(\mathbf{b}\)? Insomma quello per cui si ha \(A\mathbf{x} =\mathbf{b}\)?

[oleg.fresi]

Penso che sia l'array

float x[10] 

[vict85]

Come non detto: nel codice inserisce una colonna in più nella matrice \(A\). Dall'ultimo ciclo è evidente che \(\mathbf{b}\) è nella \((n+1)\)-esima. È un po' strano che i cicli partano da 1 e alcune scelte computazionali sono strane (per esempio continua a sommare cose che ha teoricamente azzerato[nota]L'azzeramento non è sicuro per questioni di errori di calcolo, quindi quel codice aggiunge errori di approssimazione.[/nota]). Insomma ho l'impressione che tu stia guardando il codice di un principiante che ha copiato un codice scritto per matlab.

Ho riscritto una mia versione per eliminare alcuni calcoli inutili (anche se non ho aggiunti altri per questioni estetiche). Qui sotto trovi entrambe le versioni con una piccola modifica negli indici della tua (e l'aggiunta di un controllo per la divisione per 0). Ho scritto la soluzione direttamente nella colonna \(n+1\).

#include <iostream>

using namespace std;

constexpr size_t N = 3;

void
print( float a[ N ][ N + 1 ] )
{
    for ( size_t i = 0; i != N; i++ )
    {
        cout << "| ";
        for ( size_t j = 0; j != N + 1; j++ )
        {
            if ( j == N )
            {
                cout << "| " << ( i == 1 ? " = " : "   " ) << "| ";
            }
            cout << a[ i ][ j ] << " ";
        }
        cout << "|" << endl;
    }
}

bool solve_mine( float a[ N ][ N + 1 ] )  // sia b che x sono nella colonna N-esima
{
    for ( size_t j = 0; j != N; j++ )
    {
        if ( a[ j ][ j ] == 0. )
        {
            return false;  // se risolvibile, serve permutare le righe per risolverlo
        }

        for ( size_t i = 0; i != N; i++ )
        {
            if ( i == j )
                continue;  // non voglio modificare la riga j-esima

            float const b
                = a[ i ][ j ] / a[ j ][ j ];  // fattore che serve per azzerare l'elemento a[i][j]
            a[ i ][ j ] = 0.;                 // azzero piu' per questioni estetiche che per altro
            for ( size_t k = j + 1; k < N + 1; k++ )
            {
                // sottrae la riga j-esima moltiplicata per b alla riga i-esima
                // parto da j+1 perche' nelle colonne precedenti solo una delle due righe è non
                // nulla
                a[ i ][ k ] -= b * a[ j ][ k ];
            }
        }
    }

    // calcolo soluzione, ora la matrice è un sistema di 3 equazioni ad una sola incognita
    for ( size_t i = 0; i != N; i++ )
    {
        a[ i ][ N ] /= a[ i ][ i ];  // scrivo soluzione
        a[ i ][ i ] = 1.;            // questo e' per questioni estetiche
    }

    return true;
}

bool solve_yours( float a[ N ][ N + 1 ] )  // ho cambiato gli indici per conformita' con l'altro
{
    for ( size_t j = 0; j < N; j++ )
    {
        if ( a[ j ][ j ] == 0 )
            return false;  // aggiunto per evitare divisione per 0

        for ( size_t i = 0; i < N; i++ )
        {
            if ( i != j )
            {
                float const b = a[ i ][ j ] / a[ j ][ j ];
                for ( size_t k = 0; k < N + 1; k++ )
                {
                    a[ i ][ k ] -= b * a[ j ][ k ];
                }
            }
        }
    }
    for ( size_t i = 0; i < N; i++ )
    {
        a[ i ][ N ] /= a[ i ][ i ];
        a[ i ][ i ] = 1.0;  // aggiunto per questioni estetiche
    }
    return true;
}

int
main( )
{
    /*
     *   x + y + z = 6
     *  2x + y - z = 1
     *  2x -3y + z = -1
     *
     *  soluzione x = 1, y = 2, z = 3
     *
     *  L'ho trovata online, volevo un sistema con una soluzione per controllare che gli algoritmi
     *  fossero corretti
     */
    float a1[ N ][ N + 1 ] = {{1., 1., 1., 6.}, {2., 1., -1., 1.}, {2., -3., 1., -1.}};
    float a2[ N ][ N + 1 ] = {{1., 1., 1., 6.}, {2., 1., -1., 1.}, {2., -3., 1., -1.}};
    float reference_result[ N ] = {1., 2., 3.};

    print( a1 );
    cout << "solving mine..." << endl;

    if ( solve_mine( a1 ) )
    {
        cout << "La funzione ha prodotto un risultato. " << endl;
        float max_error = 0.0;
        for ( size_t i = 0; i != N; ++i )
        {
            max_error = max( max_error, abs( a1[ i ][ N ] - reference_result[ i ] ) );
        }
        cout << "La soluzione " << ( max_error < 1.E-4 ? "" : "non " )
             << "e' corretta, con un errore massimo di " << max_error << endl;
    }
    else
    {
        cout << "La funzione ha fallito" << endl;
    }

    cout << "La matrice e' ora cosi':" << endl;
    print( a1 );

    cout << endl << endl << "solving yours..." << endl;

    if ( solve_yours( a2 ) )
    {
        cout << "La funzione ha prodotto un risultato. " << endl;
        float max_error = 0.0;
        for ( size_t i = 0; i != N; ++i )
        {
            max_error = max( max_error, abs( a2[ i ][ N ] - reference_result[ i ] ) );
        }
        cout << "La soluzione " << ( max_error < 1.E-4 ? "" : "non " )
             << "e' corretta, con un errore massimo di " << max_error << endl;
    }
    else
    {
        cout << "La funzione ha fallito" << endl;
    }

    cout << "La matrice e' ora cosi':" << endl;
    print( a2 );
}

Sul mio pc le due soluzioni hanno lo stesso errore numerico, la ragione è che quel sistema è molto stabile numericamente e la maggior parte delle operazioni sono esatte. Pensa che avevo inizialmente fatto una ottimizzazione che consisteva nel precalcolare \(1/a[j][j]\) e aveva raddoppiato l'errore...

[oleg.fresi]

Ti ringrazio tanto, ora ho capito come funziona. Nella tua versione del codice, che sarebbe la funzione solvemine, però non capisco dove avvenga lo scambio tra righe se succede che la if restituisce false.

if ( a[ j ][ j ] == 0. ) 
{
        return false;  
} 

[vict85]

Lo scambio delle righe non c'è in quel codice. Comunque consiste in due fai: la ricerca della riga e lo scambio. Scambiare manualmente le righe non è efficiente per matrici di grandi dimensioni e potrebbe essere meglio usare puntatori alle righe e scambiare solo i puntatori.

Ma dovresti provare ad aggiungerlo da solo.

[oleg.fresi]

Vorrei capire meglio cosa fa il tuo codice. La funzione print stampa la matrice, poi ho fatto caso che la funzione solv_mine non è void, ma restituisce bool. Non capisco perchè questa scelta. Il risultato true o false è in base al fatto se il sistema ammette o no soluzione? Il primo ciclo for in pratica scandisce le colonne. Il primo if controlla se i numeri sulla diagonale sono diversi da zero. Ma se è zero, viene restituito false e cosa accade? Termina la funzione o va vanti? Il secondo for va a calclare il fattore che permette di azzerare. Il terzo ciclo for va a modificare la riga, azzerando l'elemento voluto. L'ultimo ciclo for se non sbaglio va a risolvere il sistema determinando tutte le incognite. Poi restituisce true, anche se non mi è chiaro perchè. Potresti gentilmente chiarirmi questi passaggi?

[vict85]

La funzione ritorna vero o falso a seconda che l'algoritmo sia capace di trovare una soluzione. Il sistema può avere soluzioni ma quel particolare algoritmo potrebbe non essere in grado di trovarlo. Ad esempio, il sistema \[\left(\begin{array}{ccc|c} 0 & 1 & 0 & 1\\ 1 & 0 & 0 & 2\\ 0 & 0 & 1 & 3\end{array}\right)\] ha ovviamente soluzioni ma l'algoritmo proposto non riesce a trovare la soluzione. Per trovarla devi prima scambiare la prima e la seconda riga (elemento che manca dal mio codice).
Esistono ovviamente sistemi che non possono essere risolti da questo tipo di algoritmi. Per esempio il sistema \[\left(\begin{array}{ccc|c} 1 & 2 & 1 & 5\\ 2 & 1 & 2 & 7\\ 1 & 1 & 1 & 4\end{array}\right)\]
Questo aspetto, però, dovrebbe essere stato trattato dal professore o dal libro. Insomma rientra nell'analisi dell'algoritmo stesso.

L'algoritmo non fa altro che cercare di produrre una matrice diagonale, insomma vuole che la colonna \(i\) abbia un singolo elemento non nullo nella riga \(i\). Per farlo sottrae la riga \(i\) a tutte le altre moltiplicandola per un fattore appropriato. Se \(a_{ii}= 0\) allora non puoi andare avanti.

[oleg.fresi]

per risolvere il sistema relativo alla prima matrice che hai scritto, basterebbe aggiungere una funzione che scambia le righe, ma nel secondo caso che algoritmo servirebbe. Ti premetto che l'algoritmo l'ho studiato in se e per se, se vogliamo l studio accorpato all'algebra lineare, ma tu mi pare stia alludendo all'analisi numerica che ancora non ho affrontato.

[Super Squirrel]

Per caratterizzare l'insieme delle soluzioni di un sistema di equazioni lineari basta ricorrere al teorema di Rouché-Capelli. Per farlo è necessario conoscere il rango della matrice completa ed incompleta, i quali possono essere facilmente calcolati trasformando la matrice completa in una matrice a scalini mediante l'utilizzo dell'algoritmo di Gauss. Fatte le opportune considerazioni, se si desidera calcolare le soluzioni del sistema non resta che applicare le mosse di Gauss alla matrice completa al fine di azzerare anche gli elementi (della matrice incompleta) alla destra della diagonale principale.

[oleg.fresi]

Perfetto, ora è chiaro, grazie mille!

[vict85]

"ZfreS":per risolvere il sistema relativo alla prima matrice che hai scritto, basterebbe aggiungere una funzione che scambia le righe, ma nel secondo caso che algoritmo servirebbe. Ti premetto che l'algoritmo l'ho studiato in se e per se, se vogliamo l studio accorpato all'algebra lineare, ma tu mi pare stia alludendo all'analisi numerica che ancora non ho affrontato.

Questo tipo di domande si fanno generalmente per gli esami di algebra lineare numerica. Per capire come implementare questo tipo di algoritmi non è sufficiente conoscere la teoria dell'algebra lineare. Detto questo, non ho fatto riferimento ad aspetti prettamente numerici ma a cose che si dovrebbero vedere anche in algebra lineare (tranne quando ho parlato di errori numerici).

Aspetti puramente numerici sono cose come la scelta del pivoting per ridurre gli errori numerici e il concetto di mal condizionamento di un sistema. Vi sono poi questioni puramente implementative come la scelta dell'ordine dei cicli o l'uso di versioni a blocchi per aumentare la cache locality del codice.