Inversione matrice in c++

[Keethra]

Salve a tutti!
E' da giorni che cerco disperatamente di implementare un codice che faccia l'inversione di una matrice. Io ho scritto alcune funzioni, messe in un file sorgente a parte, che vengono richiamate nel main. Una calcola la sottomatrice ottenuta eliminando colonna e riga dell'elemento di cui si vuole calcolare il complemento algebrico. Questa funzione viene richiamata per ricorsione dalla funzione determinante. La funzione inversa richiama a sua volta determinante e sottomatrice. Purtroppo questo codice non funziona. Sapete dirmi dove sta l'errore?

Questo è ciò che ho scritto io nel file sorgente dove contengono le funzioni che avete implementato:

double** sottomatrice(double** A, int i, int j, int dim){
// Allocazione della matrice
double** Aij = new double*[dim - 1];
for(int n = 0; n < dim - 1; n++)
{
Aij[n] = new double[dim - 1];
}
int h = 0;
int k;
for(int r = 0; r < dim; r++)
{
k=0;
//indice delle colonne della matrice Aij
if(r != i)
{
for(int c = 0; c < dim; c++)
if(c != j)
{
Aij[h][k] = A[r][c];
k++;
}
h++;
}
}
return Aij;
}


double determinante(double** A, int dim)
{
if(dim==1)
return A[0][0]; // Ritorna l'unico elemento della matrice
else
{
double ris=0.0;
for(int col=0, row=0; col<dim; col++)
{ // esegue la somma sugli elementi della prima riga.
// Chiama ricorsivamente la funzione determinante sulla matrice Aij
ris += pow(-1.0,(row+col)) * A[row][col] * determinante(sottomatrice(A, row, col, dim), dim-1);
}
return ris;
}
}


double** inversa(double** A, int dim)
{
// Allocazione della matrice inversa
double** invA;
invA = new double*[dim];
for(int i=0; i<dim; i++)
{
invA=new double[dim];
}
// Calcolo dell'inversa
if(dim == 1)
invA[0][0] = A[0][0];
else
{
double detA = determinante(A, dim);
for(int i=0; i < dim; i++)
{
for(int j=0; j < dim; j++)
{
invA[j]=(1/detA) * pow(-1.0,(i+j)) * determinante(sottomatrice(A,i,j,dim),(dim-1));
}
}
}
return invA;
}


questo è cio che ho scritto nel main:

double **a; //dichiarazione della matrice che contiene la copia della matrice da invertire

a = new double*[dim_ridotta];
if (a == NULL)
{
cout << "Non è stata allocata memoria per pt.\n";
exit(EXIT_FAILURE);
} ;
for (int i = 0; i < dim_ridotta; i++)
{
a= new double[dim_ridotta]; //copio kr e affianco la matrice identita'
if (a == NULL)
{
cout << "Non è stata allocata memoria per a.\n";
exit(EXIT_FAILURE);
}
}



for (int i = 0; i < dim_ridotta; i++)
{
for (int j = 0; j < dim_ridotta; j++)
{
a[j]=k_r[j];
}
}





double **inv; //dichiarazione la matrice che contiene l'inversa di k_rl
inv = new double*[dim_ridotta];
if (inv == NULL)
{
cout << "Non è stata allocata memoria per pt.\n";
exit(EXIT_FAILURE);
}
for (int i = 0; i < dim_ridotta; i++)
{
inv= new double[dim_ridotta]; //copio kr e affianco la matrice identita'
if (inv == NULL)
{
cout << "Non è stata allocata memoria per a.\n";
exit(EXIT_FAILURE);
}
}

double detA = determinante(a, dim_ridotta);
cout << "\nDeterminante: " << detA << "\n";

// Se la matrice non é singolare é possibile calcolare la sua inversa

if(detA != 0.0){
double** inva = inversa(a,dim_ridotta);
cout << "\nInversa:\n";
for (int a = 0; a < dim_ridotta; a++) {
for (int b = 0; b < dim_ridotta; b++)
cout << "[" << inva[a] <<"]";
cout << "\n";
}
for (int i = 0; i < dim_ridotta; i++)
{
for (int i = 0; i < dim_ridotta; i++)
{
inv[j]=inva[j];
}
}
}
else
cout << "\n Matrice singolare, impossibile invertirla!";

inserendo i seguenti dati questi sono i risultati:

2.1213 0.7071 -0.7071
0.7071 2.1213 -0.7071
-0.7071 -0.7071 0.7071

ottengo un determanante pari a 1.4142 (giusto)
e questa matrice inversa:

1.7247*e^-350 0 0
1.7247*e^-350 0 0
1.7247*e^-350 0 0

che ovviamente non è giusto.

Aiutatemi, son disperata!!!!!

[Ziosilvio]

Prendi un libro di Algebra Lineare, studia l'algoritmo di inversione di Gauss, e implementa quello.

[Keethra]

Con Gauss ci avevo già provato ma senza grossi risultati...

[Ziosilvio]

Si vede che non ci hai provato abbastanza

Scherzi a parte, implementare l'inversione di matrici usando l'algoritmo di Cramer è un suicidio in termini di tempo di esecuzione, e probabilmente anche di complessità della struttura dati necessaria per rappresentare "bene" le matrici.
Forse a questo punto un libro di Analisi Numerica può essere di aiuto più di uno di Algebra Lineare.
Puoi provare con le Numerical Recipes, che esistono in versione C, C++, Fortran 77 e Fortran 90.
Leeggi bene la licenza d'uso, che è piuttosto limitativa: in particolare, non puoi includere parti del libro sotto forma di codice sorgente.

[Keethra]

Facendo le giuste correzioni ora funziona quello che vi ho scritto qui

[Keethra]

Ecco con l'altro metodo, ma funziona fino ad un certo punto. Perchè?

int main()
{

int dim_ridotta;
cout << "Inserisci la dimensione della matrice: ";
cin >> dim_ridotta;


double s;
double inversa[dim_ridotta][dim_ridotta];
double k_r[dim_ridotta][dim_ridotta];
double a[dim_ridotta][2*dim_ridotta];


cout << "\n\nInserisci la matrice da invertire:\n";
for (int i = 0; i < dim_ridotta; i++) //scorro le righe
{ //inizio for i
for (int j = 0; j < dim_ridotta; j++)
{ //inizio for j
cout <<"["<< i + 1 <<"]"<<"["<< j + 1<<"]"<< ": ";
cin >> k_r[i][j]; //copio kr nella parte sinistra della matrice
} //fine for j
} //fine for i
cout <<endl;

//copio kr
for (int i = 0; i < dim_ridotta; i++) //scorro le righe
{ //inizio for i
for (int j = 0; j < dim_ridotta; j++)
{ //inizio for j
a[i][j] = k_r[i][j]; //copio kr nella parte sinistra della matrice
} //fine for j
} //fine for i


//affianco l'identita'
for (int l = 0; l < dim_ridotta; l++) //scorro le righe
{ //inizio for l
for (int k = dim_ridotta; k < 2*(dim_ridotta); k++)
{ //inizio for k
if(l == k - dim_ridotta)
a[l][k] = 1; //assegno valore 1 lungo la diagonale
else a[l][k] = 0; //assegno valore 0 a tutti gli altri elementi
} //fine for k
} //fine for l


//inversione della matrice di rigidezza ridotta
//inizio dell'applicazione del teorema di Jordan
for (int i = 0; i < dim_ridotta; i++)
{ //inizio for i
if(a[i][i]==0) //inizio scambio_righe
{
double u[2*dim_ridotta];
for (int j = i+1; j < dim_ridotta; j++)
{ //inizio for m
if (a[j][i] != 0)
{ //inizio if
for (int k = i; k < 2*dim_ridotta; k++)
{ //inizio for k
u[k] = a[i][j];
a[i][k] = a[j][k];
a[j][k] = u[k];
}
j = dim_ridotta+1;
} //fine if
if (j != dim_ridotta+1)
{ //inizio if
printf(" MATRICE NON INVERTIBILE!! ");
printf(" Premere un tasto per continuare ");
scanf("%c");
} //fine if
} //fine for j
} //fine scambio_righe
s=a[i][i];
for (int j = i; j < 2*dim_ridotta; j++)
{//inizio for j
a[i][j]=a[i][j]/s;
} //fine for j
for (int k = i+1; k < dim_ridotta; k++)
{ //inizio for k
s=a[k][i];
for (int j = i; j < 2*dim_ridotta; j++)
{//inizio for j
a[k][j] = a[k][j]-a[i][j]*s;
} //fine for j
} //fine for k
}//fine for i

if (a[(dim_ridotta - 1)][(dim_ridotta - 1)] == 0)
{ //inizio if
printf(" MATRICE NON INVERTIBILE!! ");
printf(" Premere un tasto per continuare ");
scanf("%c");

} //fine if
for (int k = dim_ridotta; k < 2*dim_ridotta; k++)
{//inizio for k
a[(dim_ridotta-1)][k] = a[(dim_ridotta-1)][k]/a[(dim_ridotta-1)][(dim_ridotta-1)] ;
a[(dim_ridotta-1)][(dim_ridotta-1)] = 1;
}//fine for k

for (int i = 1; i < dim_ridotta; i++)
{ //inizio for i
for (int k = 0; k < i; k++)
{ //inizio for k
s=a[k][i];
for (int j = i; j < (2*dim_ridotta); j++)
{ //inizio for j
a[k][j] = a[k][j]-a[i][j]*s;
} //fine for j
} //fine for k
} //fine for i
//fine applicazione della jordan_form



//creazione della matrice inversa
for (int i = 0; i < dim_ridotta; i++)
{ //inizio for i
for(int j = dim_ridotta; j < 2*dim_ridotta; j++)
{
inversa[i][(j - dim_ridotta)]= a[i][j];
}
}//fine for i




printf( "\n\nLa matrice inversa e' \n\n");
for (int i = 0; i < dim_ridotta; i++)
{ //inizio for i
for(int j = 0; j < dim_ridotta; j++)
{
cout << inversa[i][j] << "\t ";
} //fine for j
cout <<"\n";
}//fine for i



system("PAUSE");
return 0;
}