_BlackTornado_
06-11-2010, 16:02
Salve a tutti...
Purtroppo, nella mia niubbagine, torno da voi, perchè mi trovo con un problema che non so risolvere.
Sto facendo un programma che risolve numericamente con Newton-Raphson un particolare sistema di equazioni non lineari.
Si tratta quindi di prendere matrici, fare jacobiani, invertirle e poi moltiplicarle e sottrarle per qualche vettore. Il tutto, iterato per tot volte.
In pratica, già dopo pochi cicli, alcuni valori mi schizzano verso più infinito, altri verso meno infinito.
Dopo 5 cicli mi ritrovo come soluzione inf inf -inf -inf.
Se lo mando avanti per 100 cicli, mi trovo nan nan nan nan.
Ora, nella mia completa niubbaggine, non mi ero mai trovato ad affrontare casi del genere (quando il programma mi ha dato "nan" ho dovuto cercare su google il significato), ed è un comportamento che non sto in alcun modo ricercando (il programma dovrebbe restituire dei numeri).
Mi sono accorto che la funzione che "inverte lo Jacobiano" mi restituisce, in alcuni valori, un "-0" invece che 0. Può essere questo che mi sputtana tutto? Come faccio ad evitarlo?
Un altra possibile causa potrebbe essere che alcuni determinanti che "dovrebbero" venire 0 vengono in realtà un numero molto piccolo, ma non 0.
L'invertiMatrice quindi va a dividere per questi "quasi 0" e ottiene numeri enormi.
Ho provato ad inserire un controllo sul determinante, dicendo di non invertire se il determinante è minore di 0.0001, ma apparentemente non è cambiato niente.
PS: Mathematica risolve lo stesso sistema con Newton-Raphson senza incontrare determinanti uguali a 0, o almeno credo.
PS: Ho dimenticato di dire che sono ragionevolmente sicuro che le varie funzioni "funzionino" in linea di massima, perchè ho fatto numerose prove.
Quando però mi trovo con tanti "0" nella matrice, quando vado a invertire, in questo caso diventano "-0", e il programma "impazzisce".
Questi sono i sottoprogrammi in questione:
bool risoluzioneNewtonRaphson (float pEnne[4], float gdl, float tolleranza1, float tolleranza2){
int i, j;
float funzione[4], Jacobiano[4][4], JacInverso[4][4], prodotto[4], diffSoluzioni[4], pEnnepUno[4];
for (i=0; i<5; i++){
creaFunzione(pEnne, gdl, funzione);
//cout << "funzione" << endl;
//stampaVettore(funzione);
creaJacobiano(pEnne, Jacobiano);
cout << endl << endl << "jacobiano";
stampaMatrice(Jacobiano);
if (invertiMatrice(Jacobiano, JacInverso)==false){
return false;
}
cout << endl << endl <<"JacInverso";
stampaMatrice(JacInverso);
calcolaProdotto(Jacobiano, funzione, prodotto);
//cout << endl <<"Prodotto" << endl;
//stampaVettore(prodotto);
sottraiVettori(pEnne, prodotto, pEnnepUno);
//cout << endl << "pEnnepUno" << endl;
//stampaVettore(pEnnepUno);
sottraiVettori(pEnnepUno, pEnne, diffSoluzioni);
creaFunzione(pEnnepUno, gdl, funzione);
for (j=0; j<4; j++){
pEnne[j] = pEnnepUno[j];
}
if (moduloVettore(diffSoluzioni) <= tolleranza2 && moduloVettore(funzione) <= tolleranza1){
return true;
}
}
return false;
}
void creaFunzione(float pEnne[], float gdl, float funzione[4]){
funzione[0] = l2*cos(gdl) + pEnne[0]*cos(pEnne[1]);
funzione[1] = l1 + l2*sin(gdl) + pEnne[0]*sin(pEnne[1]);
funzione[2] = l2*cos(gdl) + (l3l4 - pEnne[0])*cos(pEnne[1]) + pEnne[2];
funzione[3] = l2*sin(gdl) + (l3l4 - pEnne[0])*sin(pEnne[1]) + pEnne[3] + l6;
return;
}
void creaJacobiano(float pEnne[], float Jacobiano[4][4]){
Jacobiano[0][0] = cos(pEnne[1]);
Jacobiano[1][0] = sin(pEnne[1]);
Jacobiano[2][0] = -1*cos(pEnne[1]);
Jacobiano[3][0] = -1*sin(pEnne[1]);
Jacobiano[0][1] = -1*pEnne[0]*sin(pEnne[1]);
Jacobiano[1][1] = pEnne[0]*cos(pEnne[1]);
Jacobiano[2][1] = -1*(l3l4 - pEnne[0])*sin(pEnne[1]);
Jacobiano[3][1] = (l3l4 - pEnne[0])*cos(pEnne[1]);
Jacobiano[0][2] = 0;
Jacobiano[1][2] = 0;
Jacobiano[2][2] = 1;
Jacobiano[3][2] = 0;
Jacobiano[0][3] = 0;
Jacobiano[1][3] = 0;
Jacobiano[2][3] = 0;
Jacobiano[3][3] = 1;
return;
}
bool invertiMatrice(float matr[4][4], float inv[4][4]){
int i, j;
float matrid[4][4], complAlg[4][4];
// In questo IF, il determinante non e' zero, perchè con le approssimazioni float e' praticamente impossibile che
// det venga esattamente 0, ma in genere viene un numero molto piccolo.
if (det(4, matr)<0.0001){
cout << "il determinante e':" << det(4, matr) << endl;
cout << endl << "La matrice non e' invertibile!" << endl;
return false;
}
complAlgebrici(matr, matrid, complAlg);
calcolaTrasposta(complAlg, inv);
for(i=0; i<4; i++){
for(j=0; j<4; j++){
inv[i][j]=inv[i][j]/(det(4, matr));
}
}
return true;
}
void complAlgebrici(float matr[4][4], float matrid[4][4], float complAlg[4][4]){
int i, j;
for (i=0; i<4; i++){
for (j=0; j<4; j++){
riduciMatrice(matr, matrid, i, j);
if (pd(i-j)==true){
complAlg[i][j] = det(4-1, matrid);
}
else {
complAlg[i][j] = -1*det(4-1, matrid); // Attenzione: qui è meglio mettere -1* invece di -det per evitare roba che diventi -0.
}
}
}
return;
}
float det(int dim, float matr[4][4]){
int i,j,k;
float determinante = 0;
float matrid[4][4];
if (dim==1){
determinante = matr[0][0];
}
for (i=0; i<dim && dim>1; i++){
for (j=0; j<dim; j++){
for (k=1; k<dim && j!=i; k++){
if (j<i){matrid[j][k-1]=matr[j][k];}
if (j>i){matrid[j-1][k-1]=matr[j][k];}
}
}
determinante=determinante + pot(i) * matr[i][0] * det(dim-1,matrid);
}
return determinante;
}
int pot(int i){
if(i%2==0){
return 1;
}
else {
return -1;
}
}
void calcolaTrasposta(float matr[4][4], float trasp[4][4]){
int i, j;
for (i=0; i<4; i++){
for (j=0; j<4; j++){
trasp[j][i] = matr[i][j];
}
}
return;
}
void riduciMatrice(float matr[4][4], float matrid[4][4], int i, int j){
int r, c, contarighe, contacolonne;
contarighe = 0;
contacolonne = 0;
for (r=0; r<4; r++){
for(c=0; c<4; c++){
if (r!=i && c!=j){
matrid[contarighe][contacolonne] = matr[r][c];
contacolonne++;
if (contacolonne == 4-1){
contarighe ++;
contacolonne=0;
}
}
}
}
return;
}
Purtroppo, nella mia niubbagine, torno da voi, perchè mi trovo con un problema che non so risolvere.
Sto facendo un programma che risolve numericamente con Newton-Raphson un particolare sistema di equazioni non lineari.
Si tratta quindi di prendere matrici, fare jacobiani, invertirle e poi moltiplicarle e sottrarle per qualche vettore. Il tutto, iterato per tot volte.
In pratica, già dopo pochi cicli, alcuni valori mi schizzano verso più infinito, altri verso meno infinito.
Dopo 5 cicli mi ritrovo come soluzione inf inf -inf -inf.
Se lo mando avanti per 100 cicli, mi trovo nan nan nan nan.
Ora, nella mia completa niubbaggine, non mi ero mai trovato ad affrontare casi del genere (quando il programma mi ha dato "nan" ho dovuto cercare su google il significato), ed è un comportamento che non sto in alcun modo ricercando (il programma dovrebbe restituire dei numeri).
Mi sono accorto che la funzione che "inverte lo Jacobiano" mi restituisce, in alcuni valori, un "-0" invece che 0. Può essere questo che mi sputtana tutto? Come faccio ad evitarlo?
Un altra possibile causa potrebbe essere che alcuni determinanti che "dovrebbero" venire 0 vengono in realtà un numero molto piccolo, ma non 0.
L'invertiMatrice quindi va a dividere per questi "quasi 0" e ottiene numeri enormi.
Ho provato ad inserire un controllo sul determinante, dicendo di non invertire se il determinante è minore di 0.0001, ma apparentemente non è cambiato niente.
PS: Mathematica risolve lo stesso sistema con Newton-Raphson senza incontrare determinanti uguali a 0, o almeno credo.
PS: Ho dimenticato di dire che sono ragionevolmente sicuro che le varie funzioni "funzionino" in linea di massima, perchè ho fatto numerose prove.
Quando però mi trovo con tanti "0" nella matrice, quando vado a invertire, in questo caso diventano "-0", e il programma "impazzisce".
Questi sono i sottoprogrammi in questione:
bool risoluzioneNewtonRaphson (float pEnne[4], float gdl, float tolleranza1, float tolleranza2){
int i, j;
float funzione[4], Jacobiano[4][4], JacInverso[4][4], prodotto[4], diffSoluzioni[4], pEnnepUno[4];
for (i=0; i<5; i++){
creaFunzione(pEnne, gdl, funzione);
//cout << "funzione" << endl;
//stampaVettore(funzione);
creaJacobiano(pEnne, Jacobiano);
cout << endl << endl << "jacobiano";
stampaMatrice(Jacobiano);
if (invertiMatrice(Jacobiano, JacInverso)==false){
return false;
}
cout << endl << endl <<"JacInverso";
stampaMatrice(JacInverso);
calcolaProdotto(Jacobiano, funzione, prodotto);
//cout << endl <<"Prodotto" << endl;
//stampaVettore(prodotto);
sottraiVettori(pEnne, prodotto, pEnnepUno);
//cout << endl << "pEnnepUno" << endl;
//stampaVettore(pEnnepUno);
sottraiVettori(pEnnepUno, pEnne, diffSoluzioni);
creaFunzione(pEnnepUno, gdl, funzione);
for (j=0; j<4; j++){
pEnne[j] = pEnnepUno[j];
}
if (moduloVettore(diffSoluzioni) <= tolleranza2 && moduloVettore(funzione) <= tolleranza1){
return true;
}
}
return false;
}
void creaFunzione(float pEnne[], float gdl, float funzione[4]){
funzione[0] = l2*cos(gdl) + pEnne[0]*cos(pEnne[1]);
funzione[1] = l1 + l2*sin(gdl) + pEnne[0]*sin(pEnne[1]);
funzione[2] = l2*cos(gdl) + (l3l4 - pEnne[0])*cos(pEnne[1]) + pEnne[2];
funzione[3] = l2*sin(gdl) + (l3l4 - pEnne[0])*sin(pEnne[1]) + pEnne[3] + l6;
return;
}
void creaJacobiano(float pEnne[], float Jacobiano[4][4]){
Jacobiano[0][0] = cos(pEnne[1]);
Jacobiano[1][0] = sin(pEnne[1]);
Jacobiano[2][0] = -1*cos(pEnne[1]);
Jacobiano[3][0] = -1*sin(pEnne[1]);
Jacobiano[0][1] = -1*pEnne[0]*sin(pEnne[1]);
Jacobiano[1][1] = pEnne[0]*cos(pEnne[1]);
Jacobiano[2][1] = -1*(l3l4 - pEnne[0])*sin(pEnne[1]);
Jacobiano[3][1] = (l3l4 - pEnne[0])*cos(pEnne[1]);
Jacobiano[0][2] = 0;
Jacobiano[1][2] = 0;
Jacobiano[2][2] = 1;
Jacobiano[3][2] = 0;
Jacobiano[0][3] = 0;
Jacobiano[1][3] = 0;
Jacobiano[2][3] = 0;
Jacobiano[3][3] = 1;
return;
}
bool invertiMatrice(float matr[4][4], float inv[4][4]){
int i, j;
float matrid[4][4], complAlg[4][4];
// In questo IF, il determinante non e' zero, perchè con le approssimazioni float e' praticamente impossibile che
// det venga esattamente 0, ma in genere viene un numero molto piccolo.
if (det(4, matr)<0.0001){
cout << "il determinante e':" << det(4, matr) << endl;
cout << endl << "La matrice non e' invertibile!" << endl;
return false;
}
complAlgebrici(matr, matrid, complAlg);
calcolaTrasposta(complAlg, inv);
for(i=0; i<4; i++){
for(j=0; j<4; j++){
inv[i][j]=inv[i][j]/(det(4, matr));
}
}
return true;
}
void complAlgebrici(float matr[4][4], float matrid[4][4], float complAlg[4][4]){
int i, j;
for (i=0; i<4; i++){
for (j=0; j<4; j++){
riduciMatrice(matr, matrid, i, j);
if (pd(i-j)==true){
complAlg[i][j] = det(4-1, matrid);
}
else {
complAlg[i][j] = -1*det(4-1, matrid); // Attenzione: qui è meglio mettere -1* invece di -det per evitare roba che diventi -0.
}
}
}
return;
}
float det(int dim, float matr[4][4]){
int i,j,k;
float determinante = 0;
float matrid[4][4];
if (dim==1){
determinante = matr[0][0];
}
for (i=0; i<dim && dim>1; i++){
for (j=0; j<dim; j++){
for (k=1; k<dim && j!=i; k++){
if (j<i){matrid[j][k-1]=matr[j][k];}
if (j>i){matrid[j-1][k-1]=matr[j][k];}
}
}
determinante=determinante + pot(i) * matr[i][0] * det(dim-1,matrid);
}
return determinante;
}
int pot(int i){
if(i%2==0){
return 1;
}
else {
return -1;
}
}
void calcolaTrasposta(float matr[4][4], float trasp[4][4]){
int i, j;
for (i=0; i<4; i++){
for (j=0; j<4; j++){
trasp[j][i] = matr[i][j];
}
}
return;
}
void riduciMatrice(float matr[4][4], float matrid[4][4], int i, int j){
int r, c, contarighe, contacolonne;
contarighe = 0;
contacolonne = 0;
for (r=0; r<4; r++){
for(c=0; c<4; c++){
if (r!=i && c!=j){
matrid[contarighe][contacolonne] = matr[r][c];
contacolonne++;
if (contacolonne == 4-1){
contarighe ++;
contacolonne=0;
}
}
}
}
return;
}