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