W kontekście rozmaitych forumowych "specjalistów" od metod numerycznych -
zabawne, jak tym razem będą wyjaśniać trywialny w istocie rezultat...

exact 0.5022749400837604
tapez 0.5022749194207212
simpson 0.5022738634614975

/*
* tq.c
* 2012 (C) slawek
*/

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

#define NPTS 10000 /* NPTS must be even! */

double x[NPTS+1];
double y[NPTS+1];

double f(double x)
{
return sin(x)*exp(-x);
}

double integrate_exact(double a, double b)
{
return ((cos(a) + sin(a))*exp(-a) - (cos(b) + sin(b))*exp(-b))/2.0;
}

double generate_points(double x[], double y[], double a, double b, int n)
{
int i;

for(i = 1; i <= n; i++)
{
x[i] = a + (i-1) * (b - a)/(n - 1);
y[i] = f(x[i]);
}
}

double integrate_trapez(double x[], double y[], int n)
{
int i;
double s = 0.;

for(i = 2; i <= n; i++)
s += (y[i] + y[i-1]) * (x[i] - x[i-1]) * 0.5; // an naive approach

return s;
}

double integrate_simpson(double x[], double y[], int n)
{
int i;
double odd = 0., even = 0.;

for(i = 1 ; i < n ; i+=2) odd += y[i];
for(i = 2 ; i <= n ; i+=2) even += y[i];

return (2.0*odd-y[1]-y[n] + 4.0*even)*(x[3]-x[1])/6.0;
}

int main(int argc, char* argv[])
{
const double a = 0.;
const double b = 5.;

generate_points(x,y,a,b,NPTS);
printf("exact %.16lg\n", integrate_exact(a,b));
printf("tapez %.16lg\n", integrate_trapez(x,y,NPTS));
printf("simpson %.16lg\n", integrate_simpson(x,y,NPTS));

puts("press enter...");
getchar();
return 0;
}

do góry

Użytkownik "slawek" <s...@h...pl> napisał w wiadomości grup
dyskusyjnych:509ee300$0$26682$6...@n...neostrad
a.pl...
> W kontekście rozmaitych forumowych "specjalistów" od metod numerycznych -
> zabawne, jak tym razem będą wyjaśniać trywialny w istocie rezultat...
>
> exact 0.5022749400837604
> tapez 0.5022749194207212
> simpson 0.5022738634614975

Autopoprawka - dorzuciłem jeszcze dwa algorytmy - całkowanie spline'em i
jeszcze jeden algorytm (ale nie wyjaśniam jaki, ot "niespodzianka").

exact 0.5022749400837604
trapez 0.5022749194207212 -0.04114 ppm
simpson 0.5022738634614975 -2.14349 ppm
spline 0.5022749607347478 +0.04111 ppm
best 0.5022749400777344 -0.00001 ppm

Dla łatwiejszego porównania - podane są względne (parts per milion)
odchylenia od wyniku analitycznego.

Co widać? Ano, metoda Simpsona jest NAJGORSZA i daje w porównaniu ze
zwykłymi trapezami odchylenie około 50 razy większe.

W porównaniu z "thebestciakiem" Simpson nie ma szans - thebeściak jest o
pięć rzędów dokładniejszy - i to chyba najlepiej podsumowuje wasze "miszcze"
doświadczenia z numerycznym całkowaniem.

Trochę inne dane - znowu to samo - metoda Simpsona daje najmniej dokładny
wynik.

exact 0.4993567179497454
trapez 0.4993567163164947 -3.3e-003 ppm
simpson 0.4993567459083201 +5.6e-002 ppm
spline 0.4993567195827356 +3.3e-003 ppm
best 0.4993567179496152 -2.6e-007 ppm

do góry

W dniu sobota, 10 listopada 2012 18:28:01 UTC-5 użytkownik slawek napisał:

> s += (y[i] + y[i-1]) * (x[i] - x[i-1]) * 0.5; // an naive approach

Ja w kwestii pobocznej: "an" zamiast "a" i "the" wymawiane "thi" zamiast
"the" istnieje w jezyku angielskim dla oddzielenia samoglosek:
"an approach", "a naive approach". To dla nas trudne, bo nie mamy
w jezyku polskim tych slowek i nigdy sie nie wie czy "a" czy "the" czy nic.

--
Edek

do góry

W dniu niedziela, 11 listopada 2012 14:33:24 UTC+1 użytkownik slawek napisał:
> U�ytkownik "slawek" <s...@h...pl> napisa� w wiadomo�ci grup
>
> dyskusyjnych:509ee300$0$26682$6...@n...neostrad
a.pl...
>
> > W kontek�cie rozmaitych forumowych "specjalist�w" od metod numerycznych -
>
> > zabawne, jak tym razem b�d� wyja�nia� trywialny w istocie rezultat...
>
> >
>
> > exact 0.5022749400837604
>
> > tapez 0.5022749194207212
>
> > simpson 0.5022738634614975
>
>
>
> Autopoprawka - dorzuci�em jeszcze dwa algorytmy - ca�kowanie spline'em i
>
> jeszcze jeden algorytm (ale nie wyja�niam jaki, ot "niespodzianka").
>
>
>
> exact 0.5022749400837604
>
> trapez 0.5022749194207212 -0.04114 ppm
>
> simpson 0.5022738634614975 -2.14349 ppm
>
> spline 0.5022749607347478 +0.04111 ppm
>
> best 0.5022749400777344 -0.00001 ppm
>
>
>
> Dla �atwiejszego por�wnania - podane s� wzgl�dne (parts per milion)
>
> odchylenia od wyniku analitycznego.
>
>
>
> Co wida�? Ano, metoda Simpsona jest NAJGORSZA i daje w por�wnaniu ze
>
> zwyk�ymi trapezami odchylenie oko�o 50 razy wi�ksze.
>
>
>
> W por�wnaniu z "thebestciakiem" Simpson nie ma szans - thebe�ciak jest o
>
> pi�� rz�d�w dok�adniejszy - i to chyba najlepiej podsumowuje wasze
"miszcze"
>
> do�wiadczenia z numerycznym ca�kowaniem.
>
>
>
> Troch� inne dane - znowu to samo - metoda Simpsona daje najmniej dok�adny
>
> wynik.
>
>
>
> exact 0.4993567179497454
>
> trapez 0.4993567163164947 -3.3e-003 ppm
>
> simpson 0.4993567459083201 +5.6e-002 ppm
>
> spline 0.4993567195827356 +3.3e-003 ppm
>
> best 0.4993567179496152 -2.6e-007 ppm

a jaka to wlasnie metoda? ja taka 'numeryką' się nie zajmuje (i pewnie moze nigdy nie
bede zajmowac) ale to ciekawe co tak dobrze dziala - ja bym pewnie robil wlasnie
trapezami, chyba ze jest kjakis trick ktory dobrze dziala

do góry

On Sun, 11 Nov 2012 07:42:42 -0800 (PST), e...@g...com
wrote:
> "the" istnieje w jezyku angielskim dla oddzielenia samoglosek:

Nie zawsze - ale to ntg.

do góry

Dnia 10.11.2012 slawek <s...@h...pl> napisał/a:
> W kontekście rozmaitych forumowych "specjalistów" od metod numerycznych -
> zabawne, jak tym razem będą wyjaśniać trywialny w istocie rezultat...
>
> exact 0.5022749400837604
> tapez 0.5022749194207212
> simpson 0.5022738634614975
>
[...]
Akurat w tym przypadku tak, ale weź inny przedział całkowania (na
przykład -10..10 tak żeby wpływ funkcji wykładniczej był nieco
bardziej znaczący) i już simpson może wypaść lepiej.

do góry

Użytkownik "slawek" <s...@h...pl> napisał:

> Trochę inne dane - znowu to samo - metoda Simpsona daje najmniej dokładny wynik.
>
> exact 0.4993567179497454
> trapez 0.4993567163164947 -3.3e-003 ppm
> simpson 0.4993567459083201 +5.6e-002 ppm
> spline 0.4993567195827356 +3.3e-003 ppm
> best 0.4993567179496152 -2.6e-007 ppm

Dorzuc jeszcze wiecej !
Po prostu osmieszasz sie doszczetnie czlowieku :)

PS: Ze juz nie wspomne o "jakosci" kodu (a raczej nacpanego
antykodu) ktory bez wstydu zamiesciles.

AK

do góry

Użytkownik "Baranosiu" <r...@w...pl> napisał:

> [...]
> Akurat w tym przypadku tak, ale weź inny przedział całkowania (na
> przykład -10..10 tak żeby wpływ funkcji wykładniczej był nieco
> bardziej znaczący) i już simpson może wypaść lepiej.

Glownie nawet nie o to chodzi.
Chodzi glownie o zbieznosc (a wiec kosz ale i dokladnosc).
PS: Dla 1000 000 punktow wszyskie metody dadza juz jeszcze bardziej do siebie
zblizone wyniki :). o ale dla slawusie liczenie funkcji jest bezkosztowe wiec
w/g niego
hulaj dusza :) Dlaczego "tylko" 10 000 a ie 1 000 000 ? A co ! :)
PS1: Odwoluje nieniejszym me bledne mniemanie ze slawek to faktycznie gosc z
uczelni i kiepski ale jednak numeryk.
Jestem przekonany ze ani z uczelnia (no chyba ze tak bardzo juz zpsialy
uczelbie polskie?)
ani z numeryka ma malo wspolnego. Ot dostal ten neiwatpliwy juz dla mnie
palant
palant w rece Fortran i mysli, ze to starczy aby szpanowac na "merytorycznego"
:))
PS2: Masz racje. Przedzial i funkcje tez sobie dobral "pod przyklad"

AK

do góry

Użytkownik "Baranosiu" <r...@w...pl> napisał w wiadomości grup
dyskusyjnych:k7olf5$rpm$...@n...task.gda.pl...
> Akurat w tym przypadku tak, ale weź inny przedział całkowania (na
> przykład -10..10 tak żeby wpływ funkcji wykładniczej był nieco
> bardziej znaczący) i już simpson może wypaść lepiej.

Oczywiście, że dla /pewnych/ przedziałów lub /pewnych/ funkcji może być
tak... albo może tak nie być.

Jednakże mit o wyższości metody Simpsona nad metodą trapezów jest obalony -
nie można a priori założyć, że wyniki otrzymane metodą Simpsona będą
dokładniejsze.

I jeszcze drobiazg - większość ludzi, jak usłyszy "całkowanie", to kojarzy
to z zadaną w postaci w wzoru funkcją podcałkową. W przykładowym programie
taka funkcja, f, była wyłącznie dla niezaśmiecania forum tablicą parunastu
tysięcy wartości. Bo istotą rzeczy jest - w tym do czego mi są potrzebne
całki - że są dane pary (x,y), nie ma jawnie postaci funkcji. Można to sobie
np. wyobrazić jako zapis PCM dźwięku - przy próbkowaniu 44 kHz będzie to
44000 punktów (x,y) na każdą sekundę - i teraz trzeba to scałkować - dane
są - funkcji zapisanej wzorkiem nie ma.

Teraz trochę teorii. Jeżeli jest dość dużo sampli, czyli par (x,y) jest n+1,
a to n+1 wynosi kilkaset tysięcy lub więcej, to nie ma znaczenia, czy urwie
się jeden z sampli gdzieś z przodu lub z tyłu - wartość całki powinna być
/niemal/ taka sama. Stosując wzór Simpsona dla danych "urwanych z przodu"
otrzymuje się zupełnie inne sumowanie niż dla "urwanych z tyłu" (inne
elementy są są wtedy parzyste/nieparzyste). Można potem dodać "urwany"
przedział (ale to nie jest konieczne). Oczywiście dla otrzymania
dokładniejszych wyników można i należy wziąć średnią arytmetyczną wyników
całkowania "urwanego z przodu" i "urwanego z tyłu". Tylko że, surprise,
surprise, dostaje się wtedy po prostu wzór trapezów (kwadraturę
Newtona-Cotesa rzędu 1). Dlaczego tak jest? Między innymi dlatego, że wzór
na całkowanie numeryczne musi mieć "symetrię translacyjną" - tj. zastąpienie
przedziału indeksów [1, n] przez [2, n+1] powinno dawać ten sam wynik z
dokładnością do końcówek [1,2] i [n, n+1]. Tego warunku wzór Simpsona nie
spełnia, a wzór trapezów tak.

Kwadratury Newtona-Cotesa, Gaussa-Czebyszewa itd. powstawały w czasach
"przedkomputerowych". Służyć miały rozwiązaniu problemu: mamy funkcję f i
przedział [a, b] - jak tylko kilkakrotnie ewaluując f obliczyć wartość całki
z f, wiedząc że f jest taka a taka? Kilkakrotnie, tzn. 3-5 razy, bo każde
wyliczenie f trwa długo - papier i ołówek, ew. arytmometr (ale to już raczej
nie czasy Newtona i Gaussa, tylko trochę później).

Jeszcze nieco o całkowaniu metodą Romberga - to nic innego niż zastosowanie
ekstrapolacji Richardsona do całkowania. A ta jest, podobnie jak
ekstrapolacja Aitkena, należy do metod w których rezultaty kilku, coraz
dokładniejszych, obliczeń ekstrapoluje się do wyniku jaki prawdopodobnie
otrzymanoby przy jeszcze większej staranności. Praktycznie, w przypadku
danych stablicowanych oznaczałoby to policzenie metodą trapezów dla co
drugiego punktu (krok 2h) i normalnie (h = x[2] - x[1]), a następnie próbę
odgadnięcia ile to miałoby być dla h/2. Trzeba dodatkowego założenia o
odstępstwie wartości oszacowanych numerycznie od dokładnej. Tymczasem dla
danych stablicowanych takie założenie /może/ nie być prawdziwe - oraz co też
istotne - nie da się obliczyć w razie potrzeby np. wartości f((x[k] +
x[k+1])/2), gdyż formuła dla funkcji f nie jest znana, a interpolacji nie
należy stosować (tzn. można, lecz nie da ona nowych informacji, więc nie ma
sensu robić interpolacji w danych skoro efektywniej jest robić to dla sum).

Praktycznie? W niezłej bibliotece CERNLIB jest TRAPER (D108, Trapezoidal
Rule Integration with an Estimated Error), który - jak pamiętam - wbrew swej
nazwie używa interpolacji parabolicznej pomiędzy dostarczonymi odciętymi.
Lepiej jest używać używać funkcji sklejanych z wielomianów m-tego stopnia,
ale wymaga to rozwiązywania równania trójdiagonalnego, a to przy dużej
ilości danych może prowadzić do dużych błędów z zaokrągleń. Wracamy do
punktu wyjścia - metody Newtona-Cotesa rzędu 1.

Ogólnie sytuacja jest dość nieciekawa i to w tak trywialnie prostych
zagadnieniach, jak obliczanie RMS sygnału audio. Nic lepszego niż trapezy, a
w zasadzie nawet i to nie - bo z wzoru na trapezy wychodzi zwykłe sumowanie
wszystkiego co jest w środku i jeszcze doliczenie tylko połowy końcówek.

slawek

do góry

Użytkownik "AK" <n...@n...com> napisał w wiadomości grup
dyskusyjnych:k7oo6p$3ut$...@n...task.gda.pl...
> PS: Dla 1000 000 punktow wszyskie metody dadza juz jeszcze bardziej do
> siebie
> zblizone wyniki :). o ale dla slawusie liczenie funkcji jest
> bezkosztowe wiec w/g niego
> hulaj dusza :) Dlaczego "tylko" 10 000 a ie 1 000 000 ? A co ! :)

Oczywiście że jest "bezkosztowe" - są np. dane z przetwornika analog-cyfra -
czyli nie są liczone, ale pobierane z zewnątrz. Przy próbkowaniu 44 kHz to
88000 sampli samo przychodzi co sekundę kanałem lewym i prawym. Jakieś
problemy ze słuchem?

> ani z numeryka ma malo wspolnego. Ot dostal ten neiwatpliwy juz dla
> mnie palant
> palant w rece Fortran i mysli, ze to starczy aby szpanowac na
> "merytorycznego" :))

Ciekawe jak nazwać takiego kogoś jak AK, który udaje jakiegoś specjalistę, a
nie potrafi odróżnić kody źródłowego programu w języku C od programu w
języku Fortran?

> PS2: Masz racje. Przedzial i funkcje tez sobie dobral "pod przyklad"

Złość bierze? Coś nie tego z rozumieniem na czym polega dedukcja? Przecież
wystarczy dowolny przykład - a skoro dowolny - to może być i taki, który
sobie starannie wybiorę.

do góry