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

do góry

On 2012-11-12, slawek <h...@s...pl> wrote:
> Forumowe trolle-specjaliści na wszystkim się znają... tylko akurat na
> pl.comp.programming nie odróżniają Fortranu od C, nie potrafią policzyć
> samodzielnie numerycznie całki, itd. itp.

Akurat nie brakuje speców od numeryki, którym wszyscy na tej grupie
moglibyśmy najwyżej temperować ołówki, a którzy mogliby C od Fortrana
nie odróżnić, więc w czym problem?

pozdrawiam,
PK

do góry

W dniu 2012-11-11 18:44, AK pisze:
> PS2: Masz racje. Przedzial i funkcje tez sobie dobral "pod przyklad"

Nie. Schrzanił implementację.

Jego przykład:
http://c.wrzuta.pl/wi4888/9c4c1562000019e750a147f7/s
in_exp

Trapez to zielona *, simpson to czerwona *.
Wszytko zgodnie z teorią.

pzdr
bartekltg

do góry

W dniu 2012-11-12 21:48, bartekltg pisze:
>
> Wątkowy przykład, sin(x) exp(-x) na [0,5]
> http://w18.wrzuta.pl/obraz/70M5xduXyOU/sin_exp


Na w szelki wypadek, bo dokopanie do rozsądnej
rozdzielczości zrobiło się nieintuicyjne :

http://c.wrzuta.pl/wi4888/9c4c1562000019e750a147f7/s
in_exp


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


http://c.wrzuta.pl/wi19003/0b599d75000065a450a147f6/
sin10_exp

pzdr
bartekltg

do góry

Przedstawie moze dla jasnosci modelową 'sytuację
z dresem', modelowa sytuacja dresiarska jest calkiem prosta i wyglada mw tak:

1. przychodzi ktos i mowi cos na temat, czesto z sensem, czasem z mniejszy sensem a
czasem zadaje tylko pytania, ale trzyma sie tamatu i generalnie dyskusja jest
przedmiotowa i sklanialaby do dyskusji utrzymanej w przedmiotowym nurcie

2. teraz pojawia sie dres,

dres rzuca cos w stylu 'ty jestes głupi czy tylko udajesz' ew cos jeszcze bardzij
durnego,
uzywa tez oczywiscie slowa wytrycha pt 'troll';
dres nie jest zainteresowany normalnymi zagadnieniami, ew sygnalizuje tylko ze
jestes glupi i to co mowisz jest glupie
nie jest jednak w stanie sie do tego odniesc

Dresy produkują odzielny nurt na grupie, obok
tego prawdziwego i merytorycznego produkują
swoisty nurt jełopienia wypełniony jełopieniem
i dowalactwem

do góry

>
> Dresy produkują odzielny nurt na grupie, obok
> tego prawdziwego i merytorycznego produkują
> swoisty nurt jełopienia wypełniony jełopieniem
> i dowalactwem

Moze pozniej jeszcze dopracuje wyraźniejsze
kryterium (mnie cale te nurty dresiarskie
glownie meczą) w każdym razie dresa
łatwo jest poznac po własnie zmianie nurtu
z normalnego merytorycznego na dowałki,
czyli poznac mozna dosyc prosto

Nie wszystko oczywiscie co nie jest scisle programistyczne nie jest merytoryczne, bo
jakies dygresje nt wlasnych doswiadczen
itp tez są merytoryczne, jednak to co jest
ufundowane na dowale jest jełopieniem i dresiarstwem w 100% (mz wynika tio z tego
ze czesc ludzi jest po prostu dresiarzami
i lubią przyjełopić, na durnote ciezko cos poradzic)

do góry

Użytkownik "kenobi" <p...@g...com> napisał:

> Moze pozniej jeszcze dopracuje wyraźniejsze kryterium

Spojrz w lustro. Tam masz prawzor dresa.

AK

do góry

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

>W dniu 2012-11-11 18:44, AK pisze:
>> PS2: Masz racje. Przedzial i funkcje tez sobie dobral "pod przyklad"
>
> Nie. Schrzanił implementację.

Faktycznie (nie do wiary).
Nie zauwazylem. Prawde mowiac nie wglebialem sie. Szkoda mi bylo czasu.
Tym wieksze uznanie dla Ciebie za tak porzadne podejscie do tematu.
Zreszta bardzo przydatne i uczace.
Szczery szacunek i podziekowanie.

PS: Wlasciwie po Twym poscie powinno z wiadomej strony nastapic milczenie.
Nie liczylby jednak na to. Te duet (slawek, kenobi) jest wyjatkowo odporny na
fakty :)

AK

do góry

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?

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

Niestety, pomijając dość nonszalancki kod źródłowy (np. "karma_dla_plota"
jako identyfikator, używanie iloczynu Kroneckera choć można prościej, brak
komentarzy, niezbyt przemyślane to i owo... oczywiście /trochę/ się
czepiam), takie podejście ma jedną istotną wadę/zaletę - ukrywa samo
sumowanie jako liczenie długości wektora przez Matlab/Octave. Czyli
zasadnicza część algorytmu jest poza kontrolą - choć może to być i zaleta
("sumowanie słupka").

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

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. Mając takie - i tylko takie - dane masz obliczyć
możliwie dokładnie całkę - np. po to, aby mieć wartość skuteczną (tj.
średnią całkową). Jawna postać f(x) jest tylko do wygenerowania danych i
ewentualnie analitycznego oszacowania "prawdziwej" wartości - algorytm
całkujący nie ma i nie może mieć do niej dostępu.

Ok, jest z tym pewien /ekstra/ problem - wartość całki liczonej analitycznie
nie jest aż tak dobrym oszacowaniem - przecież jest to nic innego, niż całka
z funkcji interpolującej dane, a ta interpolacja może nie odpowiadać
rzeczywistości (fakt, że dla każdego n zachodzi y[n] = f(x[n]) niczego nie
przesądza). W istocie rzeczy "naprawdę prawdziwa" wartość całki jest
niemożliwa do obliczenia - bo dyskretne dane są określone na zbiorze o
mierze zero i takie tam. Ale to już jest trochę "przefilozofowane".

>Wartość dokładna wyliczona przez wolfrem alpha;)

Nie "wolfrem alpha", ale ale Wolfram Alpha - że złośliwie skomentuję. FYI,
ja preferuję normalną Mathematicę, bo Alfa się komercjalizuje (nic dziwnego)
i przestaje być tak wygodna jak była - silnik do całkowania mają zapewne ten
sam. Sympatyczne w Alpha było/jest to, że można wygenerować szczegółowy opis
rozwiązania.

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

A nie da się ich liczyć analitycznie, tj. na liczbach int? Mathematica nie
ma z tym trudności.

do góry

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?

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.

Kod ponizej (C++):

#include <cmath>
#include <iostream>
#include <iomanip>

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

template <class F>
double integrate_simpson(F f, double x0, double x1, size_t k)
{
const size_t n = 2 * k;
const double h = (x1 - x0) / n;
double sum = f(x0) + f(x1);
for (size_t j = 1; j < k; ++j) {
sum += 2 * f(x0 + 2 * j * h);
}
for (size_t j = 1; j <= k; ++j) {
sum += 4 * f(x0 + (2 * j - 1) * h);
}
return h * sum / 3.0;
}

int main()
{
std::cout << std::setprecision(16) << integrate_simpson(function, 0.0, 5.0, 10000/2)
<< std::endl;
}


RW

do góry