[C++] NaN e inf... Come li evito?

[_BlackTornado_]

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;
     }

[Tommo]

Originariamente inviato da _BlackTornado_

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.

...

Niubbezza eh

premesso che non so niente del dominio del tuo problema, sicuramente il tuo è un algoritmo che richiede una precisione che i float non possono fornirti, da cui l'errore che si propaga e sballa tutti i risultati.
Infatti non credo proprio che Mathematica utilizzi float normali in calcoli del genere, userà un suo qualche tipo di dato molto ma molto più preciso.

Hai provato ad usare dei double invece dei float?

[_BlackTornado_]

Originariamente inviato da Tommo

Niubbezza eh

premesso che non so niente del dominio del tuo problema, sicuramente il tuo è un algoritmo che richiede una precisione che i float non possono fornirti, da cui l'errore che si propaga e sballa tutti i risultati.
Infatti non credo proprio che Mathematica utilizzi float normali in calcoli del genere, userà un suo qualche tipo di dato molto ma molto più preciso.

Hai provato ad usare dei double invece dei float?

Ho provato... L'unica cosa che ottengo è che ci mette di più a passare da "numeroni" a "inf" e poi "nan".

Il mio riferimento a Mathematica era unicamente per dire che non sto chiedendo al programma di dividere qualcosa per 0, altrimenti Mathematica me l'avrebbe detto. Tutto qui.

Non pretendo e non voglio ottenere la precisione di Mathematica, soltanto una serie che converga.

[Tommo]

Beh, se le iterazioni raddoppiano coi double, allora vuol dire che i numeri in gioco aumentano la loro risoluzione di un'ordine di grandezza ad ogni iterazione;
ad una certa il float fa overflow e diventa inf;
Quindi viene utilizzato come denominatore di una divisione inf/inf, e diventa NaN.
Evidentemente a lungo andare tutti i campi della matrice diventano prima inf e poi NaN.

Direi proprio che devi capire come modificare l'algoritmo per tenere i valori entro range accettabili!
Visto che i numeri sono così enormi, potresti provare a svolgere i calcoli in scala logaritmica?

[_BlackTornado_]

Originariamente inviato da Antonio23

se in un algoritmo iterativo come newton-raphson i numeri non fanno altro che crescere (l'algoritmo divergere) o è sbagliata la tua implementazione oppure è sbagliata la scelta dei parametri (degli stepsize ecc...)

La prima...

Facevo fare l'inverso dello Jacobiano... E poi nella procedura successiva ci mandavo lo Jacobiano invece che l'inverso :P

Adesso converge... A numeri apparentemente sbagliati, ma questa è un altra storia.

Grazie a tutti.