Nadesłany przez Codename, 08 kwietnia 2011 03:26
Kod przedstawiony poniżej przedstawia główną część rozwiązania problemu.Pobierz pełne rozwiązanie.
Jeżeli nie odpowiada Ci sposób formatowania kodu przez autora skorzystaj z pretty printer'a i dostosuj go automatycznie do siebie.
gs_2_c.c:
/* Metoda Gaussa-Seidela */ /* Jacek Czarniecki */ /* www.algorytm.org */ #include <stdio.h> #define MAX 100 const int MAX_ITER = 200; // maksymalna liczba iteracji const float EPS = 0.000001; // przybliżenie const float INIT_X = 0.0; // wartości wektora poczatkowego float L[MAX][MAX],D[MAX],U[MAX][MAX],b[MAX],x[MAX],tmp[MAX]; int n; float I[MAX][MAX],TL[MAX][MAX]; float min(float x, float y){ return x<y?x:y; } float f_abs(float x){ return x>0?x:-x; } int read(){ int i, j; scanf("%i",&n); // wczytanie rozmiaru macierzy for(i=0;i<n;i++){ for(j=0;j<n;j++){ // wczytywanie bezpośrednio do L,D,U if (i < j) {scanf("%f",&U[i][j]);} else if (i > j) {scanf("%f",&L[i][j]);} else {scanf("%f",&D[i]); } } } //na przekątnej nie powinno być zer for(i=0;i<n;i++) scanf("%f",&b[i]); return 0; } int copy(float *a, float *b,int n){ while(n--) b[n] = a[n]; return 0; } float cmp(float *a, float *b,int n, float eps){ float w,v=0; int flag = 0; // sprawdzenie warunku zakonczenia: ||x(k)-x(k-1)|| <= epsilon while(n--) { if(f_abs(a[n]-b[n]) > eps) {flag = 1; break;} } return flag; } int test_zero(float *a, int n){ while(n--) if(a[n]==0) return 0; return 1; // czy na przekatnej wystepuja zera? } float norma1(){ float n1,h; int i,j; for (n1=i=0; i < n; i++) { for (h=j=0; j < n; j++) h += f_abs(TL[j][i]); if (h > n1) n1 = h; } // norma "jeden" return n1; } float norma2(){ float n2,h; int i,j; for (n2=i=0; i < n; i++) { for (h=j=0; j < n; j++) h += f_abs(TL[i][j]); if (h > n2) n2 = h; } // norma "nieskonczonosc" return n2; } int zbiega(){ int i,j,k; float h; // wyznaczenie (L+D)^-1 * N for (i=0;i<n;i++) for (j=0;j<i;TL[i][j]=L[i][j++]); // TL <- L - kopia macierzy L for (i=0;i<n;TL[i][i]=D[i++]); // TL <- L + D for (i=0;i<n;i++) for (j=0;j<n;j++) I[i][j]=(i!=j?0:1); for (i=0;i<n;i++){ // TL <- (L + D)^-1 for (j=0;j<=i;j++) I[i][j]/=D[i]; for (k=i+1;k<n;k++) for (j=0;j<=k;j++) I[k][j]-=(I[i][j]*TL[k][i]); } for (i = 0; i < n; i++) { // TL <- (L + D)^-1 * N for (j = 0; j < n; j++) { TL[i][j] = 0; for (k = 0; k <= i; k++) TL[i][j] += (I[i][k] * U[k][j]); } } return (min(norma1(), norma2()) <= 1)?1:0; } int solve(){ int i,j,k,counter; read(); // wczytywanie macierzy if(test_zero(D, n)) { // sprawdzanie przekątnej if(zbiega()){ // warunek zbieżności for (i=0; i<n; i++) { D[i] = 1/D[i]; // D <- D^-1 b[i] *= D[i]; // b <- D^-1 * b for (j=0; j<i; j++) L[i][j] *= D[i]; // L <- D^-1 * L for (j=i+1; j<n; j++) U[i][j] *= D[i]; // U <- D^-1 * U x[i] = INIT_X; // x <- [0,...,0] } counter = 0; do{ // kolejne iteracje copy(x,tmp,n); // tmp[] <- x[] for (i=0; i<n; i++) x[i] = b[i]; for (i=0; i<n; i++) { // x = D^-1*b - for (j=0; j<i; j++) x[i] -= L[i][j]*x[j]; // D^-1*L * x - for (j=i+1; j<n; j++) x[i] -= U[i][j]*tmp[j]; // D^-1*U * x } counter++; }while( cmp(x,tmp,n,EPS) && counter < MAX_ITER); // warunki zakonczenia iteracji for(i=0;i<n;i++) printf("%.0f ",x[i]); putchar('\n'); // wypisanie wyniku } else printf("Nie spelniono warunku zbieznosci\n"); } else printf("Na przekatnej wystepuja zera\n"); return 0; } int main(){ int t; scanf("%i",&t); while(t--) solve(); return 0; }