Metoda trapezów

Przeszukałem “cały” internet wpisując hasło całkowanie metodą trapezów, a potem w języku angielskim: trapezoidal rule in numerical integration. Nie znalazłem jednej podstawowej interesującej mnie informacji, którą mam czarno na białym, w jednym z moich “zasobów” papierowych. Możliwe, że to jest tak podstawowe i oczywiste, że nigdzie indziej tego nie znalazłem.

We wszystkich źródłach i gotowych algorytmach, to użytkownik programu ma podać przedział i ilość trapezów, dzielących ten przedział. A co jak biedny użytkownik [i programista] nie ma o tym zielonego pojęcia [albo nawet bladego]. Jest na to prosta rada. Użytkownik podaje 1 przedział, i program błyskawicznie “wypluwa” odpowiedź. Super. Zadowolony uzer [użytkownik] szybko zapisuje wynik na kartce i chytrze podaje dwa razy większą wartość - ilość trapezów. Znowu to samo, Uzer spisuje wynik i porównuje z poprzednim. Duża różnica. Może więc od razu dam podział na 1000000 trapezów, myśli sobie chytrze. Komputer teraz dostał zadyszki, a uzer w między czasie policzył wynik ręcznie lub na wolframiealpha i gdy w między czasie program skończył, okazało się, że wynik programu jest zadawalająco dokładny.

Ale , żeby nie zabierać Ci już więcej twojego cennego czasu drogi czytelniku, przechodzę do mojego problemu i do konkretów.

Rozwiązaniem jest oczywiście, aby to program porównywał kolejne przybliżenia i zdecydował, kiedy skończyć. Użytkownik, może ewentualnie zadecydować jaka dokładność go interesuje lub ile czasu jest skłonny czekać, na może mniej dokładny, ale jakiś nawet przybliżony wynik. Może też nic nie podawać oprócz przedziału, a program sam zadecyduje o max dokładności.

Może to jest tak oczywiste i było na pierwszej lekcji wstępu do metod numerycznych, że nigdzie w internecie nie znalazłem tak “skonstruowanego” programu ani nawet o tym wzmianki.

METODA TRAPEZÓW

Nie będę opisywał w całości i od początku, bo można znaleźć opisy i filmy w necie, np: http://www.algorytm.edu.pl/algorytmy-maturalne/metoda-trapezow.html3

Końcowy wzór z tej strony, pomijając wartości bezwzględne można przekształcić do postaci:

Pn ~= h*(1/2 (f(a) +f(b)) + f(x1) +..+ f(xn))

Wyobraźmy sobie przedział a=0, b = 16; i 15 punktów x w tym przedziale oznaczonych odpowiednio x1, x2, x3, x4, x5, x6 … x15

Pierwsze przybliżenie - 1 trapez h0 = 16
P0 ~= h0 * (0.5 (f(a) + f(b))
Drugie przybliżenie - 2 trapezy h1 = 8
P1 ~= h1 *(0.5(f(a) +f(b)) + f(x8))
Trzecie przybliżenie - 4 trapezy h2 = 4
P2 ~= h2 *(0.5(f(a) +f(b)) + f(x8) + f(x4) + f(x12))
Czwarte przybliżenie - 8 trapezów h3 = 2
P3 ~= h3 *(0.5(f(a) +f(b)) + f(x8) + f(x4) + f(x12) + f(x2) + f(x6) + f(x10) + f(x14))
Piąte przybliżenie - 16 trapezów h4 = 1
P4 ~= h4 *(0.5(f(a) +f(b)) + f(x8) + f(x4) + f(x12) + f(x2) + f(x6) + f(x10) + f(x14) + f(x1) + f(x3 + f(x5) +f(x7…x15))

I dalej jeżeli dzielimy przez dwa - podwajamy ilość trapezów możemy robić podobnie. - jeżeli już to widzisz.
Co nam to może dać w zadaniu Metody Trapezów na Spoju.
Możemy szybko policzyć górną granicę późniejszego przeszukiwania binarnego:
Dodatkowo, należy uwzględnić dokładność, z jaką możemy lub chcemy liczyć na liczbach typu double, np EPS = 0.00000001.

while (S*h + EPS < Ppodane)
   n *= 2
   h /= 2
   S = S_old + S_poprawka
}
 jeżeli S*h + - EPS == Ppodane) to mamy rozwiązanie, 
w przeciwnym wypadku dolną granicą jest  n/2, a górną n

Tyle z gróbsza, jeżeli jesteś zainteresowany lub masz pytania, śmiało pytaj.

Jeszcze dwie sprawy.
W zależności od znaków [położenie] wykresu funkcji, możemy mieć sytuację, że przedział granic jest 1 … n lub n … 1

Albo inaczej. Jeden przypadek, kiedy kolejne przybliżenia są z coraz mniejszym nadmiarem, a drugi przypadek, inne ułożenie wykresu, kiedy kolejne przybliżenia są z coraz mniejszym niedomiarem. Czyli:

while (S*h - EPS > Ppodane)

...```