Obliczanie wyznacznika macierzy
-
- 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
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
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
-
- 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
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.
-
- 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
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?
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?
- Mariusz M
- 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
Metoda permutacyjna (rekurencyjna)
Metoda rozwinięć Laplace
Metoda eliminacji Gaussa
Metoda rozkładu LU
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;
}
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;
}
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;
}
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);
}
-
- 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
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)
}