W dniu 2012-11-11 17:58, Baranosiu pisze:
> 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
>>

Hmm... Nawet z tak prostą operacją jak dodanie słupka iloczynów
nie dał sobie rady? Niemożliwe. Chociaż...

> [...]
> 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.

Sławek schrzaniał algorytm i tyle. Pewnie źle dobrał
epsilon maszynowy;)



Wrzucam porównanie kilku podstawowych metod.
Newtona_Cotesa do 6 punktowych,
(za stopień 0 robi 'midpoint', dwupunktowa to trapezy,
trzypunktowa - simpson), kwadratury Gaussa 2-6 punktowe
i kwadratury Lobatto, (czyli takie, które mają wezły
na brzegach, ale pozostałe są tak dobrane, by wyciągnąć
największy rząd) operte na 4 do 7 węzłów.
(znów 2 to trapez, 3 to simpson).

Węzły N-C są równoodległe, pozostałe to pierwiastki pewnych
wielomianów. Wielomiany te wyznaczane są numerycznie z definicji
rekurencyjnej. Tak samo są wyszukiwane pierwiastki. Wszystko
za pomocą double, dlatego też nie można przesadzać z rzędem
metody - prowadzi to do błędów, o których potem.

Wagi dobierane są wprost z definicji. W_j to całka na rozpatrywanym
przedziale z wielomianu który zeruje się na wszystkich węzłach,
poza x_j, gdzie przyjmuje wartość 1. Całka z wielomianu liczona
algebraicznie;)

Następnie tworzymy węzły i wagi dla kwadratury złożonej z określonej
liczby powtórzeń na przedziale do całkowania. Jeśli podstawowa
kwadratura miała węzły na obu brzegach, w złożonej będą one zespolone
w jeden (o wadze równej sumie).


Takie podejście umożliwia testowanie wielu różnych metod bez rozbijania
się na przypadki.

Na wykres wrzucamy błąd względny w funkcji ilości wywołań
funkcji podcałkowej (nie ilości przedziałów, a rzeczywista
ilość wywołań). Skala log log.
Wartość dokładna wyliczona przez wolfrem alpha;)


Wątkowy przykład, sin(x) exp(-x) na [0,5]
http://w18.wrzuta.pl/obraz/70M5xduXyOU/sin_exp

Wersja zagęszczona, sin(10*x) exp(-x) na [0,5]
http://w18.wrzuta.pl/obraz/0SDBiKFZnw0/sin10_exp


Co widać?

Simpson i 3/8 oraz midpoint i trapezy czy metody N-C 5 i 6 punktowe
sa tego samego rzędu.
Oczywiście pierwsza para zbiega znacznie szybciej. Dokładność
20 punktowego złożonego Simpsona trapez osiąga dla 200 punktów.

Coś więcej o tym napisać? Niska kwadratura oczko wyższa jest
oczko lepsza. Dość oczywiste, eksperyment potwierdził, tyle.


Dla tego przykładu "lżejsze" metody z każdej pary
mają ciut lepszą stałą. W sumie nie takie dziwne, podobne
oszacowanie na błąd, a rozkłada się on na 3 lub 4 punkty.

Co jeszcze? Dwupunktowy Gauss działa jak N-C na 3 i 4 punktach,
teoria znów górą.

O, i ciekawostka nawiązująca do niedawnego wątku o dokładności
numerycznej. Zamiast pracowicie przepisywać współczynniki
węzłów i wag zbudowałem maszynkę, która je wyznacza,
a to obejmowało rekurencyjne wyznaczanie współczynników
wielomianów (odejmowanie!) i wyliczanie wartości tych
wielomianów (schemat hornera, ale znów możliwa strata
dokładności przez odejmowania) wyznaczone wagi mogą być
obarczone błędem większym niż dokładność numeryczna double.
Tym większym, im wyższy rząd.

I doskonale widać to na wykresie. Metody oparte na większej
liczbie punktów zbiegają szybciej, ale wynik stabilizuje się
na gorszej dokładności.


Kody
https://skydrive.live.com/redir?resid=25C9E1E8909658
1!2201&authkey=!AGMwPjeQhRXoDe0


A jak już ciekawostkami. Była taka funkcja użyta
przez Rungego do pokazania kłopotów z interpolacją
wielomianową:
http://en.wikipedia.org/wiki/Runge%27s_phenomenon

No to hop:
testuj_calke(@(x)1./(1+25*x.^2),-1,1,0.5493603067780
0634434450877,20,500)

http://w18.wrzuta.pl/obraz/8lIeKjKRUs4/runge

:-)



pzdr
bartekltg