Całkowanie numeryczne metodą trapezów

Wstęp

Całkowanie metodą trapezów stanowi udoskonalenie metody prostokątów, w której wykorzystywana jest jedynie pojedyncza wartość ze środka przedziału. W metodzie trapezów pole pod wykresem przybliżane jest trapezami prostokątnymi o wysokości h każdy oraz podstawach długości f(xi), f(xi+h) dla i-tego prostokąta:

Wizualizacja całkowania numerycznego. Źródło: wikipedia

Tym samym całka przybliżana jest jako suma pól N trapezów w przypadku podziału przedziału całkowania na N podobszarów:

\int_{a}^{b}f(x)dx\approx \sum_{i=0}^{N-1}\frac{h \cdot{}(f(x_{i}) + f(x_{i+1}))}{2}

Ponadto zakładamy:

h = \frac{b-a}{N}\hspace{40mm}x_{i} = a + hi \hspace{40mm}b > a

Ponieważ obliczenie całki wprost z tego wzoru wymagałoby (pomijając obliczanie xi) 2N - 1 dodawań, N mnożeń oraz N dzieleń, możemy go przekształcić do poniższej, równoważnej postaci, wymagającej już tylko N dodawań, 1 mnożenia i 1 dzielenia.

\int_{a}^{b}f(x)dx\approx h \cdot{} (\frac{f(x_{0})+f(x_{N})}{2} + \sum_{i=1}^{N-1}f(x_{i}))

Implementacja w Pythonie

Dla przejrzystości pominiemy walidację argumentów i obsługę błędów.

Postać oryginalna

def trapeze_integration_simple(my_func, a, b, n):
  # Wysokość trapezu
  delta_x = (b-a)/n
  # Sumujemy pola trapezów, x_i liczymy każdorazowo od nowa
  # by uniknąć kumulowania się błędów numerycznych - kosztem dużej ilości mnożeń
  total = 0  
  for i in range(n):
    x = a + delta_x * i
    total += delta_x * (my_func(x) + my_func(x + delta_x)) / 2
  return total

# Wywołujemy naszą funkcję całkującą, przekazując wzór całkowanej funkcji jako lambdę
integral = trapeze_integration_simple(lambda x: x**2, 0.0, 1.0, 100)

Postać przekształcona

def trapeze_integration_modified(my_func, a, b, n):
  # Wysokość trapezu
  delta_x = (b-a)/n
  subtotal = my_func(a) + my_func(b)
  # Tym razem wyeliminujemy mnożenia, co krok dodając do x nasze delta_x
  # Musimy pamiętać, że błędy zaokrągleń mogą się skumulować i istotnie
  # zniekształcić wynik
  x = a
  for _ in range(1, n):
    subtotal += my_func(x)
    x += delta_x
  # Na końcu wykonujemy mnożenie i dzielenie
  return subtotal * delta_x / 2

# Wywołujemy naszą funkcję całkującą, przekazując wzór całkowanej funkcji jako lambdę
integral = trapeze_integration_modified(lambda x: x**2, 0.0, 1.0, 100)

Przykłady

Przykład 1: funkcja kwadratowa

TODO: przedstawić przykład krok po kroku z wynikami

Przykład 2: funkcja cosinus

TODO: przedstawić przykład krok po kroku z wynikami

0 komentarzy