Elektroda.pl
Elektroda.pl
X

Search our partners

Find the latest content on electronic components. Datasheets.com
Please add exception to AdBlock for elektroda.pl.
If you watch the ads, you support portal and users.

Interpolcaja trygonometryczna - algorytm Goertzla

torcek 14 Feb 2009 15:43 2104 10
  • #1
    torcek
    Level 10  
    Mam problem z programem!
    Wszystko ładnie kompiluje ale wyświetla mi złe wartości więc pewnie jest coś źle jest w programie. Męczę się z tym od tygodnia i nie potrafię tego zrobić. Jeśli macie jakieś propozycje to napiszcie.
    Używam kompilatora Dev c++

    O to kod programu:

    Code:
    #include <stdio.h>
    
    #include <string>
    #include <stdlib.h>
    #include <math.h>

    #define MAXN 101
    // MAXN - rozmiar tablic węzłów i wartości funkcji w węzłach interpolacji
    //        zakładamy, że n<MAXN
       int n=6; // maksymalny indeks węzła interpolacji, (n+1) - lb. węzłów
       int m=20; // współczynnik m we wzorze (24)
       int mn; // lb. pośrednich punktów między węzłami interpolacji w drugim pliku   
       int ddelta; // współczynnik delta we wzorze (24)
       double a, b; // granice przedziału zawierającego węzły interpolacji
       double pi2=8*atan(1.0);  // funkcja na podstawie której wyznaczane sš warto?ci w węzłach interpolacji
       double x[MAXN]; // tablica węzłów na osi x
       double y[MAXN]; // tablica węzłów na osi y
       double f[MAXN]; // tablica wartości funkcji w węzłach interpolacji
       char *plkNazwaWezly, *plkNazwa; // nazwy pierwszego i drugiego pliku
       double aa[MAXN], bb[MAXN]; // tablice współczynników ai i bi we wzorze (24)




    double gfun(double y)
    {
       double x;
       // I
       if (y < ((b - a)/2))
          for(x=a; x<=((a+b)/2); x++)
          {
          return (x - a);
          }
       // II
       else if (y == ((b - a)/2))
                    for(x=((a+b)/2); x<=(a+(3*((b-a)/2)/4)); x++)
          {
          return (b - a)/2;
          }
       // III
       else if (y < (b - a)/2)
          for(x=(a+(3*((b-a)/2)/4)); x<=b; x++)
          {
          return (-2*(x - b));
          }
    }
    double g(double t)
    {
       double tt;
       tt=fmod(t-a, b-a);
       return gfun(tt<0?b+tt:a+tt);
    }
    // Wyznaczenie węzłów interpolacji i wartości funkcji w węzłach
    void generWezlyTryg()
    {
       double delta=(b-a)/(n+1);
       for (int i=0; i<=n; i++)
       {
          y[i]=a+delta*i;
          //x[i]=y2x(y[i]);
          f[i]=gfun(y[i]);
          printf("%f    %f \n", y[i], f[i]);
       }
       y[n+1]=b;
       m=n/2; //wyznaczenie współczynnika m ze wzorze (24)
       ddelta=n-2*m; // wyznaczenie współczynnika delta ze wzoru (24)
    }
    // Generacja parametrów ai oraz bi ze wzoru (24)
    void genAi_Bi()
    {
       double u1, u2, u, z;
       int k, j;
       z = 2 / ((double)n + 1);
       
       for(j = 0; j <= m+1; j++)
       {
              u2 = u1 = 0;       
              for(k = n; k >= 1; k--)
              {
                    u = f[k] + 2 * u1 * cos(x[j]) - u2;
                    u2 = u1;
                    u1 = u;
              }             
          aa[j] = z * (f[0] + u1 * cos(x[j]) - u2);
          bb[j] = z * u1 * sin(x[j]);      
       }

       return;
    }
    double Tryg(double t)
    {
       double tt;
       
       tt = 0.5 * aa[0];
       
       for(int j = 1; j <= m; j++)
       {
          tt += (aa[j] * cos(j * t)) + (bb[j] * sin(j * t));
       }
       tt += (0.5 * ddelta) * aa[m+1] * cos((m+1) * t);
       return tt;
    }
    // Wyprowadzenie danych do pierwszego pliku
    void generWezlyPlik()
    {
       FILE *plk;
       plk = fopen( plkNazwaWezly, "w" );
    // wydruk nagłówka dla pliku .rex
       fprintf(plk, "yi\t g(yi)\n");
    // wydruk trójek węzeł  na osi y, węzeł na osi x i wartość funkcji w węźle interpolacji
       for (int i=0; i<=n; i++)
       {
          fprintf(plk, "%lf\t  %lf\n", y[i],  f[i]);
       }
    // wydruk stopki dla pliku .rex
       fprintf(plk, "*InfKon*\nV 0 0 0 %d",n+1 );
       fclose(plk);
    }
    //Wyprowadzenie wiersza do drugiego pliku
    void drukWiersz(double ty, FILE *plk)
    {
       double eps=1.0e-20;
       double wyn, gg;
       wyn=Tryg(ty);
       gg=g(ty);
       fprintf(plk, "%lf\t %lf\t %lf\t %lf\t %lf\n", ty, gg, wyn, gg-wyn, (gg-wyn)/(fabs(gg)+eps));
    }
    // Wyprowadzenie danych do drugiego pliku y należy do przedziału <a, b>
    void generWyn()
    {
       double delta;
       int i, j;
       FILE *plk;
       plk = fopen( plkNazwa, "w" );
    // wydruk nagłówka dla pliku .rex
       fprintf(plk, "y\t g(y)\t tt(y)\t g(y)-tt(y)  (g(y)-wyn)/(abs(g(y))+eps)\n");
       for (i=0; i<=n; i++)
       {
          delta=(y[i+1]-y[i])/mn;
          for(j=0; j<mn; j++)
          {
             drukWiersz(y[i]+delta*j, plk);
          }
       }
       drukWiersz(b, plk);
    // wydruk stopki dla pliku .rex
       fprintf(plk, "*InfKon*\nV 0 0 0 %d", (n+1)*mn+1);
       fclose(plk);
    }
    // Wyprowadzenie danych do drugiego pliku y należy do przedziału <a-(b-a), b+(b-a)>
    void generWyn1()
    {
       double ty, delta;
       int i, j;
       FILE *plk;
       plk = fopen( plkNazwa, "w" );
    // wydruk nagłówka dla pliku .rex
       fprintf(plk, "y\t g(y)\t tt(y)\t g(y)-tt(y)  (g(y)-tt(y))/(abs(g(y))+eps)\n");
       for (i=0; i<=n; i++)
       { //y należy do przedziału <a+(b-a), a)
          delta=(y[i+1]-y[i])/mn;
          for(j=0; j<mn; j++)
          {
             ty=-(b-a)+y[i]+delta*j;
             drukWiersz(ty, plk);
          }
       }
       for (i=0; i<=n; i++)
       {
          delta=(y[i+1]-y[i])/mn;
          for(j=0; j<mn; j++)
          { //y należy do przedziału <a, b)
             ty=y[i]+delta*j;
             drukWiersz(ty, plk);
          }
       }
       for (i=0; i<=n; i++)
       { //y należy do przedziału <b, b+(b-a))
          delta=(y[i+1]-y[i])/mn;
          for(j=0; j<mn; j++)
          {
             ty=(b-a)+y[i]+delta*j;
             drukWiersz(ty, plk);
          }
       }
       //y=b+(b-a)
             drukWiersz(ty, plk);
    // wydruk stopki dla pliku .rex
       fprintf(plk, "*InfKon*\nV 0 0 0 %d", 3*(n+1)*mn+1);
       fclose(plk);
    }

    int main(int argc, char* argv[])
    {
    // Odczyt danych z linii polecenia
       a=atof(argv[1]); //lewa granica przedziału
       b=atof(argv[2]); //prawa granica przedziału
       n=atoi(argv[3]); //maksymalny indeks węzła interpolacji
       mn=atoi(argv[4]); //liczba pośrednich punktów między węzłami interpolacji w drugim pliku
       plkNazwa=argv[5];  //nazwa drugiego pliku (pliku z wynikami działania funkcji interpolującej)
       plkNazwaWezly=argv[6]; //nazwa pierwszego pliku (pliku z węzłami)
    // Wyprowadzenie na ekran danych z linii polecenia
       printf("a= %lf\n", a);
       printf("b= %lf\n", b);
       printf("n= %d\n", n);
       printf("mn= %d\n", mn);
    // Wyznaczenie węzłów interpolacji i wartości funkcji w węzłach
       generWezlyTryg();
    // Wyprowadzenie danych do pierwszego pliku
       generWezlyPlik();
    // Generacja parametrów ai oraz bi ze wzoru (24)
       genAi_Bi();
    // Wyprowadzenie danych do drugiego pliku y należy do przedziału <a, b>
       generWyn();
    // Wyprowadzenie danych do drugiego pliku y należy do przedziału <a-(b-a), b+(b-a)>
    //   generWyn1();
       system("pause");
       return 0;
    }



    Dane a, b, n i mn podajemy w parametrach w devie

    Z góry dzięki za szybką odpowiedź

    Kod proszę umieszczać w znacznikach code - krzychoocpp
  • #2
    mykhaylo
    Level 14  
    Jeśli wszystko się ładnie kompiluje, to nie znaczy że będzie działać. Wklejając masę kodu nie spodziewaj się, że ktoś to będzie analizował. Też kiedyś pisałem mnóstwo algorytmów i według mnie najlepszą metodą na znalezienie błędu jest debugowanie. Stawiasz breakpoint-a po czym śledzisz krok po kroku co się wykonuje i jakie przyjmuje wartości. Czasochłonne ale skuteczne.
  • #3
    torcek
    Level 10  
    Wszystko niby ok fajnie pięknie ale jedne co jest na pewno źle to w funkcji gfun() tak gdzie jest drugi przedział ze wzoru na funkcję przechodzącą przez 2 punkty wychodzi że gfun= (b-a)/2 dla y= (b-a)/2 bo y w tym jest miejscu jest stały, w forze daje krok dla x w jakim przedziale ma się zmieniać a w wynikach mam tylko wyniki z przedziału I i III. Nie wiem ja to zrobić żeby mi to dobrze policzył.

    to jest pierwszy problem drugi to to żeby wyskakuje mi okienko z wynikami razem z errorem i nie zapisuje mi wyników do pliku, i tu nie mam już zupełnie pojęcia dlaczego tak jest
  • #4
    mykhaylo
    Level 14  
    Jeśli chodzi o porównywanie double, to różnie to bywa. Może w taki sposób sprawdzaj czy obie wartości są równe:
    Code:

    bool doubleEquals(double left, double right, double epsilon) {
      return (fabs(left - right) < epsilon);
    }

    epsilon, jest to stała "tolerancji", jeśli obie wartości mają różnicę mniejszą od tej stałej, to uważamy że są równe.

    Co to za błąd wyskakuje z okienkiem z wynikami. Ma jakąś konkretną treść?
  • #5
    torcek
    Level 10  
    "Wystąpił problem z aplikacją Projekt.exe i zostanie ona zamknięta" jest to błąd windowsowski
    Sygnatura błędu:
    AppName: projekt1.exe AppVer: 0.0.0.0 ModName: ntdll.dll
    ModVer: 5.1.2600.5512 Offset: 00001278
  • #6
    mykhaylo
    Level 14  
    Wydaje mi się że są jakieś błędy z operacjami na pamięci. Z tego powodu ten program może się wywalać. Musisz wykorzystać opcję Debug. Innej rady nie mam.
  • #7
    torcek
    Level 10  
    Słuchaj jak używam debugera to wyskakuje mi taki warning:
    "W twoim programie zgłoszono naruszenie praw dostępu (błąd segmentacji)"

    i na czerwono mam zaznaczoną linijkę gdzie jest warunek że y ma być równe ((b-a)/2

    Wiecie może co z tym zrobić ?
  • #8
    Dr.Vee
    VIP Meritorious for electroda.pl
    1) kompiluj bez optymalizacji
    2) upewnij się, że podajesz wszystkie 6 parametrów do programu

    Pozdrawiam,
    Dr.Vee
  • #9
    torcek
    Level 10  
    Podaje tylko 4 wartości:
    początek przedziału i koniec czyli a i b
    liczbę węzłów n
    i licze punktów pośrednich mn

    W kodzie jest mały błąd bo tam jest przypisane do n i m wartość która podaje w parametrach i do m nie mam nic przypisywać
  • #10
    Dr.Vee
    VIP Meritorious for electroda.pl
    Powinieneś podać 6 parametrów - dodatkowo 2 nazwy plików.

    Pozdrawiam,
    Dr.Vee
  • #11
    torcek
    Level 10  
    Zaczęło nawet działać :)
    Zmieniłem wszystko w funkcji gfun na
    Code:

    double gfun(double y)
    {
     double t;
     for (t=0; t<=(b-a)/2; t++)
    {
      if (y<(a+b)/2)     // I
       return (t-a);
      else if (y<(a+(3*((b-a)/2)/4)))     // II
       return (b-a)/2;
      else
        return (-2*(t-b));
    }
    }

    reszta pozostała taka sama

    Wszystkie wyniki w miare ok ale nie wiem dlaczego dla funkcji f[i] dla 6 węzłów wyświetla mi wynik 2 2 2 4 4 4 przedział od -2 do 2

    Dodano po 22 [minuty]:

    Dobra panowie gfun poprawiona ponownie i wszystko wylicza poprawnie o prócz y który liczy mi do 5 a nie do 2 przy założeniu że a=-2 i b=2

    ============

    Panowie już ostatnie pytanie program działa bez problemu
    Code:

    if (x < (a+b)/2)                 
             return (x - a);             
       else if (x < (a+(3*((b-a)/2)/4)))
             return (b - a)/2;
       else
           return (-2*(x - b));

    i teraz mam tutaj problem
    czyli funkcja w przedziale miedzy (a+b)/2 a (a+(3*((b-a)/2)/4)) ma być linią prostą (równoległą do osi X :) ) czyli tj równą 1 a mi ona leci hen do góry
    Czy mam w tym else if dac jeszcze jeden warunek tj =1 czy cos w tym stylu?

    Proszę nie pisać postu pod postem tylko używać funkcji ZMIEŃ.
    Posty scaliłem.
    [Dr.Vee]


    No i jak panowie dacie jakiś pomysł żeby funkcją w tym przedziale między (a+b)/2 a (a+(3*((b-a)/2)/4)) czyli to gdzie jest else if była linią równoległą do x.
    Programik mi działa ok wyświetla mi dobrze punkty w I i III przedziale ale w II jest źle ponieważ punkty które powinny tworzyć linię równoległa do osi x idą do góry. Wykres ma wyglądać mniej więcej tak linia się wznosi , jest równoległa do x i na koniec opada

    Macie jakieś sugestie ?