Nawiązując do wcześniejszego wątku nt. całkowanie numerycznego (dane są jako
tablica y[1],...,y[n], całka jest od a do b, przy czym x = (j - 1) * h,
gdzie j jest indeksem do tablicy y-ków; ważne że n jest w miarę duże -
kilka-kilkanaście tysięcy, oraz że wartości funkcji nie można obliczyć gdzie
się chce - tzn. są zadane w tablicy i nie da się mieć jakiś dodatkowych
wartości.)

1. Tu macie wykresy porównujące trzy różne algorytmy, jakie zastosowałem -
http://drop.io/vokdmhp - procedurki simple, tpz, zint.

Simple to po prostu łopatologia stosowana - jeżeli punkt x jest wewnątrz
przedziału (a,b), to wchodzi do całki jako f(x)*h. Czyli po prostu pętla po
j ze sprawdzaniem czy (j-1)*h jest większe od a lecz mniejsze niż b. Jeżeli
równe - to jeszcze waga 1/2 . Ku mojemu zdumieniu, całkiem dokładna metoda -
w zasadzie asymptotycznie powinna być, ale okazało się, że nawet dla nie za
dużej liczby punktów też działa prawie tak dobrze jak bardziej wyrafinowane.
To NIE JEST metoda trapezów, bo a,b są rzeczywiste (float point), i
niekoniecznie przypadają tak, że a/h ewentualnie b/h jest całkowite.

Tpz to po prostu trapeziki.

Zint to całkowanie przez naturalne spliny kubiczne plus końcówki
aproksymowane wielomianem stopnia trzeciego (ale już nie zszywanym ze
splinem).

Teraz ciekawostka (całkowana była funkcja sin(x) w przedziale od 0 do pi,
powinno "chyba" wyjść 2.000000000000000000000000000000000000 jako wynik :)
: zarówno zint jak i tpz dają niemal taki sam błąd w obliczeniach, co jest
ździebko szokujące (patrz test2b.png) . Tpz daje niedomiar (funkcja jest
wypukła, tego się można było spodziewać), zint daje nadmiar (a to ciekawe,
znaczy się spline nieco się "wybrzusza" bardziej niż sama funkcja). Jak brać
logartym z modulu błędu względnego - wychodzi - surprise, surprise - że
10-krotne zwiększenie kroku daje 100 krotne zwiększenie dokładności - czyli
"doświadczalnie" mamy błąd rzędu O(h^2).

Aż się prosi, aby brać średnią arytmetyczną zint i tpz - patrz rysunek
test1a.png.


2. Jak to draństwo można ulepszyć? A. Zamiast naturalnych funkcji sklejanych
brać takie, które dla wielomianu stopnia 3-ciego będą go wiernie odtwarzać.
B. zszyć wielomian do liczenia "końcówek" z całym długim splinem. C. Zamiast
"cubic spline" (App. Num. Analysis., Curtis F. Gerald et al.) wziąć
monotoniczne funkcje sklejane Hermite'a; D. zastosować podejście zbliżone do
całkowania według Romberga (jeżeli wynik to I+epsilon*h^alfa, to przeliczyć
dla trzech różnych h i ekstrapolować metodą Aitkena do h = 0).

Moim zdaniem A/B niewiele zmienią (tylko końce, a to można pominąć gdy
(b-a)/h jest dostatecznie duże). Natomiast C może być ciekawe. D. jest chyba
najprostsze do zastosowania. Powstaje też ciekawe pytanie jak to wszystko
rozszerzyć do całek podwójnych (nie przez iterację rzecz jasna, bo to
trywialne).

3. Znalazłem małego wrednego buga w zint (zgubił się współczynnik, testy
oczywiście na wersji poprawionej). Dość dużo trzeba jeszcze wymęczyć w
zint - np. obsługę przypadków gdy są tylko 3 punkty, albo są 4 punkty które
leżą poza przedziałem (a,b). Takie tam.

4. Dla innych funkcji tabelka poniżej. Wyniki oznaczone jako "mathematica"
są wyliczone z analitycznych wartości całek w arytmetyce arbitralnej
precyzji. Krok h celowo "nieokrągły".

5. Dalej aktualne jest pytanie o najlepszy sposób całkowania danych
doświadczalnych itp. (jest mnoooostwo zastosowań dla tego rodzaju
algorytmu).

slawek

nmax = 1048576
step = 0.0000030386438397372824

Itegrate[Sin[x],{x,0,Pi}]
simple : 0.19999999999990321076D+01 -0.4839D-12
trapez : 0.19999999999984523491D+01 -0.7738D-12
zint : 0.20000000000015352164D+01 0.7676D-12
mathematica : 0.20000000000000000000D+01

Itegrate[Sin[x],{x,0,Pi/2}]
simple : 0.99999870039620486484D+00 -0.1300D-05
trapez : 0.10000000437106320028D+01 0.4371D-07
zint : 0.10000000437121594477D+01 0.4371D-07
mathematica : 0.10000000000000000000D+01

Itegrate[Sin[100*x]/x,{x,0,Pi}]
simple : 0.15676132924346011244D+01 0.4028D-11
trapez : 0.15676132924528567436D+01 0.1567D-10
zint : 0.15676132924034666960D+01 -0.1583D-10
mathematica : 0.15676132924282872860D+01

Itegrate[1000/Pi+Sin[100*x]/x,{x,0,Pi}]
simple : 0.10015672417271147197D+04 -0.3710D-06
trapez : 0.10015676132924576223D+04 0.2929D-13
zint : 0.10015676132924033936D+04 -0.2486D-13
mathematica : 0.10015676132924282911D+04

Itegrate[Log[x+1]/x,{x,0,1}]
simple : 0.38629503701836448437D+00 0.1750D-05
trapez : 0.38629436111950254951D+00 -0.1005D-11
zint : 0.38629436112027532024D+00 0.9959D-12
mathematica : 0.38629436111989062796D+00

do góry

slawek <s...@h...pl> napisał(a):

> Teraz ciekawostka (całkowana była funkcja sin(x) w przedziale od 0 do pi,
> powinno "chyba" wyjść 2.000000000000000000000000000000000000 jako wynik :)
> : zarówno zint jak i tpz dają niemal taki sam błąd w obliczeniach, co jest
> ździebko szokujące

Nie jest szokujące. Błąd na części rosnącej jest równy błędowi na części
malejącej z przeciwnym znakiem. Metoda prostokątów jest obarczona
większym błędem, tyle że w tym konkretnym przypadku błędy wzajemnie
się znoszą.


> 2. Jak to draństwo można ulepszyć?
W kilku eksperymentach które przeprowadziłem (mogę dać kod w C++)
aproksymacja wielomianem 5-go stopnia dawała ten sam błąd w czasie
do 10 razy krótszym niż 3-go stpnia.

Pozdrawiam

--
Wysłano z serwisu Usenet w portalu Gazeta.pl -> http://www.gazeta.pl/usenet/

do góry

On Sat, 25 Jul 2009 19:08:07 +0000 (UTC), "Mariusz Marszałkowski"
<b...@N...gazeta.pl> wrote:

>slawek <s...@h...pl> napisał(a):
>
>> Teraz ciekawostka (całkowana była funkcja sin(x) w przedziale od 0 do pi,
>> powinno "chyba" wyjść 2.000000000000000000000000000000000000 jako wynik :)


In your dreams

A.L.

do góry

Użytkownik "Mariusz Marszałkowski" <b...@N...gazeta.pl> napisał w
wiadomości grup dyskusyjnych:h4fl6n$fm4$...@i...gazeta.pl...
> Nie jest szokujące. Błąd na części rosnącej jest równy błędowi na części
> malejącej z przeciwnym znakiem. Metoda prostokątów jest obarczona
> większym błędem, tyle że w tym konkretnym przypadku błędy wzajemnie
> się znoszą.

Jest szokujące - ponieważ *znaki* błędów są przeciwne.

Co pozwala zrobić przez proste uśrednienie znacznie lepszą metodę, tu masz
obrazek:http://drop.io/aqtdw2y

> W kilku eksperymentach które przeprowadziłem (mogę dać kod w C++)
> aproksymacja wielomianem 5-go stopnia dawała ten sam błąd w czasie
> do 10 razy krótszym niż 3-go stpnia.

a. Dla ilu punktów?

b. Funkcje takie jak sinus czy logarytm ładnie się rozwijają w szereg
Taylora (kryterium Abela jest spełnione od pierwszego wyrazu). W docelowym
zastosowaniu tak nie będzie - a co najgorsze - wzrost stopnia wielomianu źle
zadziała gdy y[n] >> y[k] dla k = 1,2,...,n-1

slawek

do góry

Użytkownik "A.L." <a...@a...com> napisał w wiadomości grup
dyskusyjnych:73nm6558hfckfanddpjvev20i8d5nojhh6@4ax.
com...
> On Sat, 25 Jul 2009 19:08:07 +0000 (UTC), "Mariusz Marszałkowski"
> <b...@N...gazeta.pl> wrote:
>
>>slawek <s...@h...pl> napisał(a):
>>
>>> Teraz ciekawostka (całkowana była funkcja sin(x) w przedziale od 0 do
>>> pi,
>>> powinno "chyba" wyjść 2.000000000000000000000000000000000000 jako wynik
>>> :)
>
>
> In your dreams
>
> A.L.

Zwykle tak nie piszę, ale dla A.L. zrobię wyjątek.

Primo, policz sobie analitycznie całkę, łosiu.

Secundo, jak tego jeszcze nie potrafisz, to tak z pamięci, sinus ma
pierwotną cosinus plus stała, całka od 0 do pi to dwie całki od zero do 90
stopni, cosinus z zera jest 1, cosinus z 90 stopni to 0, czyli wychodzi
okrągłe 2. Jak widzisz, nawet jeszcze przed śniadaniem i ziewając potrafię
taką pierdołę policzyć. A jakby mi się nie chciało liczyć to po prostu
przypomniałbym sobie że amplituda jest sqrt(2) większa niż napięcie
skuteczne, a moc to kwadrat napięcia przez opór. Czy jakoś tak. Powinno być
2. Tak, nawet jak śpię, to wiem ile mi się powinno przyśnić.

Tertio, cudzysłowy widział? Smiley'a widział? Potrafił przeliterować
"arbitralna precyzja"? No to co trzaska dziobem (o klawiaturę)?

slawek

do góry

Użytkownik "slawek" <s...@h...pl> napisał w wiadomości grup
dyskusyjnych:h4gu5f$h3e$...@a...news.neostrada.pl
...
> b. Funkcje takie jak sinus czy logarytm ładnie się rozwijają w szereg
> Taylora (kryterium Abela jest spełnione od pierwszego wyrazu). W docelowym
> zastosowaniu tak nie będzie - a co najgorsze - wzrost stopnia wielomianu
> źle zadziała gdy y[n] >> y[k] dla k = 1,2,...,n-1

Dwie autopoprawki: dla małych wartości argumentów (jeżeli x = 1000 to
oczywiście nie, ale jeżeli Abs[x] < 1 to tak); w konkretnym przypadku n jest
rzędu 10 tysięcy. Już dla 20 tysięcy nie chce mi się czekać aż obliczenia
generujące y-ki się zakończą (algorytm z kosztem O(n^3), trwa to parę
godzin.

slawek

do góry

slawek <s...@h...pl> napisał(a):

> Jest szokujące - ponieważ *znaki* błędów są przeciwne.
Szczerze to też doznałem szoku. Ale gdy zrozumiałem że tyle samo traci
"pod górkę" co zyskuje "z górki" to szok minął. :)

> a. Dla ilu punktów?
Różne funkcje, różne ilości punktów, wszystkich wyników nie pamiętam.
Było różnie, ogólnie było lepiej, a czasami ta sama dokładność nawet
w 10 razy krótszym czasie. Można było zredukować ilość punktów i dlatego
krótszy czas.

> zastosowaniu tak nie będzie - a co najgorsze - wzrost stopnia wielomianu
> źle zadziała gdy y[n] >> y[k] dla k = 1,2,...,n-1
Nie wiem, dlaczego wtedy jest źle?



--
Wysłano z serwisu Usenet w portalu Gazeta.pl -> http://www.gazeta.pl/usenet/

do góry

Użytkownik "Mariusz Marszałkowski" <b...@N...gazeta.pl> napisał w
wiadomości grup dyskusyjnych:h4h98b$q6j$...@i...gazeta.pl...
> Szczerze to też doznałem szoku. Ale gdy zrozumiałem że tyle samo traci
> "pod górkę" co zyskuje "z górki" to szok minął. :)

Nie, to nie tak. Po prostu algorytm z trapezami (I) i algorytm ze splinami
(II) dają dla danego przedziału (a,b) liczby różne od wyniku analitycznie
otrzymanego U o odpowiednio -epsilon i +epsilon, gdzie epsilon jest funkcją
kroku h. Oczywiście, że dla (0,pi/2) jest trochę np. nadmiar a dla (pi/2,pi)
jest niedomiar w danej metodzie, więc się znosi - to też jest, ale nie o to
chodzi. Chodzi o to, że dwa zupełnie różne algorytmy dają wartość odchylenia
od wyniku takie same co do wartości bezwzględej (z różnicą do około 1%) dla
różnych kroków itd. - tylko znaki są różne. Czyli średnia z obu metod jest
dużo dokładniejsza i błąd jej malej znacznie szybciej niż O(h^2). Patrz
http://drop.io/aqtdw2y - dla małych h, czyli dużych n wynik jest 10 tysięcy
razy dokładniejszy! Błąd jest rzędu O(h^3) - wykres to ładnie prezentuje.

> Nie wiem, dlaczego wtedy jest źle?

Niekoniecznie źle być musi. Po prostu "cały wielki program" liczy y[j+1]
jako funkcję poprzednich y[1],y[2],...,y[j] - całki i takie tam. Możliwe
jest, że y[j] leci do nieskończoności. Biorąc dużo punktów opóźnia się
moment wykrycia takiej rozbieżności - procedura wygładza zachowanie funkcji,
która w pobliżu osobliwości mało przypomina grzeczny wielomian. To trochę
bardziej skomplikowane - bo to i układ równań jest, i jeszcze y jest
zespolone i parę jeszcze zagwozdek po drodze.

slawek

do góry

On 25 Lip, 14:26, "slawek" <s...@h...pl> wrote:

> 5. Dalej aktualne jest pytanie o najlepszy sposób całkowania danych
> doświadczalnych itp. (jest mnoooostwo zastosowań dla tego rodzaju
> algorytmu).
>
> slawek
Dalej aktualne jest pytanie o to co rozumiesz pod pojęciem
najlepszej metody całkowania? Przedstawione wyniki, nie bronią
Twojej tezy, że metoda oparta o splajny jest dużo lepsza od trapezów.

Z poważaniem,
Wit Jakuczun

do góry

Wit Jakuczun <w...@g...com> napisał(a):

> Dalej aktualne jest pytanie o to co rozumiesz pod poj=EAciem
> najlepszej metody ca=B3kowania? Przedstawione wyniki, nie broni=B1
> Twojej tezy, =BFe metoda oparta o splajny jest du=BFo lepsza od trapez=F3w.

Jakie w ogóle stosuje się kryteria oceny metod całkowania?

Pozdrawiam




--
Wysłano z serwisu Usenet w portalu Gazeta.pl -> http://www.gazeta.pl/usenet/

do góry