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