Czy w przypadku danych zaszumionych lub o drobnej strukturze są lepsze kwadratury niż reguła punktu środkowego?

Niezbędne są tylko dwie pierwsze części tego długiego pytania. Pozostałe służą wyłącznie ilustracji.

tło

Zaawansowane kwadratury, takie jak kompozyt wyższego stopnia Newtona-Cotesa, Gaußa-Legendre'a i Romberga, wydają się być głównie przeznaczone do przypadków, w których można dokładnie próbkować funkcję, ale nie integrować analitycznie. Jednak w przypadku funkcji o strukturach większych niż interwał próbkowania (patrz przykład A dodatek) lub szumu pomiarowego nie mogą one konkurować z prostymi podejściami, takimi jak punkt środkowy lub reguła trapezowa (demonstracja znajduje się w dodatku B).

Jest to nieco intuicyjne, ponieważ np. Złożona reguła Simpsona zasadniczo „odrzuca” jedną czwartą informacji, przypisując jej niższą wagę. Jedynym powodem, dla którego takie kwadratury są lepsze dla wystarczająco nudnych funkcji, jest to, że prawidłowe zarządzanie efektami granicznymi przeważa nad efektem odrzucenia informacji. Z innego punktu widzenia intuicyjnie jest dla mnie jasne, że w przypadku funkcji o drobnej strukturze lub szumie próbki, które są oddalone od granic domeny integracji, muszą być prawie w równej odległości i mieć prawie taką samą wagę (w przypadku dużej liczby próbek ). Z drugiej strony kwadratura takich funkcji może skorzystać na lepszej obsłudze efektów granicznych (niż w metodzie punktu środkowego).

Pytanie

Załóżmy, że chcę zintegrować numerycznie jednowymiarowe dane z zaszumieniem lub drobną strukturą.

Liczba punktów próbkowania jest stała (ze względu na kosztowną ocenę funkcji), ale mogę je dowolnie umieszczać. Jednak ja (lub metoda) nie mogę umieszczać punktów próbkowania interaktywnie, tj. W oparciu o wyniki z innych punktów próbkowania. Nie znam też wcześniej potencjalnych regionów problemowych. Tak więc coś w stylu Gauß – Legendre (nierównoodległe punkty próbkowania) jest w porządku; kwadratura adaptacyjna nie jest, ponieważ wymaga interaktywnie rozmieszczonych punktów próbkowania.

• Czy w takim przypadku zaproponowano metody wykraczające poza metodę punktu środkowego?

• Lub: Czy istnieje dowód, że metoda punktu środkowego jest najlepsza w takich warunkach?

• Mówiąc bardziej ogólnie: czy są jakieś prace nad tym problemem?

Dodatek A: Konkretny przykład funkcji o drobnej strukturze

Chcę oszacować dla: z i . Typowa funkcja wygląda następująco:$\int_0^1f(t)\, \mathrm{d}t$φi∈[0,2π]logωi∈[1,1000

$$f(t) = \sum_{i=1}^{k} \frac{\sin(ω_i t-φ_i)}{ω_i},$$ $φ_i∈ [0,2π]$$\log{ω_i} ∈ [1,1000]$

Wybrałem tę funkcję dla następujących właściwości:

• Można go zintegrować analitycznie w celu uzyskania wyniku kontroli.
• Ma drobną strukturę na poziomie, który uniemożliwia uchwycenie wszystkiego przy użyciu liczby próbek, których używam ( ).$<10^2$
• Nie jest zdominowany przez swoją delikatną strukturę.

Załącznik B: Benchmark

Dla kompletności, oto test porównawczy w Pythonie:

import numpy as np
from numpy.random import uniform
from scipy.integrate import simps, trapz, romb, fixed_quad

begin = 0
end   = 1

def generate_f(k,low_freq,high_freq):
    ω = 2**uniform(np.log2(low_freq),np.log2(high_freq),k)
    φ = uniform(0,2*np.pi,k)
    g = lambda t,ω,φ: np.sin(ω*t-φ)/ω
    G = lambda t,ω,φ: np.cos(ω*t-φ)/ω**2
    f = lambda t: sum( g(t,ω[i],φ[i]) for i in range(k) )
    control = sum( G(begin,ω[i],φ[i])-G(end,ω[i],φ[i]) for i in range(k) )
    return control,f

def midpoint(f,n):
    midpoints = np.linspace(begin,end,2*n+1)[1::2]
    assert len(midpoints)==n
    return np.mean(f(midpoints))*(n-1)

def evaluate(n,control,f):
    """
    returns the relative errors when integrating f with n evaluations
    for several numerical integration methods.
    """
    times = np.linspace(begin,end,n)
    values = f(times)
    results = [
            midpoint(f,n),
            trapz(values),
            simps(values),
            romb (values),
            fixed_quad(f,begin,end,n=n)[0]*(n-1),
        ]

    return [
            abs((result/(n-1)-control)/control)
            for result in results
        ]

method_names = ["midpoint","trapezoid","Simpson","Romberg","Gauß–Legendre"]

def med(data):
    medians = np.median(np.vstack(data),axis=0)
    for median,name in zip(medians,method_names):
        print(f"{median:.3e}   {name}")

print("superimposed sines")
med(evaluate(33,*generate_f(10,1,1000)) for _ in range(100000))

print("superimposed low-frequency sines (control)")
med(evaluate(33,*generate_f(10,0.5,1.5)) for _ in range(100000))

(Używam tutaj mediany, aby zmniejszyć wpływ wartości odstających z powodu funkcji, które mają tylko zawartość wysokich częstotliwości. Dla średniej wyniki są podobne.)

Mediany względnych błędów integracji to:

superimposed sines
6.301e-04   midpoint
8.984e-04   trapezoid
1.158e-03   Simpson
1.537e-03   Romberg
1.862e-03   Gauß–Legendre

superimposed low-frequency sines (control)
2.790e-05   midpoint
5.933e-05   trapezoid
5.107e-09   Simpson
3.573e-16   Romberg
3.659e-16   Gauß–Legendre

Uwaga: po dwóch miesiącach i jednej nagrodie bez rezultatu opublikowałem to na MathOverflow .

Po pierwsze, myślę, że źle rozumiesz pojęcie kwadratury adaptacyjnej. Kwadratura adaptacyjna nie oznacza „interaktywnego umieszczania punktów próbki”. Cała idea kwadratury adaptacyjnej polega na opracowaniu schematu, który zintegruje określoną funkcję z określonym (szacowanym) błędem bezwzględnym lub względnym przy możliwie najmniejszej ocenie funkcji.

Druga uwaga: piszesz „Liczba punktów próbkowania jest stała (ze względu na kosztowną ocenę funkcji), ale mogę je dowolnie umieszczać”. Myślę, że powinna istnieć idea, aby liczba punktów próbkowania (lub ocen funkcji w terminologii kwadraturowej) była jak najmniejsza (tj. Nie ustalona).

Jaki jest więc pomysł na kwadraturę adaptacyjną zaimplementowaną na przykład w QUADPACK ?

1. Podstawowym składnikiem jest „zagnieżdżona” zasada kwadraturowa: jest to kombinacja dwóch reguł kwadraturowych, w których jedna ma wyższą kolejność (lub dokładność) jak druga. Dlaczego? Na podstawie różnicy między tymi regułami algorytm może oszacować błąd kwadraturowy (oczywiście algorytm użyje najbardziej dokładnego jako wyniku odniesienia). Przykładami mogą być reguła trapezowa z węzłami i węzłami. W przypadku QUADPACK są to reguły Gaussa-Kronroda. Są to interpolacyjne reguły kwadraturowe, które wykorzystują regułę kwadraturową Gaussa-Legendre'a określonego rzędu 2 n + 1 N N 2 N - $2^{n}$$2^{n+1}$$N$i optymalne rozszerzenie tej zasady. Oznacza to, że można uzyskać wyższą kolejność kwadraturową, ponownie wykorzystując węzły Gaussa-Legendre'a (tj. Kosztowne oceny funkcji) o różnych wagach i dodając szereg dodatkowych węzłów. Innymi słowy, oryginalna reguła Gaussa-Legendre'a rzędu zintegruje dokładnie wszystkie wielomiany stopnia dokładnie, podczas gdy rozszerzona reguła Gaussa-Kronroda dokładnie zintegruje jakiś wielomian wyższego rzędu. Klasyczną zasadą jest G7K15 (Gauss-Legendre 7. rzędu z Gauss-Kronrod 15. rzędu). Magia polega na tym, że 7 węzłów Gaussa-Legendre'a jest podzbiorem 15 węzłów Gaussa-Kronroda, więc z 15 ocenami funkcji mam ocenę kwadratury wraz z oszacowaniem błędu!$N$$2N-1$

2. Kolejnym składnikiem jest strategia „dziel i rządź”. Załóżmy, że puściłeś G7K15 na swoim integrandzie i zaobserwowałeś błąd kwadratury, który według twojego gustu jest zbyt duży. QUADPACK podzieli pierwotny interwał na dwie równomiernie rozmieszczone podinterwały. A następnie dokona ponownej oceny dwóch subintegrałów przy użyciu podstawowej zasady, G7K15. Teraz algorytm ma globalne oszacowanie błędów (które powinno być, być może, niższe niż pierwsze), ale także dwa lokalne oszacowania błędów. Wybiera przedział z największym błędem i dzieli go na dwa. Oszacowano dwie nowe całki i zaktualizowano błąd globalny. I tak dalej, dopóki błąd globalny nie spadnie poniżej żądanego celu lub nie zostanie przekroczona maksymalna liczba podziałów.

Wzywam więc do aktualizacji powyższego kodu za pomocą tej scipy.quadmetody. Być może w przypadku integrandu z dużą liczbą „drobnych struktur” może być konieczne zwiększenie maksymalnej liczby podziałów ( limitopcja). Możesz także grać z parametrami epsabsi / lub epsrel.

Jeśli jednak masz tylko dane eksperymentalne, widzę dwie możliwości.

1. Jeśli masz możliwość wybrania punktów pomiarowych, tj. Wartości , wybrałbym je jednakowo i najlepiej jako potęgę , abyś mógł zastosować zagnieżdżoną regułę trapezoidalną (i czerpać zyski z ekstrapolacji Romberga).$t$$2$
2. Jeśli nie masz możliwości wyboru węzłów, tzn. Pomiary są wykonywane losowo, najlepszą moim zdaniem nadal jest reguła trapezowa.

Nie jestem przekonany, że twój kod pokazuje coś fundamentalnego w różnych regułach kwadratury i jak dobrze radzą sobie z hałasem i drobną strukturą, i wierzę, że jeśli wybierzesz różne struktury drobnych kar, znajdziesz coś innego. Oto twierdzenie:

Żadna metoda kwadraturowa nie może dać niskiego błędu bezwzględnego lub względnego względem funkcji z nieograniczoną zmiennością całkowitą. W systemie zmiennoprzecinkowym z zaokrągleniem jednostek mamy oszacowanie$\mu$ gdzie jest sumą kwadraturową działającą na implementacji numerycznej z .

$$\left| \int_{a}^{b} f \, \mathrm{d}x - \hat{Q}[\hat{f}] \right| \le \left| \int_{a}^{b} f \, \mathrm{d}x - Q[f] \right| + \mu\left[ 4\int_{a}^{b} |f| \, \mathrm{d}x + \int_{a}^{b} |xf'| \, \mathrm{d}x \right]$$ P$\hat{Q}$f f$\hat{f}$$f$

Dowód: niech węzły kwadraturowe będą mieć a (nieujemne) wagi kwadraturowe będą wynosić i określają przybliżone liczby zmiennoprzecinkowe przez i . Załóżmy, że spełnia gdzie gdzie jest jednostkowym zaokrągleniem. Następnie $\{x_i\}_{i=0}^{n-1}$$\{w_i\}_{i=0}^{n-1}$$\hat{w}_{i}$$\hat{x}_i$$\hat{f}$$\hat{f}(x) = f(x)(1+2\delta)$$|\delta| \le \mu$$\mu$

$$\begin{align*} \hat{Q}[\hat{f}] &= \sum_{i=0}^{n-1} \hat{w}_i \otimes \hat{f}(\hat{x}_i) \\ &= \sum_{i=0}^{n-1} w_i (1+\delta^w_i)f(x_i + \delta_i^x x_i)(1+2\delta_i^{f})(1+\delta_i^{*}) \\ &\approx \sum_{i=0}^{n-1} w_i \left[f(x_i)+ \delta_i^x x_i f'(x_i) \right] (1+\delta^w_i + 2\delta_i^{f} + \delta_i^*) \\ &\approx \sum_{i=0}^{n-1} w_i f(x_i) + \sum_{i=0}^{n-1}\delta_i^x w_i x_i f'(x_i) + w_i f(x_i) (\delta^w_i + 2\delta_i^{f} + \delta_i^*) \\ \end{align*}$$ , aby Zakłada się, że suma jest obliczana bez błędów; pomnóż przez aby odrzucić to założenie. $$\begin{align*} |\hat{Q}[\hat{f}] - Q[f]| &\le \mu \sum_{i=0}^{n-1}w_i(|x_i f'(x_i)| + 4|f(x_i)|) \\ &\approx 4\mu \int |f| \, \mathrm{d}x + \mu \int |xf'| \, \mathrm{d}x \end{align*}$$ $n$

Mutatis mutandis można również pokazać, że wynik jest zachowany w arytmetyki punktu stałego.