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
1). 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\).
- Codice:
#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...
