25/05/12 15:30
Puffetta
Ciao a tuttti!
ho scritto questo programma che deve risolvere un sistema lineare del tipo Ax=b con il metodo di gauss o di cholesky(metodo scelto dall'utente) ho però un grosso problema: il metodo di gauss funziona perfettamente(a quanto pare), mentre se uso il metodo di cholesky il programma va in loop. mi potreste aiutare a trovare l'errore? io non riesco a vederlo:-(
grazie a tutti!!!!
ho scritto questo programma che deve risolvere un sistema lineare del tipo Ax=b con il metodo di gauss o di cholesky(metodo scelto dall'utente) ho però un grosso problema: il metodo di gauss funziona perfettamente(a quanto pare), mentre se uso il metodo di cholesky il programma va in loop. mi potreste aiutare a trovare l'errore? io non riesco a vederlo:-(
grazie a tutti!!!!
/*Risoluzione di sistemi lineari usando metodi diretti quali metodo di eliminazione di gauss, cholesky e thomas*/
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<time.h>
typedef double matrice[300][300];
typedef double vettore[300];
int sceglimatrice();
void creamatrice(int mat, matrice A, vettore b, int n);
double U();
void prodotto(matrice A, matrice B, int n);
int sceglimetodo(int mat);
void gauss(matrice A, vettore b, vettore x, int n);
void riduci(matrice, vettore, int);
int cercapivot(vettore, int, int);
void scambiarighe(matrice, int, int, int);
void cholesky(matrice A, vettore b, vettore x, int n);
void risolvi_sup(matrice A, vettore b,vettore x, int n);
void risolvi_inf(matrice A, vettore b,vettore x, int n);
int main()
{
int n, metodo, mat;
matrice A;
vettore x, b;
printf("Risolviamo il sistema lineare Ax=b.\nInseriamo tutti i dati riguardanti la matrice A\n");
mat=sceglimatrice();
printf("\nQuante righe e colonne ha la matrice A? n=");
scanf("%d", &n);
creamatrice(mat, A, b, n);
metodo=sceglimetodo(mat);
/*
for(int i=1; i<=n; i++)
{
for(int j=1; j<=n; j++)
printf(" %lf ", A[i][j]);
printf("\n");
}
for(int i=1; i<=n; i++) printf("\t%lf", b[i]);
printf("\nmetodo=%d\n", metodo);
*/
switch(metodo)
{
case 1: gauss(A, b, x, n);
break;
case 2: cholesky(A, b, x, n);
break;
}
printf("\n\n Le soluzioni del sistema sono:");
for(int i=1; i<=n; i++) printf("\t%lf", x[i]);
system("PAUSE");
return 0;
}
double U()
{
return (double)rand()/RAND_MAX;
}
int sceglimatrice()
{
int mat=0;
printf("Scegliere il tipo da matrice da utilizzare:\n");
printf("\nMATRICE RANDOM\n\n");
printf(" 1) Matrice randomica generica\n\n");
printf(" 2) Matrice randomica simmetrica e definita positiva\n\n");
printf(" 3) Matrice randomica tridiagonale\n\n");
printf(" 4) Matrice randomica diagonale dominante\n");
printf("\nMATRICE INSERITA DA TASTIERA\n\n");
printf(" 5) Matrice generica\n\n");
printf(" 6) Matrice simmetrica e definita positiva\n\n");
printf(" 7) Matrice tridiagonale\n\n");
printf(" 8) Matrice diagonale dominante\n\n");
printf("Quale tipo di matrice devo usare? ");
scanf("%d", &mat);
while((mat<=0)||(mat>8))
{
printf("\nIl valore inserito non e' corretto. Inserire nuovamente il valore ");
scanf("%d", &mat);
}
return mat;
}
//Questa funzione crea la matrice in base alle esigenze dell'utente
void creamatrice(int mat, matrice A, vettore b, int n)
{
srand(time(0));
switch(mat)
{
//Matrice randomica generica e vettore randomico
case 1:
{
for(int i=1; i<=n; i++)
{
for(int j=1; j<=n; j++)
A[i][j]=U();//matrice
b[i]=U();//vettore
}
}
break;
//Matrice randomica simmetrica e definita positiva e vettore randomico
case 2:
{
matrice B;
for(int i=1; i<=n; i++)
{
for(int j=i; j<=n; j++)
{
B[i][j]=U();
B[j][i]=B[i][j];
}
b[i]=U();
}
//Uso la matrice B per creare la matrice A richiesta facendo A=(Bt)*B con B simmetrica
prodotto(A, B, n);
}
break;
//Matrice randomica tridiagonale e vettore randomico
case 3:
{
for(int i=1; i<=n; i++)
{
A[i][i]=U();
if(i!=n)
A[i+1][i]=U();
if(i!=1)
A[i-1][i]=U();
b[i]=U();
}
}
break;
case 4:
{
//diagonale dominante!!!!!!
}
break;
}
//Matrici inserite da tastiera. uso dunque lo scanf per acquisire le componenti
if((mat==5)||(mat==6)||(mat==7)||(mat==8))
{
printf("\nInserire i coefficienti della matrice A:");
for(int i=1; i<=n; i++)
for(int j=1; j<=n; j++)
scanf("%lf", &(A[i][j]));
printf("\nInserire il vettore dei termini noti");
for(int i=1; i<=n; i++)
{
printf("Inserire il termine noto b[%d]=", i);
scanf("%lf", &b[i]);
}
}
}
int sceglimetodo(int mat)
{
int metodo=0;
printf("\nQuale metodo devo usare per risolvere il sistema lineare quadrato Ax=b?\n1.Gauss 2.Cholesky 3.Thomas\n");
scanf("%d", &metodo);
while((metodo<=0)||(metodo>4))
{
printf("\n Il valore inserito non e' corretto!!! Scegliere di nuovo il metodo ");
scanf("%d", &metodo);
}
if((metodo==2)&&((mat==2)||(mat==6)))
return metodo;
if((metodo==3)&&((mat==3)||(mat==7)))
return metodo;
while((metodo!=1)&&(mat!=1)&&(mat!=5))
{
printf("\nIl metodo scelto non puo' essere usato sulla matrice A.\nscegliere un altro metodo ");
scanf("%d", &metodo);
if((metodo==2)&&((mat==2)||(mat==6)))
break;
if((metodo==3)&&((mat==3)||(mat==7)))
break;
}
return metodo;
}
void gauss(matrice A, vettore b, vettore x, int n)
{
riduci(A, b, n);
risolvi_sup(A, b, x, n);
}
void riduci(matrice A, vettore b, int n)
{
double mik; vettore kcolonna; int kmax;
for(int k=1; k<=n-1; k++)
{
for(int h=k; h<=n; h++)
{
kcolonna[h]=A[h][k];
kmax=cercapivot(kcolonna, k, n);
if(kmax!=k) scambiarighe( A, k, kmax, n);
}
for(int i=k+1; i<=n; i++)
{
mik=A[i][k]/A[k][k];
for(int j=k+1; j<=n; j++) A[i][j]=A[i][j]-mik*A[k][j];
b[i]=b[i]-mik*b[k];
}
}
return;
}
int cercapivot(vettore kcolonna, int k, int n)
{
double maxmodulo=fabs(kcolonna[k]);
int imax=k;
for(int i=k+1; i<=n; i++)
if(!maxmodulo) exit (1);
return imax;
}
void scambiarighe(matrice A, int k, int kmax, int n)
{
double aux;
for(int j=k; j<=n; j++)
{
aux=A[j][k];
A[j][k]=A[j][kmax];
A[j][kmax]=aux;
}
return;
}
void cholesky(matrice A, vettore b, vettore x, int n)
{
printf(" ********** ");
matrice H, HT;
vettore y;
printf(" ++++++++++ ");
//creo la matrice H
for(int i=1; i<=n; i++)
for(int j=1; j<=n; j++)
{
if(i==j==1)
H[i][j]=sqrt(A[i][i]);
if(j<i)
H[i][j]=0;
if(j>i)
{
double somma=0;
for(int k=1; k<=j-1; k++)
somma+=H[i][k]*H[j][k];
H[i][j]=(A[i][j]-somma)/H[j][j];
}
if(j==i)
{
double somma=0;
for(int k=1; k<=i-1; k++)
somma+=pow(H[i][k], 2);
H[i][j]=A[i][i]-somma;
}
}
for(int i=1; i<=n; i++)
for(int j=1; j<=n; j++)
HT[i][j]=H[j][i];
risolvi_inf(HT, b, y, n);
risolvi_sup(H, y, x, n);
}
//Questa funzione mi risolve un sistema la cui matrice e' triangolare inferiore
void risolvi_inf(matrice A, vettore b,vettore x, int n)
{
x[1]=b[1]/A[1][1];
for(int i=2; i<=n; i++)
{
for(int j=1; j<=i; j++) b[i]=b[i]-A[i][j]*x[j];
x[i]=b[i]/A[i][i];
}
return;
}
//Questa funzione mi risolve un sistema la cui matrice e' triangolare superiore
void risolvi_sup(matrice A, vettore b,vettore x, int n)
{
x[n]=b[n]/A[n][n];
for(int i=n-1; i>=1; i--)
{
for(int j=i+1; j<=n; j++) b[i]=b[i]-A[i][j]*x[j];
x[i]=b[i]/A[i][i];
}
return;
}
//Questa matrice mi fa il prodotto (Bt)*B e pone il risultato della nuova matrice A
void prodotto(matrice A, matrice B, int n)
{
for(int i=1; i<=n; i++)
for(int j=1; j<=n; j++)
{
A[i][j]=0;
for(int k=1; k<=n; k++)
A[i][j]=A[i][j]+B[i][k]*B[j][k];
}
return;
}
aaa
