Obliczanie wyznacznika macierzy

movax1
Użytkownik
Użytkownik
Posty: 64
Rejestracja: 3 paź 2009, o 12:15
Płeć: Mężczyzna
Podziękował: 16 razy
Pomógł: 2 razy

Obliczanie wyznacznika macierzy

Post autor: movax1 »

Witam,
Jako projekt piszę program obliczający wyznacznik macierzy dowolnego rozmiaru.
Do stopnia 10 oblicza on wyznacznik w czasie poniżej sekundy
Stopień 11: 5 sekund
Stopień 12: 20 sekund
Stopień 13: 12 minut
Stopień 14: 2 godziny

Oczywiście dotyczy to pesymistycznego przypadku, gdy w macierzy nie ma zer ;)
Wyznacznik jest liczony przy użyciu reguły Sarrusa i rozwinięcia Laplace'a.

I tak się zastanawiam: Jak można to zoptymalizować? Czy są może jakieś własności wyznacznika, z których można skorzystać?

Rozważam też jakieś operacje elementarne, by wyzerować co nieco - jednak pytanie: czy opłaca się dorzucać funkcję sprowadzająca jakieś elementy do zera? Selekcjonowanie, które wiersze należy dodać/odjąć również pochłania nieco obliczeń.


Pozdrawiam
Awatar użytkownika
kadiii
Użytkownik
Użytkownik
Posty: 642
Rejestracja: 20 gru 2005, o 21:04
Płeć: Mężczyzna
Lokalizacja: Wrocław
Pomógł: 130 razy

Obliczanie wyznacznika macierzy

Post autor: kadiii »

... ementation masz przykładowe rozwiązania twojego problemu
spajder
Użytkownik
Użytkownik
Posty: 735
Rejestracja: 7 lis 2005, o 23:56
Płeć: Mężczyzna
Lokalizacja: Łódź
Podziękował: 2 razy
Pomógł: 133 razy

Obliczanie wyznacznika macierzy

Post autor: spajder »

Zdecydowanie operacje elementarne. Najpierw tak ułóż wiersze, aby na głównej przekątnej nie było zer, potem wyzeruj wszystko poniżej głównej przekątnej. Jak się nie da to wyznacznik wynosi 0. W innym wypadku wyznacznik to to iloczyn elementów głównej przekątnej.
movax1
Użytkownik
Użytkownik
Posty: 64
Rejestracja: 3 paź 2009, o 12:15
Płeć: Mężczyzna
Podziękował: 16 razy
Pomógł: 2 razy

Obliczanie wyznacznika macierzy

Post autor: movax1 »

Dzięki za linka, zaraz poczytam

Czyli nie muszę sprowadzać elementów do postaci z "wiodącymi jedynkami"? Po prostu - wystarczy aby na głównej przekątnej nie było zer, a potem zerować pozostałe elementy?
spajder
Użytkownik
Użytkownik
Posty: 735
Rejestracja: 7 lis 2005, o 23:56
Płeć: Mężczyzna
Lokalizacja: Łódź
Podziękował: 2 razy
Pomógł: 133 razy

Obliczanie wyznacznika macierzy

Post autor: spajder »

Metod jest kilka, to co napisałem jest jedną z nich i w zupełności wystarczy.
Awatar użytkownika
Mariusz M
Użytkownik
Użytkownik
Posty: 6903
Rejestracja: 25 wrz 2007, o 01:03
Płeć: Mężczyzna
Lokalizacja: 53°02'N 18°35'E
Podziękował: 2 razy
Pomógł: 1246 razy

Obliczanie wyznacznika macierzy

Post autor: Mariusz M »

Metoda permutacyjna (rekurencyjna)

Kod: Zaznacz cały

#include<stdio.h>
#include<stdlib.h>
#include<conio.h>

/*
Obliczanie wyznacznika metodą generowania permutacji
Rekurencyjnie generowane permutacje zapisywane są do pliku
Ustalany jest znak (parzystość) permutacji
następnie sumowane są wszystkie iloczyny elementów
których pierwsze indeksy tworzą permutację początkową 1...n
a drugie indeksy są permutacją liczb 1...n
Każdą permutację sumuje się tylko raz
*/

void perm(int,int,int*,int*);
double det(int,double**);

int main(){
char ch;
int i,j,n;
double** a;
do{
printf("Podaj n=");
scanf("%d",&n);

a=(double**)malloc((n+1)*sizeof(double));
for(i=1;i<=n;i++)
a[i]=(double*)malloc((n+1)*sizeof(double));

for(i=1;i<=n;i++){
for(j=1;j<=n;j++){
printf("A[%d][%d]=",i,j);
scanf("%lf",&a[i][j]);
}
printf("\n");
}
printf("det(A)=%lf\n",det(n,a));
for(i=1;i<=n;i++)
free(a[i]);
free(a);
ch=getch();
}
while (ch!=27);
return 0;
}


void perm(int poz,int n,int* t,int* s){

int i,j;
int wystapil;

FILE* fp;

if((fp=fopen("perm.txt","a"))==NULL){
printf("Nie mozna otworzyc pliku do dopisywania \n");
exit(1);
}

if(poz>n){
for(i=1;i<=n;i++)
fprintf(fp,"%d ",t[i]);
fprintf(fp,"\n");
}
else
for(j=1;j<=n;j++){
wystapil=0;
for (i=1;i<=poz-1;i++)
if(t[i]==s[j]) wystapil=1;
if(!wystapil){
t[poz]=s[j];
perm(poz+1,n,t,s);
}
}
fclose(fp);
}

double det(int n,double** a){

int *r,*s,*t;
int i,j;
FILE* fp;
double d,x;

remove("perm.txt");

s=(int*)malloc((n+1)*sizeof(int));
t=(int*)malloc((n+1)*sizeof(int));

for(i=1;i<=n;i++)
s[i]=i;
printf("\n");

perm(1,n,t,s);
free(s);
free(t);

if((fp=fopen("perm.txt","r"))==NULL){
printf("Nie mozna otworzyc pliku do dopisywania \n");
exit(1);
}
r=(int*)malloc((n+1)*sizeof(int));
d=0;
while(!feof(fp)){
for(i=1;i<=n;i++){
fscanf(fp,"%d",&r[i]);
}
fscanf(fp,"\n");
r[0]=1;
for(i=1;i<=n;i++)
for(j=i+1;j<=n;j++)
if(r[i]>r[j]) r[0]*=-1;
x=r[0];
for(i=1;i<=n;i++)
x*=a[i][r[i]];
d+=x;
}
fclose(fp);
free(r);
remove("perm.txt");
return d;
}
Metoda rozwinięć Laplace

Kod: Zaznacz cały

#include<stdio.h>
#include<stdlib.h>
#include<conio.h>

/*
Obliczanie wyznacznika metodą rozwinięć Laplace
*/

double det(double**, int);

int main(){
int i,j,n;
double** a;
char ch;
do{
printf("Podaj n=");
scanf("%d",&n);
a=(double**)malloc((n+1)*sizeof(double));
for(i=0;i<n+1;i++)
a[i]=(double*)malloc((n+1)*sizeof(double));
for(i=0;i<n;i++){
    for(j=0;j<n;j++){
      printf("A[%d][%d]=",i+1,j+1);
      scanf("%lf",&a[i][j]);
    }
printf("\n");
}
printf("det(A)=%lf\n",det(a,n));
for(i=0;i<n+1;i++)
free(a[i]);
free(a);
ch=getch();
}
while(ch!=27);
return 0;
}

double det(double** a, int size)
{
	int i, k, l, z = 1;
	double** b;
	double d = 0;
	b=(double**)malloc((size+1)*sizeof(double));
	for(i=0;i<size+1;i++)
	b[i]=(double*)malloc((size+1)*sizeof(double));
	if (size == 1)
		d = a[0][0];
	else
	{
		// rozwiniecie Laplace'a względem pierwszego wiersza
		for (i = 0; i < size; i++)
		{
			// przepisywanie wierszy
			for (k = 0; k < size - 1; k++)
			{
				// przepisywanie kolumn
		    	for (l = 0; l < size - 1; l++)
				{
						// przepisanie komórki znajdującej się w tej
						// samej kolumnie lecz w wierszu poniżej
					if (l < i)  b[k][l] = a[k + 1][l];
						// przepisanie komórki znajdującej się
						// w następnej kolumnie i wierszu poniżej
					else  b[k][l]=a[k + 1][l + 1];
				}
			}
			// wywołanie rekurencyjne
			d = d + z * a[0][i] * det(b, size-1);
			// zmiana znaku
			z = -z;
		}
	}
	for(i=0;i<size+1;i++)
	free(b[i]);
	free(b);
	return d;
}
Metoda eliminacji Gaussa

Kod: Zaznacz cały

#include<stdio.h>
#include<stdlib.h>
#include<conio.h>
#include<math.h>

/*
Obliczanie wyznacznika metodą eliminacji Gaussa
*/

double det(int n,double** a);

int main(){
char ch;
int i,j,n;
double** a;
do{
printf("Podaj n=");
scanf("%d",&n);

a=(double**)malloc((n+1)*sizeof(double));
for(i=1;i<=n;i++)
a[i]=(double*)malloc((n+1)*sizeof(double)); //allokacja pamięci na tablicę

for(i=1;i<=n;i++){
for(j=1;j<=n;j++){
printf("A[%d][%d]=",i,j);
scanf("%lf",&a[i][j]);
}
printf("\n");
}
printf("det(A)=%lf\n",det(n,a));
for(i=1;i<=n;i++)
free(a[i]);
free(a);  //zwalnianie pamięci
ch=getch();
}
while (ch!=27);
return 0;
}

double det(int n,double** a){
int i,j,k,wm;
double p,w=1.0,max;
for (i=1;i<n;i++){
  max=fabs(a[i][i]);
	wm=i;                          
for(j=i+1;j<=n;j++)
  if(fabs(a[j][i])>max){
	max=fabs(a[j][i]);      
	wm=j;                        //wyszukujemy element główny w kolumnie
  }
 if(wm>i) {
  for(k=i;k<=n;k++){
	p=a[i][k];                //gdy znajdziemy element główny to zamieniamy wiersze               
  a[i][k]=a[wm][k];
a[wm][k]=p;
  }
 w*=-1;           //zamiana wierszy zmienia znak wyznacznika 
 }
 if (a[i][i]==0) w=0;     // kolumna jest zerowa wyznacznik jest równy zero
 if (w!=0)
 for(j=i+1;j<=n;j++){
	p=(double)a[j][i]/a[i][i];
  for(k=n;k>i;k--)
  a[j][k]-=p*a[i][k];              //właściwy etap eliminacji (sprowadzanie do postaci trójkątnej)
 }
 w*=a[i][i];
 }
 w*=a[n][n];           // wyznacznik macierzy trójkątnej to iloczyn elementów na głównej przekątnej
 return w;
}
Metoda rozkładu LU

Kod: Zaznacz cały

#include<stdio.h>
#include<stdlib.h>
#include<conio.h>
#include<math.h>

#define TINY 1.0e-15

/*
Obliczanie wyznacznika metodą rozkładu LU=PA
*/

void ludcmp(double** a,int n,int* indx,double* d);

int main(){
char ch;
int i,j,n;
int* indx;
double det;
double **a;
do{
printf("Podaj n=");
scanf("%d",&n);
a=(double**)malloc((n+1)*sizeof(double));
for(i=0;i<n+1;i++)
a[i]=(double*)malloc((n+1)*sizeof(double));
indx=(int*)malloc((n+1)*sizeof(int));
for(i=1;i<=n;i++){
for(j=1;j<=n;j++){
printf("A[%d][%d]=",i,j);
scanf("%lf",&a[i][j]);
}
printf("\n");
}
ludcmp(a,n,indx,&det);
for(i=1;i<=n;i++) det*=a[i][i];
printf("det A=%lf\n",det);
for(i=0;i<n+1;i++)
free(a[i]);
free(a);
free(indx);
ch=getch();
}
while(ch!=27);
return 0;
}

void ludcmp(double** a,int n,int* indx,double* d){
int i,imax,j,k;
double big,dum,sum,temp;
double *vv;
vv=(double*)malloc((n+1)*sizeof(double));
(*d)=1.0;
for(i=1;i<=n;i++){
  big=0.0;
  for(j=1;j<=n;j++)
  if((temp=fabs(a[i][j]))>big) big=temp;
  if(big==0.0){
  (*d)=0.0;
  return;
  }
  vv[i]=1.0/big;
  }
  for(j=1;j<=n;j++){
   for(i=1;i<j;i++){
   sum=a[i][j];
   for(k=1;k<i;k++) sum-=a[i][k]*a[k][j];
   a[i][j]=sum;
   }
   big=0.0;
   for(i=j;i<=n;i++){
   sum=a[i][j];
   for(k=1;k<j;k++)
    sum-=a[i][k]*a[k][j];
	a[i][j]=sum;
	if((dum=vv[i]*fabs(sum))>=big){
	 big=dum;
	 imax=i;
	 }
	 }
	 if(j!=imax){
	 for(k=1;k<=n;k++){
	 dum=a[imax][k];
	 a[imax][k]=a[j][k];
	 a[j][k]=dum;
	 }
	 (*d)=-(*d);
	 vv[imax]=vv[j];
	 }
	 indx[j]=imax;
	 if(a[j][j]==0) a[j][j]=TINY;
	 if(j!=n){
	  dum=1.0/(a[j][j]);
	  for(i=j+1;i<=n;i++) a[i][j]*=dum;
	 }
}
free(vv);
}
movax1
Użytkownik
Użytkownik
Posty: 64
Rejestracja: 3 paź 2009, o 12:15
Płeć: Mężczyzna
Podziękował: 16 razy
Pomógł: 2 razy

Obliczanie wyznacznika macierzy

Post autor: movax1 »

Dziękuję raz jeszcze za kod; Czy mógłbym jeszcze prosić o ogólne wyjaśnienie, na jakiej zasadzie działa poniższa pętla? Poszczególne operacje oczywiście rozumiem, ale nie bardzo widzę ogólny algorytm obliczania tego.

Kod: Zaznacz cały

 
 if (w!=0)
 for(j=i+1;j<=n;j++)
 {
    p=(double)a[j][i]/a[i][i];
    for(k=n;k>i;k--)
    a[j][k]-=p*a[i][k];              //właściwy etap eliminacji (sprowadzanie do postaci trójkątnej)
  }
ODPOWIEDZ