[MATEMARTI]

Ciao a tutti! Ho un problema con il criterio di arresto a posteriori nel codice Jacobi. Gira, ma scrive nel risultato 1.#QNANO ;-(
Qualche idea per aiutarmi? … La matrice prova è 5x5 ed ha 12 sulla diagonale principale e 1 altrove. Il vettore dei termini noti è (16,…,16). Grazie mille in anticipo! Questo è il codice :

/* QUESTO PROGRAMMA REALIZZA IL METODO DI JACOBI CON UN CRITERIO DI ARRESTO A POSTERIORI */
# include <stdio.h>
# include <stdlib.h>
typedef double matrice[51][51];
typedef double vettore[51];
void letvet(int,vettore);
void letmatr(int,matrice);
double modulo(double);
double calcola_lambda(int,matrice);
void controllo(double);
void vettore_d(int,vettore,vettore,matrice);
double delta_k(vettore,vettore,int);
void Jacobi(int,vettore,vettore,vettore,matrice);
void stampa(int,int,vettore);
int main()
{
    int n,cont=1; vettore b,vettk,vettkp1; matrice A; double lambda,tol,l,deltak;
    printf("\n Inserisci le dimensioni della matrice A(nxn). n="); scanf("%d",&n);
    letmatr(n,A);
    letvet(n,b);
    lambda=calcola_lambda(n,A); controllo(lambda);
    l=lambda/(1-lambda);
    printf("\n\n Tolleranza richiesta ="); scanf("%lf",&tol);
    vettore_d(n,b,vettkp1,A);
    // criterio di arresto a posteriori:
    do {Jacobi(n,b,vettk,vettkp1,A); cont++;
       deltak=delta_k(vettk,vettkp1,n);}
    while((deltak*l)>=tol);
    stampa(n,cont,vettkp1);
    system("pause");
    return 0;
}
void letvet(int n,vettore v)
{    printf("\n\n INSERISCI LE COMPONENTI DEL VETTORE b DEI TERMINI NOTI:\n");
     for(int i=1;i<=n;i++)
        {printf(" Componente %d =",i); scanf("%lf",&v[i]);}
     return;}
void letmatr(int n,matrice a)
{    printf("\n\n INSERISCI LE COMPONENTI DELLA MATRICE A DEI COEFFICIENTI:\n"); 
     for(int i=1;i<=n;i++)
        for(int j=1;j<=n;j++)
           {printf(" Componente %d,%d =",i,j); scanf("%lf",&a[i][j]);}
     return;}
double modulo(double x)
{    if(x<0) return -x;
     else return x;}
double calcola_lambda(int n,matrice a)
{    double lambda=0, somma;
     for(int i=1;i<=n;i++)
        {somma=0;
        for(int j=1;j<=n;j++)
           if(i!=j) somma+=a[i][j];
        somma/=modulo(a[i][i]);
        if(lambda<somma) lambda=somma;}
     return lambda;}
void controllo(double l)
{    if(l>=1)
       {printf("\n Il sistema non e' risolubile. "); system("pause"); exit(1);}
     return;}
void vettore_d(int n,vettore b,vettore vettkp1,matrice a)
{    for(int i=1;i<=n;i++)
        vettkp1[i]=b[i]/a[i][i]; return;}
void Jacobi(int n,vettore b,vettore x,vettore y,matrice a)
{    for(int i=1;i<=n;i++)
        {y[i]=b[i];
        for(int j=1;j<=n;j++)
           if(i!=j) y[i]-=a[i][j]*x[j];
        y[i]/=a[i][i];}
     return;}
void stampa(int n,int k,vettore x)
{    printf("\n\n ECCO IL VETTORE SOLUZIONE DEL SISTEMA:\n");
     printf(" (");
     for(int i=1;i<=n;i++)
        printf(" %lf",x[i]);
     printf(" )");
     printf("\n\n IL VETTORE SOLUZIONI E' STATO CALCOLATO ESEGUENDO %d ITERATE.\n\n\n ",k);
     return;}
double delta_k(vettore a,vettore b,int n)
{    double deltak=modulo(a[1]-b[1]);
     for(int i=2;i<=n;i++)
        if(modulo(a[i]-b[i])>deltak) deltak=modulo(a[i]-b[i]);
     return deltak;}