Untuk data yang bising atau terstruktur dengan baik, apakah ada kuadratur yang lebih baik daripada aturan titik tengah?

Hanya dua bagian pertama dari pertanyaan panjang ini yang penting. Yang lain hanya untuk ilustrasi.

Latar Belakang

Kuadratur maju seperti komposit tingkat tinggi Newton-Cotes, Gauß-Legendre, dan Romberg tampaknya terutama ditujukan untuk kasus-kasus di mana seseorang dapat mengambil sampel fungsi secara halus tetapi tidak berintegrasi secara analitis. Namun, untuk fungsi dengan struktur yang lebih halus daripada interval pengambilan sampel (lihat Lampiran A untuk contoh) atau noise pengukuran, mereka tidak dapat bersaing dengan pendekatan sederhana seperti titik tengah atau aturan trapesium (lihat Lampiran B untuk demonstrasi).

Ini agak intuitif karena, misalnya, aturan Simpson komposit pada dasarnya "membuang" seperempat informasi dengan menetapkan bobot yang lebih rendah. Satu-satunya alasan kuadratur semacam itu lebih baik untuk fungsi-fungsi yang cukup membosankan adalah penanganan efek perbatasan yang lebih baik daripada efek informasi yang dibuang. Dari sudut pandang lain, secara intuitif jelas bagi saya bahwa untuk fungsi dengan struktur atau noise yang halus, sampel yang jauh dari batas domain integrasi harus hampir sama dan memiliki bobot yang hampir sama (untuk jumlah sampel yang tinggi). ). Di sisi lain, kuadratur fungsi-fungsi tersebut dapat mengambil manfaat dari penanganan efek perbatasan yang lebih baik (daripada untuk metode titik tengah).

Pertanyaan

Asumsikan bahwa saya ingin secara numerik mengintegrasikan data satu dimensi yang bising atau terstruktur dengan baik.

Jumlah titik pengambilan sampel ditetapkan (karena evaluasi fungsi menjadi mahal), tetapi saya dapat dengan bebas menempatkannya. Namun, saya (atau metode) tidak dapat menempatkan titik pengambilan sampel secara interaktif, yaitu berdasarkan hasil dari titik pengambilan sampel lainnya. Saya juga tidak tahu daerah masalah potensial sebelumnya. Jadi, sesuatu seperti Gauß – Legendre (titik pengambilan sampel yang tidak sama) tidak apa-apa; quadrature adaptif bukan karena membutuhkan titik pengambilan sampel yang ditempatkan secara interaktif.

• Apakah ada metode yang melampaui metode titik tengah yang disarankan untuk kasus seperti itu?

• Atau: Apakah ada bukti bahwa metode titik tengah terbaik dalam kondisi seperti itu?

• Secara umum: Apakah ada pekerjaan yang ada untuk masalah ini?

Lampiran A: Contoh spesifik dari fungsi terstruktur halus

Saya ingin memperkirakan untuk: dengan dan . Fungsi khas terlihat seperti ini:$\int_0^1f(t)\, \mathrm{d}t$

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

Saya memilih fungsi ini untuk properti berikut:

• Ini dapat diintegrasikan secara analitis untuk hasil kontrol.
• Ini memiliki struktur yang bagus pada tingkat yang membuatnya tidak mungkin untuk menangkap semuanya dengan jumlah sampel yang saya gunakan ( ).$<10^2$
• Itu tidak didominasi oleh strukturnya yang halus.

Lampiran B: Tolok Ukur

Untuk kelengkapan, berikut adalah patokan dalam Python:

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))

(Saya di sini menggunakan median untuk mengurangi pengaruh outlier karena fungsi yang hanya memiliki konten frekuensi tinggi. Untuk rata-rata, hasilnya serupa.)

Median kesalahan integrasi relatif adalah:

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

Catatan: Setelah dua bulan dan satu hadiah tanpa hasil, saya memposting ini di MathOverflow .

Pertama-tama, saya pikir Anda salah memahami konsep quadrature adaptif. Kuadratur adaptif tidak menyiratkan "menempatkan titik sampel secara interaktif". Seluruh ide di balik quadrature adaptif adalah untuk merancang skema yang akan mengintegrasikan fungsi tertentu ke kesalahan absolut atau relatif (diperkirakan) tertentu dengan evaluasi fungsi sesedikit mungkin.

Komentar kedua: Anda menulis "Jumlah titik pengambilan sampel ditetapkan (karena evaluasi fungsi menjadi mahal), tetapi saya dapat dengan bebas menempatkannya". Saya pikir idenya harus bahwa jumlah titik pengambilan sampel (atau evaluasi fungsi dalam terminologi quadrature) harus sekecil mungkin (yaitu tidak tetap).

Jadi apa ide di balik quadrature adaptif seperti yang diterapkan di QUADPACK misalnya?

1. Bahan dasar adalah aturan quadrature "bersarang": ini adalah kombinasi dari dua aturan quadrature di mana satu memiliki urutan yang lebih tinggi (atau akurasi) sebagai yang lain. Mengapa? Berdasarkan perbedaan antara aturan-aturan ini, algoritma dapat memperkirakan kesalahan quadrature (tentu saja, algoritma akan menggunakan yang paling akurat sebagai hasil referensi). Contohnya bisa berupa aturan trapesium dengan node dan node. Dalam kasus QUADPACK, aturannya adalah aturan Gauss-Kronrod. Ini adalah aturan kuadratur interpolasi yang menggunakan aturan kuadratur Gauss-Legendre dari pesanan tertentu 2 n + 1 N N 2 N - $2^{n}$$2^{n+1}$$N$dan perpanjangan optimal dari aturan ini. Ini berarti bahwa urutan quadrature yang lebih tinggi dapat diperoleh dengan menggunakan kembali node Gauss-Legendre (yaitu evaluasi fungsi yang mahal) dengan bobot yang berbeda dan menambahkan sejumlah node tambahan. Dengan kata lain, aturan orde Gauss-Legendre orde akan mengintegrasikan semua polinomial derajat persis sementara aturan Gauss-Kronrod yang diperluas akan mengintegrasikan beberapa polinomial orde tinggi dengan tepat. Aturan klasik adalah G7K15 (urutan 7 Gauss-Legendre dengan urutan 15 Gauss-Kronrod). Keajaibannya adalah bahwa 7 simpul Gauss-Legendre adalah bagian dari 15 simpul Gauss-Kronrod sehingga dengan 15 evaluasi fungsi, saya memiliki evaluasi quadrature bersama dengan perkiraan kesalahan!$N$$2N-1$

2. Bahan selanjutnya adalah strategi "memecah-dan-menaklukkan". Misalkan Anda melepaskan G7K15 ini pada integrand Anda dan Anda melihat kesalahan quadrature yang menurut selera Anda terlalu besar. QUADPACK kemudian akan membagi interval asli dalam dua sub-interval dengan jarak yang sama. Dan kemudian akan mengevaluasi kembali dua subintegral menggunakan aturan dasar, G7K15. Sekarang, algoritma memiliki estimasi kesalahan global (yang semestinya lebih rendah dari yang pertama) tetapi juga dua perkiraan kesalahan lokal. Ini mengambil interval dengan kesalahan terbesar dan membaginya menjadi dua. Diperkirakan dua integral baru dan kesalahan global diperbarui. Dan seterusnya sampai kesalahan global di bawah target yang Anda minta atau jumlah maksimum subdivisi telah dilampaui.

Jadi saya menantang Anda untuk memperbarui kode Anda di atas menggunakan scipy.quadmetode ini. Mungkin dalam kasus integand dengan banyak "struktur halus" Anda mungkin perlu meningkatkan jumlah maksimum subdivisi ( limitopsi). Anda juga bisa bermain dengan epsabsdan / atau epsrelparameter.

Namun, jika Anda hanya memiliki data eksperimental, saya melihat dua kemungkinan.

1. Jika Anda memiliki kesempatan untuk memilih titik pengukuran, yaitu nilai , saya akan memilihnya secara sama dan lebih disukai sebagai kekuatan sehingga Anda dapat menerapkan aturan trapesium bersarang (dan mendapat keuntungan dari ekstrapolasi Romberg).$t$$2$
2. Jika Anda tidak memiliki sarana untuk memilih node, yaitu pengukuran datang secara acak, opsi terbaik menurut saya masih adalah aturan trapesium.

Saya tidak yakin bahwa kode Anda menunjukkan sesuatu yang mendasar tentang berbagai aturan kuadratur dan seberapa baik mereka melawan kebisingan dan struktur halus, dan percaya bahwa jika Anda memilih berbagai struktur denda yang berbeda, Anda akan menemukan sesuatu yang berbeda. Inilah teorema:

Tidak ada metode quadrature yang dapat memberikan kesalahan absolut atau relatif rendah terhadap suatu fungsi dengan variasi total tanpa batas. Dalam sistem floating point dengan unit roundoff , kami memiliki perkiraan$\mu$ mana adalah jumlah kuadratur yang bekerja pada implementasi numerik dari .

$$\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]$$ Q$\hat{Q}$f f$\hat{f}$$f$

Bukti: Biarkan simpul quadrature menjadi dan bobot kuadratur (non-negatif) menjadi dan menunjukkan aproksimasi titik apung mereka dengan dan . Asumsikan memenuhi mana mana adalah unit roundoff. Kemudian $\{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*}$$ sehingga Ini mengasumsikan bahwa jumlah dihitung tanpa kesalahan; kalikan dengan untuk menjatuhkan asumsi itu. $$\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 Anda juga dapat menunjukkan bahwa hasilnya memegang aritmatika titik tetap.