Szereg Taylora

Odpowiedz Nowy wątek
2011-08-16 18:58
0

Witam,
Mam problem zakodawałem funkcję arctan rozwiniętą w szereg Taylora, poniżej kod

 #include <stdio.h>
#include <math.h>
 
int sgn(double x) {
    if (x > 0)
        return 1;
    else if (x < 0)
        return -1;
    else
        return 0;
}
 
double arctg(double x) {
    double suma, wyraz;
    int znak, n = 1;
 
    if (fabs(x) <= 1) {
        suma = 0;
        wyraz = x;
        znak = -1;
 
        do {
            n=1;
            suma = suma + wyraz;
            wyraz = wyraz * znak * ((2 * n - 1) * x * x) / (2 * n + 1);
            znak = znak * (-1);
            n++;
        } while (abs(suma) > 1e-7);
 
        return suma;
    }
 
    else {
        suma = sgn(x) * (M_PI / 2);
        wyraz = -1 / x;
        znak = -1;
 
        do {
            suma = suma + wyraz;
            wyraz = wyraz * znak * 1 / (((2 * n - 1) * x * x) / (2 * n + 1));
            znak = znak * (-1);
            n++;
        } while (fabs(wyraz) > 1e-7);
        return suma;
    }
 
}
 
int main(void) {
    double x = -10;
 
    printf("%10s %10s %10s %10s\n","x","arctg(x)","atan(x)","error");
 
    while (x <= 10) {
        printf("%10f %3.7f %3.7f %3.7f\n", x, arctg(x), atan(x), fabs(arctg(x) - atan(x)) / fabs(atan(x)));
        x += 1;
    }
 
    getchar();
    getchar();
    return 0;
}
 

Tak wygląda rozwinięcie dla |x|<= 1 user image a tak dla |x| > 1 sgn(x)*Pi/2 - 1/x + 1/3x^2 - 1/5x^5+...
Problem polega na tym, że błąd wychodzi za duży, powinno być co najmniej 5 zer a dopiero później jakieś cyfry.
Proszę o pomoc.

2011-08-16 20:29
0

Przyglądnąłem się tylko pierwszej części funkcji ( tej gdzie fabs(x) <= 1 ) i myślę, że warunek abs(suma) > 1e-7 jest błędny.
Ponieważ oczekujesz, że suma co do wartości bezwzględnej będzie się zbliżać do zera, a na przykład:
arctg(0.7) = ~0.61 i ile razy byś nie iterował to nie spadnie poza 1e-7 -> Twój program się zapętla.
Powinieneś też użyć fabs, a nie abs.

Ah i bym zapomniał, sądzę iż w tej lini
wyraz = wyraz znak ((2 n - 1) x x) / (2 n + 1);
niepotrzebne jest (2*n-1), nie ma go przynajmniej we wzorze z wiki, oczywiście mogę się mylić ;)

Myliłem się ;)

Napisałem tą część funkcji po swojemu i mam dobre przybliżenia:
http://pastebin.com/s574uEyr

edytowany 2x, ostatnio: Atael, 2011-08-16 21:29

Pozostało 580 znaków

2011-08-16 21:02
0
Atael napisał(a)

Powinieneś też użyć fabs, a nie abs.

Ah i bym zapomniał, sądzę iż w tej lini

Jak dam fabs to kod siada ...

Atael napisał(a)

wyraz = wyraz * znak * ((2 * n - 1) * x * x) / (2 * n + 1);
niepotrzebne jest (2*n-1), nie ma go przynajmniej we wzorze z wiki, oczywiście mogę się mylić ;)

Zauważ, że następny wyraz powstaję przez pomnożenie starego i trzeba skrócić, to co z niego nie pasuje.
Teraz podaje pełne zadanie.
http://www.image-share.com/ipng-861-44.html

edytowany 2x, ostatnio: mrfustrr, 2011-08-16 21:29
Tylko u mnie obrazka nie widać? - Afish 2011-08-16 21:23
Udostępniłem w innym sposób obrazek, powinno działać. - mrfustrr 2011-08-16 21:30
Prawda, mój błąd. Aczkolwiek abs -> działa na intach a fabs działa na float/double, więc powinieneś użyć fabs. Dalej uważam, że sam warunek jest zły. Powinieneś patrzeć na wyraz, a nie na całą sumę. - Atael 2011-08-16 21:30
Dziwi mnie to trochę bo mam przykład dla area tangesa napisany w ten sposób i działa prawidłowo. - mrfustrr 2011-08-16 21:33

Pozostało 580 znaków

2011-08-16 22:03
0
const double epsilon = 1.e-14;
 
double myATan(double x)
{
    if (x<-1)
        return -M_PI*0.5 - myATan(1.0/x);
    if (x>1)
        return M_PI*0.5 - myATan(1.0/x);
    double sum(x), a(x), x2(x*x);
    for (int n=3; a/n>epsilon || a/n<-epsilon; n+=2) {
        a*=-x2;
        sum += a/n;
    }
    return sum;
}

Edit: po małej poprawce działa mi bez zarzutu (testy przechodzą).


Jeśli chcesz pomocy, NIE pisz na priva, ale zadaj dobre pytanie na forum.
edytowany 4x, ostatnio: MarekR22, 2011-08-16 23:18
Dla -10.0 daje błąd ponad 2.0, chyba nie tak miało to działać? Na początku myślałem, że dla -1.0 się zawiesza, ale działa po prostu super długo : d - Atael 2011-08-16 22:20
zapominałem przetestować ujemne :), kolejna poprawka i działa testowane dla argumentów: 0.0, 0.2, 0.25, 0.5, 0.75, 0.9, 2, 4, 10, 20, 30, -0.2, -0.25, -0.5, -0.75, -0.9, -2, -4, -10, -20, -30 - MarekR22 2011-08-16 23:05
Błędy są już mniejsze, ale odpal u siebie dla -1, też Ci tak muli ? - Atael 2011-08-16 23:16
to akurat wada samego algorytmu. Ten szereg NIE jest zbiezny dla 1 (suma jest nieskonczona)! - MarekR22 2011-08-17 00:15
oh, tak to jest jak się bez myślenia odpala kod ;x zwracam honor - Atael 2011-08-17 01:09

Pozostało 580 znaków

2011-08-16 22:16
0
MarekR22 napisał(a)
 
const double epsilon = 1.e-10;

double myATan(double x);
{
if (x<-1)
return M_PI0.5 + myATan(-1.0/x);
if (x>1)
return M_PI
0.5 - myATan(1.0/x);
double sum(x), a(x), x2(xx);
for (int n=1; a/n>epsilon || a/n<-epsilon; n+=2) {
a
=-x2;
sum += a/n;
}
return sum;
}



Dziwny kod. Ten kawałek, to `double sum(x), a(x), x2(x*x);` rozumiem, że deklaracja trzech funkcji double bez podawania typów argumentów.
`a*=-x2;` w tej linijce kompilator wywala błąd składni.
Pokaż pozostałe 3 komentarze
Całkowicie zapomniałem, że piszesz w czystym C. - Atael 2011-08-16 22:29
w moim kodzie był mały błąd (było n=1 zmiast n=3), po przeprowadzeniu testów, wyniki dostaje poprawne. - MarekR22 2011-08-16 22:45
a=-x2; znaczy dokładnie to samo co: a = -a x2; - MarekR22 2011-08-16 22:46
Skopiowałem Twój kod, przebudowałem projekt, ale dalej mam złe wyniki. Tak samo dla -10 daje mi błąd ponad 2.0 a dla -1 działa bardzo długo, chyba coś robię nie tak oO - Atael 2011-08-16 22:52

Pozostało 580 znaków

2011-08-16 22:57
0

Jednak nie działa poprawnie:

         x   arctg(x)    atan(x)      error
-10.000000 1.6704650 -1.4711277 2.1354997
 -9.000000 1.6814535 -1.4601391 2.1515708
 -8.000000 1.6951513 -1.4464413 2.1719461
 -7.000000 1.7126934 -1.4288993 2.1986103
 -6.000000 1.7359450 -1.4056476 2.2349788
 -5.000000 1.7681919 -1.3734008 2.2874551
 -4.000000 1.8157750 -1.3258177 2.3695511
 -3.000000 1.8925469 -1.2490458 2.5151942
 -2.000000 2.0344439 -1.1071487 2.8375525

dla kodu:

#include <stdio.h>
#include <math.h>
 
int sgn(double x) {
    if (x > 0)
        return 1;
    else if (x < 0)
        return -1;
    else
        return 0;
}
 
const double epsilon = 1.e-10;
 
double arctg(double x) {
    if (x < -1)
        return M_PI * 0.5 + arctg(-1.0 / x);
    if (x > 1)
        return M_PI * 0.5 - arctg(1.0 / x);
    double sum = x, a = x, x2 = x * x;
    int n;
    for (n = 3; a / n > epsilon || a / n < -epsilon; n += 2) {
        a *= -x2;
        sum += a / n;
    }
    return sum;
}
 
int main(void) {
    double x = -10;
 
    printf("%10s %10s %10s %10s\n", "x", "arctg(x)", "atan(x)", "error");
 
    while (x <= 10) {
        printf("%10f %3.7f %3.7f %3.7f\n", x, arctg(x), atan(x),
                fabs(arctg(x) - atan(x)) / fabs(atan(x)));
        x += 1;
    }
 
    getchar();
    getchar();
    return 0;
}
 
edytowany 2x, ostatnio: mrfustrr, 2011-08-16 23:02
To jest tylko pomyłka w znakach, nie zauważyłeś, że arctg(x)=atan(x)+pi. Spróbuj to poprawić bez pomocy. - Zjarek 2011-08-16 23:15
nie rozumiem - mrfustrr 2011-08-17 14:35

Pozostało 580 znaków

2011-08-17 18:56
0

Ja mam taka uwagę do niektórych powyższych kodów: nie wiem jak tam c, ale w c++ mieszanie intów i floatów zazwyczaj prowadzi do dziwnych błędow, które potem ciężko znaleźć. :-P

Poniższy kod (C++) działa bardzo dobrze:

double arctg_taylor(double x, double eps) {
  if (std::fabs(x) >= 1) {
    throw std::runtime_error("tak se ne da");
  }
 
  double old_term = 0.0, term = x, sum = x, i = 1.0;
 
  do {
    old_term = term;
 
    term = -term * x * x;
 
    sum += (1 / (i + i + 1)) * term;
 
    i += 1.0;
  } while (std::fabs(old_term + term) > eps);
 
  return sum;
}

Kiedy zmieni się i na int (co wydaje się ok, to tylko licznik pętli) wyniki są nieprawidłowe.


"(...) otherwise, the behavior is undefined".
edytowany 1x, ostatnio: Endrju, 2011-08-17 18:57
Prawdopodobnie ma to związek z reprezentacją double w pamięci i ciągłej konwersji z int na double, ale to są tylko moje przypuszczenia. Nigdy nie musiałem zgłębiać zagadnień związanych z błędami zaokrągleń ;) - Atael 2011-08-17 19:06
Nie ma tu zadnej tajemnicy: winne jest (1 / (i + i + 1)). Kiedy i jest int, dla kazdego i > 0 wyraz ten bedzie rowny 0, co jest po prostu wynikiem obciecia czesci ulamkowej. - Endrju 2011-08-17 19:27
O racja, nie zwróciłem na to uwagi ;x - Atael 2011-08-17 19:51

Pozostało 580 znaków

Odpowiedz
Liczba odpowiedzi na stronę

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