Algorytm Gaussa dla liczb zmiennoprzecinkowych

0

Jednym z ważnych algorytmów mocno polegających na odróżnieniu liczby zerowej od niezerowej jest algorytm eliminacji Gaussa tj sprowadzania macierzy do postaci schodkowej całkowicie zredukowanej. Aby właściwie zaimplementować ten algorytm w arytmetyce zmiennoprzecinkowej należy go odpowiednio zmodyfikować
1 Przed obliczeniami należy macierz znormalizować tj obliczyć największy co do wartości bezwzględnej element macierzy  i podzielić przez niego wszystkie elementy macierzy. (otrzymujemy w ten sposób macierz o wartościach z przedziału od -1 do 1).
2 W zwykłym algorytmie Gaussa używamy dowolnego elementu niezerowego kolumny do uzyskania wartości 1 i wyzerowania reszty kolumny. Aby uzyskać algorytm numerycznie stabilny, należy zastąpić pytanie o zerowanie się elementów kolumny porównaniem z epsilonem maszynowym (mała w porównaniu z 1.0 liczbą dodatnia. Zamiast brać pierwszy lepszy element niezerowy powinniśmy wybrać wiersz zawierający w danej kolumnie element największy co do wartości bezwzględnej (oczywiście pod warunkiem, że kolumna zawiera choćby jeden element większy od epsilon - w przeciwnym wypadku kolumnę należy uznać za zerową).
3 Przykład. Załóżmy, że macierz (a właściwie  pierwsza kolumna macierzy) ma postać
4                   0.125           
5                  -0.5               
6                   0.25             
7 W przypadku zwykłego algorytmu użylibyśmy pierwszego wiersza do otrzymania 1.0 na pozycji (1,1) (w notacji algebraicznej) i wyzerowania pozostałych elementów pierwszej kolumny. W przypadku omawianej wersji należałoby najpierw przestawić wiersze pierwszy i drugi (gdyż -0.5 jest największą co do wartości bezwzględnej liczbą w pierwszej kolumnie). Analogicznie należy postępować w następnych kolumnach.
8 Przy drukowaniu postaci zredukowanej zastąp elementy macierzy wynikowej mniejsze od epsilona zerami.
W załączonym pliku zdefiniowano pewną macierz. Dopisz do tego pliku kod odpowiedzialny za eliminację Gaussa w/g przedstawionej powyżej wersji algorytmu. Pamiętaj, że kod powinien być uniwersalny, niezależny od konkretnej macierzy. Przetestuj kod na przykładach różnych macierzy (różnych funkcji zapełniających macierz liczbami.

Napisałem trochę kodu, konkretnie od komentarza /tu powinien zacząc się algorytm eliminacji/ przed tym komentarzem program jest dany do zadania. Kompletnie nie wiem co dalej zrobić, proszę o jakieś podpowiedzi.

#include<stdio.h>
#include<math.h>
/* dyrektywa preprocesora: M i N zostaną przed kompilacją zastąpione odpowiednimi napisami*/
#define M 50
#define N 100

const double eps = 1e-12; /* stała przybliżenia zera*/

int main(void){

float a[M][N], maksimum, bufor;
int m,n,i,j,w;

/*użytkownik podaje prawdziwe rozmiary macierzy*/
printf("podal liczbe wierszy macierzy <= %d\n", M);
scanf("%d", &m);
printf("podal liczbe wierszy macierzy <= %d\n", N);
scanf("%d", &n);

/*nadanie początkowych wartosci macierzy*/
for( i=0; i<m; i++)
    for( j=0; j<n; j++ )
        a[i][j] = i-j;

     maksimum = fabsf(a[0][0]);
        w=0;
    /*tu powinien zacząc się algorytm eliminacji*/

    /*szukania wartości bewzględnej największego elementu*/

    for( i=0; i<m; i++)
    for( j=0; j<n; j++ )
       if(fabsf(a[i][j])>fabsf(maksimum))
           maksimum = a[i][j];

    /*normalizacja macierzy*/
    for( i=0; i<m; i++)
       for( j=0; j<n; j++)
       {
           a[i][j] = a[i][j]/maksimum;
       }

    /*szukanie elementu największego w kolumnie*/
    for( j=0; j<n; j++)
    {
        for( i=0; i<m; i++)
        {
            if(fabsf(a[i][j])>fabsf(maksimum))
            {
                maksimum = a[i][j];
                w=i;
            }
        }

        /*zamiana  wierszy*/

        for(i=0; i<n; i++)
        {

            bufor = a[0][i];

            a[0][i] = a[w][i];

            a[w][i] = bufor;

        }

    }

/*wydruk zredukowanej macierzy*/

return 0;
}
3

Nie pisz wszystkiego w main, będzie łatwiej.
Kiedyś napisałem zwykłą eliminację Gaussa w C https://wandbox.org/permlink/ktsqg27KDrddnreY możesz zacząć od tego.

1

Albo mi się wydaje, albo tworzysz macierz prostokątną... Można też liczyć, ale tu raczej będą rozwiązania metodami przybliżonymi i może być więcej niż jedno rozwiązanie wtedy...

0
#include<stdio.h>
#include<math.h>
/* dyrektywa preprocesora: M i N zostaną przed kompilacją zastąpione odpowiednimi napisami*/
#define M 50
#define N 100

const double eps = 1e-12; /* stała przybliżenia zera*/

void normalizacja(int m, int n, float a[M][N])
{
    int i,j;
    float maksimum = a[0][0];

    /*szukania wartości bewzględnej największego elementu*/

    for( i=0; i<m; i++)
    for( j=0; j<n; j++ )
       if(fabsf(a[i][j])>fabsf(maksimum))
           maksimum = a[i][j];

    /*normalizacja macierzy*/
    for( i=0; i<m; i++)
       for( j=0; j<n; j++)
       {
           a[i][j] = a[i][j]/maksimum;
       }

}

void swap_row(int w, int n, int q, float a[M][N]) /* zamiana wierszy*/
{
    int i;
    float bufor;

    for(i=0; i<n; i++)
       {

           bufor = a[q][i];

           a[q][i] = a[w][i];

           a[w][i] = bufor;

       }
}

int main(void){

    float a[M][N], maksimum,t;
int m,n,i,j,w,k,b;

/*użytkownik podaje prawdziwe rozmiary macierzy*/
printf("podal liczbe wierszy macierzy <= %d\n", M);
scanf("%d", &m);
printf("podal liczbe wierszy kolumn <= %d\n", N);
scanf("%d", &n);

/*nadanie początkowych wartosci macierzy*/
for( i=0; i<m; i++)
    for( j=0; j<n; j++ )
        a[i][j] = i-j;

        w=0;

    /* drukowanie macierzy przed dzialaniami algorytmu- do testow */
    for( i=0; i<m; i++)
            {
            for( j=0; j<n; j++ )
            {
                printf("%.2f \t",a[i][j]);
            }
                printf("\n");
            }

     printf("\n");

    /*tu powinien zacząc się algorytm eliminacji*/

    normalizacja(m, n, a);

        for(int q=0; q<n; q++) /* postać schodkowa*/
        {
            maksimum = a[q][q];
            for( b=q; b<m; b++) /* szukam maksimum w kolumnie*/
                 {
                     if(fabsf(a[b][q])>fabsf(maksimum))
                     {
                         maksimum = a[b][q];
                         w=b;
                     }
                 }

                 swap_row(w, n, q, a); /*zamiana  wierszy*/

            for(i=q+1; i<m; i++) /* eliminacja współczynników*/

            {
                if(fabs(a[q][q])>eps)
                {
                t = -a[i][q]/a[q][q];

                    for(k=0;k<n; k++)
                    {
                        a[i][k] = t* a[q][k] + a[i][k];
                    }

                }

            }

        }

    for(i=0; i<m; i++) /* postać z 1*/
    {

        for(j=0; j<n; j++)
        {

            if(fabs(a[i][j])>eps)
            {
                t = 1/a[i][j];

                for(k=j; k<n; k++)
                {
                    a[i][k]= a[i][k]*t;
                }

                j=n-1;
            }

        }

    }

/*wydruk zredukowanej macierzy*/

        for( i=0; i<m; i++)
         {
         for( j=0; j<n; j++ )
         {
             printf("%.2f \t",a[i][j]);
         }
             printf("\n");
         }

return 0;
}

na razie mam taki kod i nie wiem jak wyzerować nad jedynkami wiodącymi o ile w ogóle dobrze ten kod liczy. Ktoś ma jakieś sugestię jak usprawnić mój kod, żeby działał poprawnie?

2

1 Nie wprost na temat: strasznie długi main() -> bo za mało funkcji. Problem ujawniają komentarze (wiele byłoby zbędne przy dobrej nazwie funkcji). Zdecydowanie ŹLE to się analizuje

  1. Formalności: jak już masz const epsylon (bardzo dobra praktyka) to M, N też.

  2. epsylon jest 1e-12, a oparłeś się na typie float, który ma zaledwie 7-8 cyfr znaczących. W tych czasach nie ma żadnych powodów, by nie używać double

0

Powie ktoś mi dlaczego ta funkcja nie działa?

Zapomniałem dodać co ona ma robić.

Ta funkcja ma usuwać, elementy nad jedynka wiodącą

void usuwanie(int m, int n, float a[M][N])
{

    int k,i,j;
    float t;
    for(i = m-1; i>0; i--)
    {

        for(j=0; j<n; j++)
        {
            if(fabsf(a[i][j])>eps)
            {
                for(k = i-1; k>=0; k--)
                {
                    t = -a[k][j]/a[i][j];

                    a[k][j] = t*a[i][j] + a[k][j];
                }
            }

            j=n-1;
        }

    }
1
Daim123 napisał(a):

Powie ktoś mi dlaczego ta funkcja nie działa?

Zapomniałem dodać co ona ma robić.

Ta funkcja ma usuwać, elementy nad jedynka wiodącą

Nie usuwa, bo niczego takiego nie zaprogramowałeś (kurtuazyjnie zakładając, ze jesteś tego autorem, w co wątpię)

0

Akurat, że sam to napisałem. : )

0

Jako dowód faza projektu funkcji :)

0

Kurczę niech ktoś mi pomoże, dla nadania wartości a[i][j]=i-j działa ok, ale już dla a[i][j]=i-j/2 nie . Ja już nie mam siły na to zadanie....

#include<stdio.h>
#include<math.h>
/* dyrektywa preprocesora: M i N zostanπ przed kompilacjπ zastπpione odpowiednimi napisami*/
#define M 50
#define N 100

const double eps = 0.000001; /* stała przybliżenia zera*/

int check_column(int m, int k, float a[M][N]) /* sprawdza czy kolumna nie jest zerowa*/
{
    int i;

    for(i=0; i<m; i++)
    {
        if(fabs(a[i][k])>eps)
        {
            return 1;

        }
    }
        return 0;

}

int check_row(int n, int k, float a[M][N]) /* sprawdza czy kolumna nie jest zerowa*/
{
    int i;

    for(i=0; i<n; i++)
        if(fabs(a[k][i])>eps) return 1;

    return 0;

}

void normalizacja(int m, int n, float a[M][N])
{
    int i,j;
    float maksimum = a[0][0];

    /*szukania wartości bewzględnej największego elementu*/

    for( i=0; i<m; i++)
    for( j=0; j<n; j++ )
       if(fabs(a[i][j])>fabs(maksimum))
           maksimum = a[i][j];

    /*normalizacja macierzy*/
    for( i=0; i<m; i++)
       for( j=0; j<n; j++)
       {
           a[i][j] = a[i][j]/maksimum;
       }

}

void swap_row(int w, int n, int q, float a[M][N]) /* zamiana wierszy*/
{
    int i;
    float bufor;

    for(i=0; i<n; i++)
       {

           bufor = a[q][i];

           a[q][i] = a[w][i];

           a[w][i] = bufor;

       }
}

void eliminacja(int i, int j, int m, int n, float a[M][N])
{
    int  k, e;
    float t;

        for(k = i+1; k< m; k++)
        {
            t = -a[k][j]/a[i][j];

            for(e = j; e < n; e++)
            {
                a[k][e] += a[i][e]*t;
            }
        }

}

void sort_col(int m, int n, float a[M][N]) /* sortowanie wierszy*/
   {
       int q,b;
       float maksimum;
       for( q=0; q<n; q++) /* postać schodkowa*/
                {
                    maksimum = a[q][q];
                    for( b=q; b<m; b++) /* szukam maksimum w kolumnie*/
                         {
                             if(fabs(a[b][q])>fabs(maksimum))
                             {
                                 maksimum = a[b][q];
                                swap_row(b, n, q, a);
                             }
                         }

                }

   }
void sort_col2(int j, int m, int n, float a[M][N]) /* sortowanie wierszy*/
{
    int q,b;
    float maksimum;
    for( q=j; q<n; q++) /* postać schodkowa*/
             {
                 maksimum = a[q][q];
                 for( b=q; b<m; b++) /* szukam maksimum w kolumnie*/
                      {
                          if(fabs(a[b][q])>fabs(maksimum))
                          {
                              maksimum = a[b][q];
                             swap_row(b, n, q, a);
                          }
                      }

             }

}

void gauss(int m, int n, float a[M][N])
{

    int i, j;

    for(j = 0; j<n; j++)
    {

       if(check_column(m, j, a)==1)
       {
           sort_col2(j, m, n, a);
           for(i=j; i<m; i++)
           {
               if(fabs(a[i][j])>0)
               {
                   eliminacja(i, j, m, n, a);
               }

           }
       }
    }

}

void jedyneczka(int m, int n, float a[M][N])
{
    int i, j, k;
    float t;
    for(i=0; i<m; i++) /* postać z 1*/
      {

          for(j=0; j<n; j++)
          {

              if(fabs(a[i][j])>eps)
              {
                  t = 1/a[i][j];

                  for(k=j; k<n; k++)
                  {
                      a[i][k]= a[i][k]*t;
                  }

                  j=n-1;
              }

          }

      }

}

void zero(int m, int n, float a[M][N])
{
    int i, j;

    for( i=0; i<m; i++)
    for( j=0; j<n; j++ )
        if(fabs(a[i][j])<eps)
            a[i][j]=0;

}

int check_column2(int m, int k, float a[M][N]) /* sprawdza czy kolumna nie jest zerowa*/
{
    int i;

    for(i=m-1; i>=0; i--)
    {
        if(fabs(a[i][k])>eps)
        {
            return 1;

        }
    }
        return 0;

}
void eliminacja2(int i, int j, int m, int n, float a[M][N])
{
    int  k, e;
    float t;

        for(k = i-1; k>= 0; k--)
        {
            t = -a[k][j]/a[i][j];

            for(e = j; e < n; e++)
            {
                a[k][e] += a[i][e]*t;
            }
        }

}

void gauss2(int m, int n, float a[M][N]) /* zeruje elementy powyżej 1*/
{

    int i, j;
    for(i = m - 1; i>=0; i--)
    {

       if(check_row(n, i, a)==1)
       {
           for(j=0; j<n; j++)
           {
               if(fabs(a[i][j])>eps)
               {
                   eliminacja2(i, j, m, n, a);
                   j=n-1;
               }

           }

       }
    }

}

int main(void){

    float a[M][N];
    int m,n,i,j;

/*uøytkownik podaje prawdziwe rozmiary macierzy*/
printf("podal liczbe wierszy macierzy <= %d\n", M);
scanf("%d", &m);
printf("podal liczbe wierszy macierzy <= %d\n", N);
scanf("%d", &n);

/*nadanie poczπtkowych wartosci macierzy*/
for( i=0; i<m; i++)
    for( j=0; j<n; j++ )
        a[i][j] = i-j;

/*tu powinien zaczπc się algorytm eliminacji*/
    normalizacja(m, n, a);
    gauss(m, n, a);
    gauss2(m, n, a);
    jedyneczka(m, n, a);
    zero(m, n, a);

    /*wydruk zredukowanej macierzy*/

   for( i=0; i<m; i++)
        {
        for( j=0; j<n; j++ )
        {
            printf("%.2f \t",a[i][j]);
        }
            printf("\n");
        }

return 0;
}

1 użytkowników online, w tym zalogowanych: 0, gości: 1, botów: 0