algorytm.org

Implementacja w C/C++



Baza Wiedzy
wersja offline serwisu przeznaczona na urządzenia z systemem Android
Darowizny
darowiznaWspomóż rozwój serwisu
Nagłówki RSS
Artykuły
Implementacje
Komentarze
Forum
Bookmarki






Sonda
Implementacji w jakim języku programowania poszukujesz?

Metoda Gaussa - Seidela - Implementacja w C/C++
Ocena użytkownikóww: *****  / 4
SłabyŚwietny
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;
}
Dodaj komentarz