Kwadratura wyjaśniona: metody numeryczne do przybliżania całek

Istnieje klasa całek oznaczonych, które po prostu nie poddają się rachunkowi całkowego na kartce papieru. Nie dlatego, że robisz coś źle. Niektóre funkcje, jak e^(−x²) czy sin(x)/x, w ogóle nie mają pierwiastków w postaci elementarnych funkcji. Gdy tak się dzieje, potrzebujesz innego podejścia: oblicz wartość funkcji w kilku starannie dobranych punktach i połącz te wartości, aby oszacować pole pod wykresem. To właśnie kwadratura.

Ten artykuł wyjaśnia, czym jest kwadratura, jak działa ogólne ramy postępowania oraz najważniejsze metody — od intuicyjnej reguły punktu środkowego po zaskakująco skuteczne podejście Gaussa. Omówimy też metody złożone, źródła błędów, strategie adaptacyjne oraz to, jak dobrać właściwą metodę do konkretnego problemu.

Czym jest kwadratura?

Słowo „kwadratura” jest stare. Wywodzi się z klasycznego problemu geometrycznego polegającego na skonstruowaniu kwadratu o tym samym polu co dany kształt — innymi słowy: „znaleźć pole”. Dziś termin ten oznacza numeryczne przybliżanie całki oznaczonej.

Sednem jest zamiana: zamiast szukać dokładnej funkcji pierwotnej, próbkujesz funkcję f(x) w zbiorze punktów zwanych węzłami, mnożysz każdą wartość przez przypisaną wagę i sumujesz wyniki:

Pojęcia kwadratura i całkowanie numeryczne znaczą zasadniczo to samo. „Kwadratura” częściej pojawia się w kontekście matematycznym; „całkowanie numeryczne” to codzienny termin w obliczeniach naukowych.

Jak działa kwadratura numeryczna

Każda reguła kwadraturowa działa na przedziale [a, b]. W tym przedziale wybiera zbiór węzłów xᵢ i przypisuje każdemu wagę wᵢ. Przybliżenie to ważona suma z powyższego wzoru. Między metodami różni się to, gdzie trafiają węzły i jak wyznacza się wagi. Tyle. Metoda trapezów, reguła Simpsona i kwadratura Gaussa to wariacje na ten sam schemat.

Błąd aproksymacji to różnica między ważoną sumą a prawdziwą wartością całki. Każda metoda wprowadza pewien błąd, a to, jak szybko on maleje, gdy dodajesz więcej węzłów (lub zmniejszasz szerokości przedziałów), odróżnia metodę dobrą od znakomitej. Pamiętaj o tym ujęciu: wybór metody kwadratury to w istocie decyzja, jak wydać „budżet” wywołań funkcji.

Popularne metody kwadratury

Poniższe metody uporządkowano od najbardziej intuicyjnych do najbardziej zaawansowanych.

Wizualne porównanie reguły punktu środkowego, trapezów i Simpsona. Obraz autora.

Reguła punktu środkowego

Weź punkt środkowy przedziału [a, b], oblicz tam f i pomnóż przez szerokość (b − a). Założenie jest takie, że wartość funkcji w środku reprezentuje średnią w całym przedziale.

Brzmi prawie zbyt prosto — i dla pojedynczego przedziału tak jest. Ale tu pierwsza niespodzianka: reguła punktu środkowego jest dokładna dla funkcji stałych i liniowych, a jej błąd w przeliczeniu na przedział skaluje się jak O(h³), gdzie h = b − a. To w rzeczywistości lepiej niż metoda trapezów dla gładkich funkcji, co często zaskakuje.

Metoda trapezów

Użyj obu końców zamiast tylko środka. Połącz je linią prostą i oblicz pole powstałego trapezu: (b − a) · [f(a) + f(b)] / 2. Dokładna dla funkcji liniowych, z błędem O(h³) na przedział — tym samym rzędem, co reguła punktu środkowego, choć ze większą stałą.

Metoda trapezów jest prawdopodobnie najczęściej nauczaną metodą kwadratury, co przyniosło jej nieco przesadzoną reputację. Dla dobrze zachowujących się funkcji jest przyzwoita. Ale złożona reguła Simpsona niemal zawsze wypada lepiej przy tej samej liczbie obliczeń.

Reguła Simpsona

Dopasuj wielomian kwadratowy przez trzy punkty: dwa końce i punkt środkowy. Całkuj ten wielomian dokładnie:

Dodatkowe obliczenie w środku kosztuje jedno wywołanie funkcji więcej, ale zysk jest realny: reguła Simpsona jest dokładna dla wielomianów do stopnia trzeciego, a błąd na przedział skaluje się jak O(h⁵), o dwa rzędy lepiej niż metoda trapezów. Dla gładkich funkcji to zwykle pierwsza metoda, po którą warto sięgnąć.

Kwadratura Gaussa

Tu zaczyna się robić naprawdę ciekawie. Reguła punktu środkowego, trapezy i Simpson używają równo rozmieszczonych węzłów. Kwadratura Gaussa stawia inne pytanie: co, jeśli dobralibyśmy zarówno położenia węzłów, jak i wagi tak, aby zmaksymalizować dokładność?

Dla n węzłów kwadratura Gaussa (konkretnie wariant Gaussa–Legendre’a) potrafi dokładnie całkować dowolny wielomian stopnia do 2n − 1. Dwa węzły dają dokładność dla wielomianów stopnia trzeciego; pięć węzłów obsługuje stopień 9. Węzły to miejsca zerowe wielomianów Legendre’a i skupiają się bliżej końców przedziału zamiast leżeć równomiernie.Dodatkowa dokładność wynika z możliwości przesuwania węzłów, a nie z tego, gdzie przypadkiem trafią. Ich skupienie przy końcach to skutek uboczny matematyki i poniżej zobaczymy dlaczego.

Tak wypada porównanie wszystkich czterech metod:

Złożone metody kwadratury

Zastosowanie reguły raz na całym przedziale [a, b] zwykle nie wystarcza. Standardowe rozwiązanie: podziel przedział na n równych podprzedziałów o szerokości h = (b − a)/n i zastosuj podstawową regułę w każdym z nich.

Złożona reguła punktu środkowego i metoda trapezów

Dla złożonej metody trapezów oblicz f w n+1 równomiernie rozmieszczonych punktach i zsumuj n trapezów. Całkowity błąd spada do O(h²). Gdy przetniesz h na pół, błąd maleje czterokrotnie. Złożona reguła punktu środkowego działa analogicznie i ma ten sam rząd błędu globalnego.

Złożona reguła Simpsona

Wersja złożona wymaga parzystej liczby podprzedziałów. Oblicz f w n+1 punktach i uzyskaj globalny błąd rzędu O(h⁴). Połowienie h zmniejsza błąd szesnastokrotnie. Dlatego złożona reguła Simpsona bywa domyślnym wyborem, gdy potrzebujesz umiarkowanej dokładności dla gładkiej funkcji.

Poniższy kod w Pythonie pokazuje, jak przybliżenie się poprawia wraz z dodawaniem podprzedziałów dla ∫[0 do π] sin(x) dx (wartość dokładna: 2.0):

import numpy as np

def composite_simpson(func, lower, upper, num_intervals):
    """Approximate integral using composite Simpson's rule."""
    if num_intervals % 2 != 0:
        raise ValueError("num_intervals must be even")
    nodes = np.linspace(lower, upper, num_intervals + 1)
    values = func(nodes)
    h = (upper - lower) / num_intervals
    weights = np.ones(num_intervals + 1)
    weights[1:-1:2] = 4   # odd-indexed interior points
    weights[2:-2:2] = 2   # even-indexed interior points
    return (h / 3) * np.dot(weights, values)

exact = 2.0
for n in [4, 8, 16, 32]:
    approx = composite_simpson(np.sin, 0, np.pi, n)
    print(f"n={n:3d}: approx={approx:.10f}, error={abs(approx - exact):.2e}")

Zaledwie 32 podprzedziały wystarczą, by błąd spadł poniżej 10⁻⁶. To zbieżność rzędu O(h⁴).

Kwadratura Gaussa

Idea stojąca za kwadraturą Gaussa zasługuje na więcej niż tylko wiersz w tabeli.

Dlaczego równomierne węzły nie zawsze są optymalne

Gdy węzły leżą w równych odstępach, wydajesz budżet obliczeń w sposób, który może nie odpowiadać trudności całkowanej funkcji. Niektóre obszary wymagają większej rozdzielczości; inne marnują obliczenia na płaskich fragmentach. Kwadratura Gaussa omija to, wybierając węzły, które są dowiedziono optymalne dla stopnia wielomianu, jaki można całkować dokładnie.

Dla n węzłów reguła Gaussa–Legendre’a wybiera n miejsc zerowych n-tego wielomianu Legendre’a Pₙ(x) jako węzły, a następnie oblicza wagi tak, by każdy wielomian stopnia ≤ 2n − 1 był całkowany bez błędu. W praktyce numpy.polynomial.legendre.leggauss(n) zwraca oba te elementy. Węzły skupiają się bliżej końców przedziału [−1, 1], standardowego przedziału odniesienia, gdzie błąd interpolacji wielomianowej bywa największy i gdzie dodatkowa gęstość próbkowania najbardziej się opłaca.

Kiedy jej używać

Kwadratura Gaussa jest najlepsza, gdy masz gładką, analityczną funkcję i potrzebujesz wysokiej dokładności przy niewielkiej liczbie wywołań funkcji. Pięć punktów Gaussa–Legendre’a może przewyższyć złożoną regułę Simpsona z dziesiątkami podprzedziałów dla dobrze zachowujących się całek. Słabo jednak nadaje się do funkcji z załamaniami, nieciągłościami lub danymi dostępnymi tylko w z góry zadanych, równych odstępach. W takich przypadkach starannie dobrane rozmieszczenie węzłów obróci się przeciwko tobie.

Błąd i dokładność kwadratury

Żadna metoda numeryczna nie jest wolna od błędów. Dwa odrębne źródła wpływają na każdy wynik kwadratury i zrozumienie obu jest kluczowe, by wiedzieć, kiedy ufać wynikowi.

Źródła błędów obejmują:

• Błąd dyskretyzacji (zwany też błędem obcięcia) to różnica między dokładną całką a wynikiem twojej reguły kwadratury przy idealnej arytmetyce. To właśnie opisują wzory O(hᵏ). Maleje, gdy dodajesz punkty lub zmniejszasz szerokości przedziałów.
• Błąd zmiennoprzecinkowy to inny problem. Sumowanie wielu małych wielkości w arytmetyce o skończonej precyzji wprowadza zaokrąglenia, które nie maleją wraz z h, a mogą wręcz rosnąć. Przy metodach złożonych z tysiącami podprzedziałów możesz dojść do progu, gdzie dominuje szum zmiennoprzecinkowy i dodawanie przedziałów przestaje pomagać. W większości zadań obliczeń naukowych osiągniesz wymaganą tolerancję dużo wcześniej. Ale taki próg istnieje.

Gładkość funkcji, szerokość przedziału i wybór metody wpływają na błąd dyskretyzacji. „Chropowata” funkcja z ostrymi cechami wymaga większej liczby przedziałów lub podejścia adaptacyjnego.

Kwadratura adaptacyjna

Jednolita siatka marnuje obliczenia, gdy całkowana funkcja jest gładka w większości miejsc, a skomplikowana w kilku.

Jak to działa

Zastosuj regułę kwadratury do przedziału, następnie zastosuj regułę wyższego rzędu (lub tę samą regułę na dwóch połowach) i porównaj. Różnica szacuje błąd lokalny. Jeśli jest poniżej tolerancji — zaakceptuj przedział. Jeśli nie — podziel go i rekurencyjnie działaj na każdej połowie. Powtarzaj, aż każdy podprzedział spełni tolerancję.

Zysk wydajności względem jednolitego zagęszczania może być duży, gdy funkcja ma odosobnione cechy: wąski pik, załamanie, szybką oscylację otoczoną gładkim obszarem. Metody adaptacyjne koncentrują obliczenia dokładnie tam, gdzie wymaga tego funkcja.

scipy.integrate.quad()

W Pythonie scipy.integrate.quad() implementuje adaptacyjną kwadraturę Gaussa (konkretnie algorytm DQAGSE z QUADPACK). Przekaż funkcję i granice; otrzymasz przybliżenie i oszacowanie błędu. Dla większości całek jednowymiarowych to właściwe narzędzie, bez ręcznego strojenia.

from scipy import integrate
import numpy as np

# Integrate a function with a peak near x=0.5
def peaked_func(x):
    return np.exp(-100 * (x - 0.5)**2)

result, error_bound = integrate.quad(peaked_func, 0, 1)
print(f"Estimate: {result:.10f}")
print(f"Error bound: {error_bound:.2e}")

quad() radzi sobie z tym wąskim pikiem bez problemu. Jednolita reguła Simpsona potrzebowałaby setek podprzedziałów, by dorównać.

Jak wybrać metodę kwadratury

Zacznij od scipy.integrate.quad() dla całek jednowymiarowych w Pythonie, chyba że masz powód, by tego nie robić. Jest adaptacyjna, dokładna i dobrze przetestowana. Wszystko poniżej jest istotne, gdy implementujesz coś samodzielnie, pracujesz w języku bez dobrej biblioteki adaptacyjnej albo optymalizujesz pod konkretny przypadek.

• Gładkie, analityczne funkcje o znanym wzorze: kwadratura Gaussa daje największą dokładność na wywołanie. Jeśli potrzebujesz umiarkowanej dokładności i prostoty, złożona reguła Simpsona to solidny domyślny wybór.
• Funkcje z załamaniami lub nieciągłościami: Metody z równymi odstępami (trapezy, Simpson) zwykle przewyższają tu kwadraturę Gaussa, bo stałe położenia węzłów nie są zaburzane przez nieprawidłowości. Podziel przedział całkowania w miejscu znanej nieciągłości i obsłuż każdą część osobno.
• Dane w stałych, równych odstępach: kwadratura Gaussa odpada. Nie możesz obliczać funkcji w dowolnych węzłach, więc użyj złożonej metody trapezów lub złożonej reguły Simpsona.
• Całki wielowymiarowe: Kwadratura szybko wpada w klątwę wymiarowości. Przy n punktach na wymiar i d wymiarach potrzebujesz łącznie nᵈ obliczeń. Dziesięć wymiarów po dziesięć punktów to 10¹⁰ obliczeń — nie do zrobienia. Monte Carlo skaluje się znacznie lepiej powyżej ~5 wymiarów, więc tam następuje zmiana metody.

Kwadratura a pokrewne metody numeryczne

Kwadratura nie istnieje w próżni. Warto zobaczyć, gdzie mieści się względem innych narzędzi, które prawdopodobnie rozważasz.

Kwadratura a całkowanie analityczne

Całkowanie analityczne znajduje dokładną funkcję pierwotną; kwadratura przybliża liczbę. Jeśli istnieje dokładna funkcja pierwotna i potrafisz ją znaleźć — użyj jej. Kwadratura ma sens, gdy całkowanie dokładne jest niemożliwe lub niepraktyczne.

Kwadratura a całkowanie Monte Carlo

Obie metody przybliżają całki oznaczone, ale w inny sposób. Kwadratura jest deterministyczna: oblicza f w stałych, z góry wybranych węzłach, a zbieżność wiąże się z gładkością funkcji i liczbą węzłów. Monte Carlo losuje próbki i je uśrednia, z błędem skaluącym się jak O(1/√N) niezależnie od wymiaru. To wolniejsza zbieżność w niskich wymiarach, ale niezależna od wymiarowości, co czyni ją praktycznym wyborem powyżej około pięciu wymiarów.

Kwadratura a interpolacja

Niektóre reguły kwadratury, w tym reguła Simpsona, działają przez dopasowanie wielomianu do punktów obliczeń i jego dokładne całkowanie. Cel jednak różni się od interpolacji. Interpolacja odtwarza funkcję; kwadratura chce tylko liczby. Wielomian jest narzędziem, nie wynikiem.

Zakończenie

W kwadraturze satysfakcjonuje mnie to, jak wiele daje jedna porządkująca idea. Każda z metod — reguła punktu środkowego, trapezy, Simpson, Gauss — to w gruncie rzeczy inna odpowiedź na to samo pytanie: mając budżet wywołań funkcji, gdzie je umieścić i jak je zważyć?

Lepsze dopasowanie wielomianów przez punkty obliczeń szybko obniża błąd. Przesunięcie do optymalnie dobranych węzłów wyciska jeszcze więcej dokładności z tego samego budżetu. Metody adaptacyjne spajają to wszystko, kierując obliczenia tam, gdzie funkcja rzeczywiście ich potrzebuje.

W Pythonie scipy.integrate.quad() to właściwy domyślny wybór dla problemów 1D. Reszta artykułu przyda się, gdy chcesz zrozumieć, co dzieje się „pod maską”, albo gdy implementujesz coś samodzielnie. Jeśli chcesz dalej rozwijać umiejętności numeryczne z SciPy, zajrzyj na nasz kurs Introduction to Optimization in Python.