W dniu 2012-11-13 11:53, Roman W pisze:
> W dniu wtorek, 13 listopada 2012 10:33:09 UTC użytkownik slawek napisał:
>> Użytkownik "bartekltg" napisał w wiadomości grup
>> dyskusyjnych:k7rnav$8eq$...@n...news.atman.pl...
>>
>>> Sławek schrzaniał algorytm i tyle. Pewnie źle dobrał
>>> epsilon maszynowy;)
>>
>> A konkretnie - poza błędem w komentarzu?

A czemu mnie pytasz, plonka dostałeś.
Najpierw zachowujesz się jak cham i prostak,
a teraz chcesz, by szukać błędów w twoim kodzie?
Skopałeś implementacje, pokazało ci to kilka osób,
debuguj sam.


> Poprawnie zaimplementowany Simpson (wg.
http://en.wikipedia.org/wiki/Simpson%27s_rule#Compos
ite_Simpson.27s_rule) daje wartosc, dla Twoich danych, 0.5022749400837603, czyli
znacznie blizej prawdziwej wartosci niz Twoje
> trapezy.

To dorzucę jeszcze jeden obrazek, który dobitnie porównuje
różne rzędy kwadratur. Tym razem pomęczyłem komputer:

http://w26.wrzuta.pl/obraz/aBYTSD3FT4A/kwadratury_sz
eroko

Hi-Res:
http://c.wrzuta.pl/wi10776/86e44c51001f621f50a27480/
kwadratury_szeroko


Załóżmy, że naszym celem jest osiągnięcie dokładności względnej.
10^-12. Daleko od dokładności numerycznej, a jednocześnie ambitnie.


Metody rzędu 1 (midpoint, trapez) potrzebują ponad _miliona_
wyliczeń funkcji podcałkowej.

Kwadratury 3 rzędu (3/8, simpson, dwuponktowy gauss) wystarczy
około 2000 punktów! (punktów, nie przedziałów)

Jak dla mnie różnica między 1000000 a 2000 jest jednoznaczna.

Dla kolejnego oczka (5 rząd, dwie następne Newtona-coesta, Gauss
na 3 punktach, lobato na 4) potrzebują między 200 a 400 punktów.

Dwóm najlepszym (rząd 11 jeśli dobrze myśle) wystarcz poniżej
20 wywołań funkcji - trzy przedziały! (7*3 = 21, ale Lobatto
używa jednego wspolnego punktu dla sąsiadów, stąd ostatecznie 19)


pzdr
bartekltg

do góry

W dniu 2012-11-13 14:24, Roman W pisze:
> W dniu wtorek, 13 listopada 2012 13:23:00 UTC użytkownik AK napisał:
>
>> No i owszem :) Pokaz mni jedno miejsce gdzie ktos w tej dyskusji twierdzil
>> ze simson jest najlepszy (np lepszy niz splajny ?).
>> My tylko obalamy twoje idiotyzmy ze trappezy sa lepsze od simpsona !
>> A jak zszyjesz te twoje trapezy to co ?
>> Chyba robi cie sie prosta (no fakt ze pod nia mozna latwo "numerycznie"
>> calke policzyc :)
>
> http://www.johndcook.com/blog/2010/12/02/three-surpr
ises-with-the-trapezoid-rule/

2. Although the trapezoid rule is inefficient in general, it can be
shockingly efficient for periodic functions.
3. The trapezoid rule can also be shockingly efficient for analytic
functions that go to zero quickly, so called double exponential functions.

W tym nie ma nic szokującego dla kogokolwiek, kto widział
oszacowanie błędu dla tego typy kwadratur;)

Błąd pojedynczego trapezu to

int_a^b f(x) - I (f,a,b) = -f'' (c) *(b-a)^3 / 12


f'' (c) to wartość drugiej pochodnej funkcji
całkowanej w punkcie c \in [a,b].

Podkreślam, to jest wzór na błąd, nie jego szacowanie.
Nie wiadomo jedynie, ile wynosi c;)

Teraz budujemy kwadraturę złożoną. Wiele n przedziałów [a_{i-1},a_i]
a_0 = a, a_n =b długości h każdy.

łączny błąd wynosi

int_a^b f(x) - I (f,a,b) = -(h)^3 / 12 sum_{i=1}^{n} f'' (c_i)
= -(h)^3 / 12 h sum_{i=1}^{n} f'' (c_i)

Z tym, że c_i \in [a_{i-1},a_i]

Kolejne c_i są dowolne, ale ograniczone do swoich równo
oddalonych pudełek...


To czym jest h sum_{i=1}^{n} f'' (c_i)

Pewnym przybliżeniem całki z f'' na [a,b]!
A ta całka wynosi dla podanych przypadków 0.
Wartości błędów w poszczególnych przedziałach
wzajemnie się znoszę, z grubsza w takim tempie,
jak zbiega to oszacowanie na int f''.


BTW, są metody całkowania opierające się na prostej
kwadraturze i analizie zachowania kolejnych pochodnych
na brzegach. Nazwa mi wyleciała, natknąłem się chyba
przy okazji problemu z całkowanie funkcji szybko
oscylujących.

pzdr
bartekltg

do góry

W dniu 2012-11-13 14:22, AK pisze:

>> Nie, nie twierdzę że same parabole "są gorsze" - ale że z parabolami
>> masz problem jak je zszyć (i tu przychodzą z pomocą modne quintic
>> splines).
>
> No i owszem :) Pokaz mni jedno miejsce gdzie ktos w tej dyskusji twierdzil
> ze simson jest najlepszy (np lepszy niz splajny ?).

Nie są. Najlepsze 'uniwersalne' metody to adaptatywne Gauss-Kronrody;-)

> My tylko obalamy twoje idiotyzmy ze trappezy sa lepsze od simpsona !
> A jak zszyjesz te twoje trapezy to co ?

Hmm, czy mi się wydaje, czy właśnie wykryłeś kolejny genialny
pomysł sławka - zszywanie kwadratur wyższego rzędu, by miały
ciągłe pochodne;)

Na tym polegają kwadratury złożone. Całkujemy odcinkami
wielomianem wysokeigo stopnia i każdy odcinek jest obliczony
z odpowiednią precyzją. Nieciągłości pochodnych nie wpływają
negatywnie na wynik.

Można użyć splajnów. Ale całkowanie funkcji w najczęstszym
przypadku na tym niewiele korzysta (są przypadki, gdzie jest
to korzystne, ale też raczej bym szukał ich w okolicach
zbiór punktów -> aproksymacja minimalnokwadratowa splajnami (nie
interpolacja!), bo wiemy, że ukłąd jest np C^2 -> całkowanie).

Splajny kubiczne na oko będą zbiegać jak inne metody 3 rzędu.


pzdr
bartekltg

do góry

W dniu 2012-11-13 13:37, AK pisze:
> Użytkownik "slawek" <h...@s...pl> napisał:
>
>> Totalne niezrozumienie problemu: tobie nadal wydaje się, że możesz sam
>> sobie określać ile razy i w jakich "węzłach" wywołasz sobie funkcję
>> f(x). A tym razem problem był i jest taki, że masz z góry zadany ciąg
>> par (x,y), dla ułatwienia x[n] = n * h.

> BTW: W przypadku danych otrzynanych z pomiarow, (a wiec obardzonych
> bledem) w ogole nie stosuje sie tego typu interppolacji, ale aproksymacje
> i jesli juz "surowymi" wielomianami to przynajmniej poprzez jakis
> nawet najprymitywniejsze "wygladzanie" danych (chcby wielomianami Gramma
> z - co bardzo wazne - "automatycznyn"/statystycznym doborem stopnia).

To dodam jeszcze, że jeśli mamy nierówno oddalone węzły,
bo np przyszły one z pomiarów (obradzonych pomijalnym
błędem), a skądinąd wiemy, że funkcja jest odpowiedniej klasy
gładkości (np to laserowe pomiary pozycji satelity:)
nadal możemy użyć metody wyższego rzędu.

Maszynkę do wyliczania wag dla kwadratur na dowolnych węzłach
zamieściłem w kodach w wątku.

Oczywiście trzeba to zrobić ostrożnie albo z głową, bo odpowiednie
układy węzłów mogą dać gorsze wyniki niż chwilowe obniżenie węzłów.


Takie rozwiązanie dopasowywania wielomianów do 'dziwnych' węzłów
stosuje się przy niektórych wielokolorowych adaptatywnych
metodach całkowanie równań różniczkowych zwyczajnych.

BTW: można nie znaczy, że zawsze jest to dobry pomysł:)


> 2. bezkosztowe liczenie funkcji, ale za to "dosc szybkie". Tak szybkie
> ze nie nadazysz z "wyliczaniem" online "poprawek" np do korekcji
> przyslowiowego narzedzia skrawajacego dla twoich 10 000 punktow.
> Tu _tez_ kluczem jest jak najlepsza zbieznosc metody po to aby
> moc maksymalnie zmniejszyc ilosc "probke" przy zachowaniu
> zalozonej dokladnosci metody. Tu _tez_(co Bartek
> dobitnie na wykresach pokazal) twoje trapezy czy prostokaty
> sa w tyle.

Milion próbek trapeza = 2000 próbek simpsona = 18 próbek gaussa
dość niskiego rzędu;)


> PS: Nie twierdze ze kwadratury Newtona-Coatsa sa super.
> Nie sa. daleko im do tego.
> Chocby dlatego, ze sa nieciagle w wezlach.

Dla całkowania to nie przeszkadza (dla ściśle wyliczonych punktów,
jak są błędy, aproksymacja splajnami zaczyna coś dawać... ale wiadomo,
nie ma metod uniwersalnych)

> PS1: tak naprawde to wiekszosc poruszanych tu rzeczy to prostota i wrecz
> podstawy/abc wrecz elementarnej numeryki.

https://usosweb.mimuw.edu.pl/kontroler.php?_action=a
ctionx:katalog2/przedmioty/pokazPrzedmiot%28prz_kod:
1000-113aMOBa%29

To przedmiot dla wszystkich, nie dla zainteresowanych:-)

pzdr
bartekltg

do góry

mi wydaje sie troche sporo jesli jeden
simpson przybliza tyle ile 500 trapezów

a jakim to dokladnie kodem bylo calkowane?
(i co zrobilo też taki wykres?)

do góry

Użytkownik "bartekltg" <b...@g...com> napisał:

> zbiór punktów -> aproksymacja minimalnokwadratowa splajnami

Dokladnie o tym wzpomnialem.
No ale slawek aproksymacje czy wygladzanie uwaza za
cytuje ""psucie wynikow" :)

AK

do góry

W dniu 2012-11-13 20:09, kenobi pisze:
> mi wydaje sie troche sporo jesli jeden
> simpson przybliza tyle ile 500 trapezów

To nie jest liniowy przelicznik!

http://c.wrzuta.pl/wi10776/86e44c51001f621f50a27480/
kwadratury_szeroko

Im więcej punktów, tym korzyści większe.

Jeśli naszym celem jest 10^-6 trzeba dla tej funkcji
~60 punktów w Simpsonie i 2000 w trapezie. Więc tylko 30 razy lepiej.

Efekt bierze się z stąd, że nasze funkcje są porządne
i coraz bardziej lokalnie przypominają prostą lub
parabolę. Tylko parabola jest znacznie lepszym przybliżeniem
i funkcja znacznie wcześniej porządnie przypomina parabolę
(właściwie to wielomian 3 stopnia) niż prostą.


> a jakim to dokladnie kodem bylo calkowane?
> (i co zrobilo też taki wykres?)

Podawałem pełne kody w pierwszym poście.
https://skydrive.live.com/redir?resid=25C9E1E8909658
1!2201&authkey=!AGMwPjeQhRXoDe0


W tym też poście masz opis użytych metod.

Kody obejmują wszytko, wyznaczanie współczynników
kwadratur, właściwe obliczenia całki, wyrysowanie
tego.

Ostatni rysunek powstał przez wklepanie
testuj_calke(@(x)sin(x).*exp(-x),0,5,0.5022749400837
603657,20,10000000)
print -dpng -r300 szeroko

Po kolei parametry funkcji to:

całkowana funkcja, dwie liczby oznaczajace przedział całkowania,
wartość dokładna obliczona gdzie indziej symbolicznie,
dwie liczby oznaczające minimalną i maksymalną ilość punktów
braną do kolejnych testów.
Ostatnią wartość do zabaw należy zmniejszyć, bo ta wymaga długich
obliczeń, a była potrzebna, by w końcu metody pierwszego rzędu
osiągnęły maksymalną dokładność.

pzdr
bartekltg

do góry

W dniu wtorek, 13 listopada 2012 20:34:58 UTC+1 użytkownik bartekltg napisał:
> W dniu 2012-11-13 20:09, kenobi pisze:
>
> > mi wydaje sie troche sporo jesli jeden
>
> > simpson przybliza tyle ile 500 trapezów
>
>
>
> To nie jest liniowy przelicznik!
>
>
>
> http://c.wrzuta.pl/wi10776/86e44c51001f621f50a27480/
kwadratury_szeroko
>
>
>
> Im więcej punktów, tym korzyści większe.
>
>
>
> Jeśli naszym celem jest 10^-6 trzeba dla tej funkcji
>
> ~60 punktów w Simpsonie i 2000 w trapezie. Więc tylko 30 razy lepiej.
>
>
>
> Efekt bierze się z stąd, że nasze funkcje są porządne
>
> i coraz bardziej lokalnie przypominają prostą lub
>
> parabolę. Tylko parabola jest znacznie lepszym przybliżeniem
>
> i funkcja znacznie wcześniej porządnie przypomina parabolę
>
> (właściwie to wielomian 3 stopnia) niż prostą.
>
>
>
>
>
> > a jakim to dokladnie kodem bylo calkowane?
>
> > (i co zrobilo też taki wykres?)
>
>
>
> Podawałem pełne kody w pierwszym poście.
>
> https://skydrive.live.com/redir?resid=25C9E1E8909658
1!2201&authkey=!AGMwPjeQhRXoDe0
>
>
>
>
>
> W tym też poście masz opis użytych metod.
>
>
>
> Kody obejmują wszytko, wyznaczanie współczynników
>
> kwadratur, właściwe obliczenia całki, wyrysowanie
>
> tego.
>
>
>
> Ostatni rysunek powstał przez wklepanie
>
> testuj_calke(@(x)sin(x).*exp(-x),0,5,0.5022749400837
603657,20,10000000)
>
> print -dpng -r300 szeroko
>
>
>
> Po kolei parametry funkcji to:
>
>
>
> całkowana funkcja, dwie liczby oznaczajace przedział całkowania,
>
> wartość dokładna obliczona gdzie indziej symbolicznie,
>
> dwie liczby oznaczające minimalną i maksymalną ilość punktów
>
> braną do kolejnych testów.
>
> Ostatnią wartość do zabaw należy zmniejszyć, bo ta wymaga długich
>
> obliczeń, a była potrzebna, by w końcu metody pierwszego rzędu
>
> osiągnęły maksymalną dokładność.
>
nie znam tego języka .m :/

coz mozna powiedziec, same oczywiste sprawy,

ogolnie jestem zwolennikiem trapezow, ale jesli ten
simpson daje taką dobrą celnosc i dla pewnej klasy
funkcji malym kosztem polepsza dokladnosc wyniku
az tak bardzo, no to warto go robic. [Zeby sie
osobiscie przekonac musialbym osobiscie klepnac
testy w c ale jako ze ja sie tym za specjalnie nie
zajmuje jestem raczej sklonny po prostu uwierzyc
temu przykladowi, pewnie moze tak byc (choc zdziwilo
mnie ze celnosc tego wielomianu jest az tak duza
jak niecelnosc lamanej aż na kilkuset punktach)

(fir)

do góry

W dniu 2012-11-13 22:51, kenobi pisze:

>>
> nie znam tego języka .m :/

Matlab. Prosty jest. Darmowe wersje to octave
i od biedy scilab.

Czy ten kod pójdzie w octave, nie wiem.
Pewne funkcje mogą mieć inne nazwy.
Ale na oko jest dobrze:
http://www.chemie.fu-berlin.de/chemnet/use/info/octa
ve/octave_10.html


> coz mozna powiedziec, same oczywiste sprawy,

Bo gadamy tu o oczywistych sprawach.
Kwadratura to policzenie funkcji w odpowiednich
miejscach, pomnożenie uzyskanych wartości
przez odpowiednie liczby i zsumowanie.

Wyznaczenie pewnych układów parametrów, które
działąją lepiej niż wersje oczywiste wymaga
kilku operacji. Prostych, ale trzeba wiedzieć,
jakich;) Np z wikipedii;)

Ok, trochę dłubaniny technicznej, żeby to się
zgrabnie działo.


Ale zawsze możesz sprawdzić, ile to wymaga.
mi to zajęło wieczorek:)

pzdr
bartekltg

do góry

W dniu wtorek, 13 listopada 2012 23:10:48 UTC+1 użytkownik bartekltg napisał:
> W dniu 2012-11-13 22:51, kenobi pisze:
>
>
>
> >>
>
> > nie znam tego języka .m :/
>
>
>
> Matlab. Prosty jest. Darmowe wersje to octave
>
> i od biedy scilab.
>
>
>
> Czy ten kod pójdzie w octave, nie wiem.
>
> Pewne funkcje mogą mieć inne nazwy.
>
> Ale na oko jest dobrze:
>
> http://www.chemie.fu-berlin.de/chemnet/use/info/octa
ve/octave_10.html
>
>
>
>
>
> > coz mozna powiedziec, same oczywiste sprawy,
>
>
>
> Bo gadamy tu o oczywistych sprawach.
>
> Kwadratura to policzenie funkcji w odpowiednich
>
> miejscach, pomnożenie uzyskanych wartości
>
> przez odpowiednie liczby i zsumowanie.
>
>
>
> Wyznaczenie pewnych układów parametrów, które
>
> działąją lepiej niż wersje oczywiste wymaga
>
> kilku operacji. Prostych, ale trzeba wiedzieć,
>
> jakich;) Np z wikipedii;)
>
>
>
> Ok, trochę dłubaniny technicznej, żeby to się
>
> zgrabnie działo.
>
>
>
>
>
> Ale zawsze możesz sprawdzić, ile to wymaga.
>
> mi to zajęło wieczorek:)
>
>

nie moja dzialka :O, mam co innego do roboty :o
(ostatnio troche czytam o shaderach - i moge powiedziec ze troche beznadziejne oneż
są)

do góry