Cvadra­tura explicată: metode numerice pentru aproximarea integralelor

Există o categorie de integrale definite care pur și simplu nu cedează în fața calculusului cu creionul pe hârtie. Nu pentru că greșești tu cu ceva. Unele funcții, precum e^(−x²) sau sin(x)/x, nu au primitive în formă închisă. Când se întâmplă asta, ai nevoie de o abordare diferită: evaluezi funcția în câteva puncte atent alese și combini acele valori pentru a estima aria de dedesubt. Asta este cvadra­tura.

Acest articol trece în revistă ce este cvadra­tura, cum funcționează cadrul general și principalele metode, de la intuitiva regulă a punctului mijlociu până la abordarea Gaussiană surprinzător de precisă. Vom acoperi și metode compuse, surse de eroare, strategii adaptive și cum să alegi metoda potrivită pentru problema pe care o ai în față.

Ce este cvadra­tura?

Cuvântul „cvadra­tură” e vechi. Vine din problema geometrică clasică de a construi un pătrat cu aceeași arie ca o formă dată, adică, în altă formulare, „găsește aria”. Astăzi, termenul înseamnă aproximarea numerică a unei integrale definite.

Ideea de bază este o substituție: în loc să găsești o primitivă exactă, eșantionezi funcția f(x) într-un set de puncte numite noduri, înmulțești fiecare valoare cu o greutate asociată și însumezi rezultatele:

Termenii cvadra­tură și integrare numerică înseamnă, în esență, același lucru. Cvadra­tura apare mai des în contexte matematice; integrarea numerică e termenul uzual în calculul științific.

Cum funcționează cvadra­tura numerică

Fiecare regulă de cvadra­tură operează pe un interval [a, b]. În interiorul acelui interval, selectează un set de noduri xᵢ și atribuie fiecăruia o greutate wᵢ. Aproximația este suma ponderată de mai sus. Ce diferă între metode este unde sunt plasate nodurile și cum sunt stabilite greutățile. Atât. Regula trapezelor, regula Simpson și cvadra­tura Gaussiană sunt toate variații pe același șablon.

Eroarea de aproximare este diferența dintre suma ponderată și integrala adevărată. Fiecare metodă introduce ceva eroare, iar cât de repede scade acea eroare pe măsură ce adaugi noduri (sau micșorezi intervalele) este ceea ce desparte o metodă bună de una excelentă. Ține minte această perspectivă: alegerea unei metode de cvadra­tură înseamnă, de fapt, alegerea modului în care îți „cheltui” evaluările funcției.

Metode comune de cvadra­tură

Metodele de mai jos merg de la cele mai intuitive la cele mai sofisticate.

Comparație vizuală a regulilor de integrare Midpoint, Trapezoidal și Simpson. Imagine de autor.

Regula punctului mijlociu

Ia punctul mijlociu al lui [a, b], evaluează f acolo și înmulțește cu lățimea (b − a). Presupunerea este că valoarea funcției în centru reprezintă media pe întregul interval.

Sună aproape prea simplu și, pentru un singur interval, chiar este. Dar iată prima surpriză: regula punctului mijlociu este exactă pentru funcții constante și liniare, iar eroarea ei scalează ca O(h³) pe interval, unde h = b − a. Asta este, de fapt, mai bine decât regula trapezelor pe funcții netede, lucru care tinde să-i ia pe nepregătite pe mulți.

Regula trapezelor

Folosește ambele capete în loc de doar centru. Le unești cu o linie dreaptă și calculezi aria trapezului rezultat: (b − a) · [f(a) + f(b)] / 2. Exactă pentru funcții liniare, cu eroare O(h³) pe interval, de același ordin ca regula punctului mijlociu, deși cu o constantă mai mare.

Regula trapezelor este probabil cea mai predată metodă de cvadra­tură, ceea ce i-a dat o reputație disproporționată. Pentru funcții bine comportate e decentă. Dar regula Simpson compusă aproape întotdeauna face mai bine cu același număr de evaluări.

Regula Simpson

Potrivești un polinom de gradul doi prin trei puncte: cele două capete și punctul mijlociu. Integrezi acel polinom exact:

Evaluarea suplimentară în centru costă un apel de funcție în plus, dar merită: regula Simpson este exactă pentru polinoame până la gradul trei, iar eroarea pe interval scalează ca O(h⁵), cu două ordine mai bună decât regula trapezelor. Pentru funcții netede, aceasta este de obicei prima metodă la care merită să apelezi.

Cvadra­tura Gaussiană

Aici lucrurile devin cu adevărat interesante. Regulile punctului mijlociu, trapezelor și Simpson folosesc noduri distanțate uniform. Cvadra­tura Gaussiană pune o altă întrebare: ce-ar fi dacă am alege atât pozițiile nodurilor, cât și greutățile pentru a maximiza acuratețea?

Cu n noduri, cvadra­tura Gaussiană (în special varianta Gauss-Legendre) poate integra exact orice polinom de grad până la 2n − 1. Două noduri îți dau exactitate pentru cubice; cinci noduri acoperă polinoame de gradul 9. Nodurile sunt rădăcinile polinoamelor lui Legendre și se aglomerează spre capetele intervalului, mai degrabă decât să stea la poziții echidistante. Precizia suplimentară vine din libertatea de a muta nodurile, nu din locul exact unde ajung. Aglomerarea lor aproape de capete este un produs secundar al matematicii, iar mai jos vom vedea de ce.

Iată cum se compară toate cele patru metode:

Metode compuse de cvadra­tură

Aplicarea unei reguli o singură dată pe întreg intervalul [a, b] de obicei nu este suficientă. Remediul standard: împarți intervalul în n subintervale egale de lățime h = (b − a)/n și aplici regula de bază pe fiecare.

Regulile compuse ale punctului mijlociu și trapezelor

Pentru regula trapezelor compusă, evaluezi f în n+1 puncte distanțate uniform și însumezi n trapeze. Eroarea totală scade la O(h²). Înjumătățești h, iar eroarea scade de patru ori. Regula punctului mijlociu compusă urmează aceeași logică și are același ordin global de eroare.

Regula Simpson compusă

Versiunea compusă necesită un număr par de subintervale. Evaluezi f în n+1 puncte și obții eroare globală O(h⁴). Înjumătățirea lui h reduce eroarea de șaisprezece ori. De aceea, regula Simpson compusă este adesea implicitul atunci când ai nevoie de acuratețe moderată pe o funcție netedă.

Codul Python de mai jos arată cum se strânge aproximația pe măsură ce adaugi subintervale pentru ∫[0 la π] sin(x) dx (valoare exactă: 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}")

Cu doar 32 de intervale, eroarea scade sub 10⁻⁶. Asta înseamnă convergență O(h⁴).

Cvadra­tura Gaussiană

Ideea din spatele cvadra­turii Gaussiene merită mai mult decât o intrare într-un tabel.

De ce nodurile distanțate uniform nu sunt mereu optime

Când nodurile sunt plasate în puncte distanțate uniform, îți cheltuiești bugetul de evaluări într-un mod care poate să nu se potrivească cu dificultatea integrandului. Unele regiuni au nevoie de mai multă rezoluție; altele irosesc evaluări pe porțiuni plate. Cvadra­tura Gaussiană evită asta alegând noduri care sunt demonstrabil optime pentru gradul de polinom pe care îl poți integra exact.

Pentru n noduri, regula Gauss-Legendre selectează cele n rădăcini ale polinomului Legendre de ordin n, Pₙ(x), ca noduri, apoi calculează greutăți astfel încât orice polinom de grad ≤ 2n − 1 să fie integrat fără eroare. În practică, numpy.polynomial.legendre.leggauss(n) le returnează pe ambele. Nodurile se aglomerează spre capetele intervalului [−1, 1], intervalul de referință standard, acolo unde eroarea de interpolare polinomială tinde să fie mai mare și unde densitatea suplimentară a eșantionării aduce cel mai mare câștig.

Când să o folosești

Cvadra­tura Gaussiană e la cel mai bun raport când ai o funcție netedă, analitică și ai nevoie de acuratețe ridicată cu puține evaluări. Cinci puncte Gauss-Legendre pot depăși regula Simpson compusă cu zeci de intervale pe integrande bine comportate. Dar e nepotrivită pentru funcții cu frânturi, discontinuități sau date disponibile doar în puncte fixate, echidistante. În astfel de cazuri, plasarea atentă a nodurilor lucrează împotriva ta.

Eroare și acuratețe în cvadra­tură

Nicio metodă numerică nu primește scutire de la eroare. Două surse distincte afectează fiecare rezultat de cvadra­tură, iar înțelegerea amândurora e esențială ca să știi când să ai încredere în răspuns.

Sursele de eroare includ:

• Eroarea de discretizare (numită și eroare de trunchiere) este diferența dintre integrala exactă și ceea ce returnează regula ta de cvadra­tură în aritmetică perfectă. Asta surprind formulele O(hᵏ). Scade pe măsură ce adaugi mai multe puncte sau reduci lățimea intervalului.
• Eroarea în virgulă mobilă este o altă problemă. Suma multor cantități mici în aritmetică cu precizie finită introduce rotunjiri care nu scad cu h și chiar pot crește. Cu reguli compuse care folosesc mii de subintervale, poți atinge un prag în care zgomotul numeric domină și adăugarea de intervale nu mai ajută. Pentru majoritatea sarcinilor de calcul științific, vei atinge toleranța dorită cu mult înainte de acel prag. Dar el există.

Netedimea funcției, lățimea intervalului și alegerea metodei contribuie toate la eroarea de discretizare. O funcție neregulată, cu variații abrupte, are nevoie de mai multe intervale sau de o abordare adaptivă.

Cvadra­tură adaptivă

O plasare uniformă a punctelor irosește evaluări când integrandul e neted în majoritatea locurilor, dar complicat în câteva.

Cum funcționează

Aplici o regulă de cvadra­tură pe un interval, apoi aplici o regulă de ordin mai înalt (sau aceeași regulă pe două jumătăți) și compari. Diferența estimează eroarea locală. Dacă e sub toleranță, accepți intervalul. Dacă nu, îl împarți și refaci recursiv pe fiecare jumătate. Repetă până când fiecare subinterval respectă toleranța.

Câștigul de eficiență față de rafinarea uniformă poate fi mare atunci când integrandul are trăsături izolate: un vârf îngust, o frântură, o oscilație rapidă înconjurată de zone netede. Metodele adaptive concentrează evaluările exact acolo unde funcția le cere.

scipy.integrate.quad()

În Python, scipy.integrate.quad() implementează cvadra­tura Gaussiană adaptivă (mai exact algoritmul DQAGSE din QUADPACK). Îi pasezi funcția și limitele; primești înapoi estimarea și o limită a erorii. Pentru majoritatea integralelor unidimensionale, acesta este instrumentul potrivit, fără reglaje manuale.

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() gestionează curat acest vârf îngust. O regulă Simpson uniformă ar avea nevoie de sute de subintervale ca să se apropie.

Cum alegi o metodă de cvadra­tură

Pornește cu scipy.integrate.quad() pentru integrale unidimensionale în Python, dacă nu ai un motiv să nu o faci. E adaptivă, precisă și bine testată. Tot ce urmează este relevant când implementezi tu însuți ceva, lucrezi într-un limbaj fără o bibliotecă adaptivă bună sau optimizezi pentru un caz specific.

• Funcții netede, analitice, cu formulă cunoscută: cvadra­tura Gaussiană oferă cea mai mare acuratețe per evaluare. Dacă ai nevoie de acuratețe moderată și simplitate, regula Simpson compusă e un implicit solid.
• Funcții cu frânturi sau discontinuități: Metodele cu puncte echidistante (trapeze, Simpson) tind să depășească cvadra­tura Gaussiană aici, deoarece plasarea fixă a nodurilor nu este perturbată de neregularități. Împarte intervalul de integrare la discontinuitățile cunoscute și tratează separat fiecare bucată.
• Date în puncte fixe, echidistante: cvadra­tura Gaussiană iese din discuție. Nu poți evalua funcția în noduri arbitrare, așa că folosește în schimb regula trapezelor compusă sau regula Simpson compusă.
• Integrale în dimensiuni mai mari: Cvadra­tura dă rapid peste blestemul dimensionalității. Cu n puncte pe dimensiune și d dimensiuni, ai nevoie de nᵈ evaluări în total. Zece dimensiuni cu zece puncte fiecare înseamnă 10¹⁰ evaluări, ceea ce nu se poate. Monte Carlo scalează mult mai bine peste aproximativ cinci dimensiuni, acolo se face trecerea.

Cvadra­tură vs. metode numerice înrudite

Cvadra­tura nu există în vid. Ajută să vezi unde se așază față de alte instrumente la care probabil te gândești.

Cvadra­tură vs. integrare analitică

Integrarea analitică găsește o primitivă exactă; cvadra­tura aproximează un număr. Dacă există o primitivă exactă și o poți găsi, folosește-o. Cvadra­tura își câștigă locul când integrarea exactă e imposibilă sau nepractică.

Cvadra­tură vs. integrare Monte Carlo

Ambele aproximează integrale definite, dar diferit. Cvadra­tura este deterministă: evaluează f în noduri fixe, prestabilite, cu convergență legată de netedimea funcției și de numărul de noduri. Monte Carlo extrage eșantioane aleatorii și le face media, cu eroare care scalează ca O(1/√N) indiferent de dimensiune. E o convergență mai lentă în dimensiuni mici, dar este independentă de dimensiune, ceea ce o face alegerea practică peste aproximativ cinci dimensiuni.

Cvadra­tură vs. interpolare

Unele reguli de cvadra­tură, inclusiv regula Simpson, funcționează potrivind un polinom prin punctele de evaluare și integrându-l exact. Dar scopul diferă de cel al interpolării. Interpolarea reconstruiește funcția; cvadra­tura vrea doar un număr. Polinomul e o unealtă, nu rezultatul.

Concluzie

Ce mi se pare satisfăcător la cvadra­tură este cât de mult câștigi dintr-o singură idee organizatoare. Fiecare metodă de aici, regula punctului mijlociu, trapezele, Simpson, Gaussiană, este de fapt un răspuns diferit la aceeași întrebare: având un buget de evaluări ale funcției, unde le pui și cum le ponderezi?

Potrivirea unor polinoame mai bune prin punctele de evaluare reduce rapid eroarea. Mutarea către noduri alese optim stoarce și mai multă acuratețe din același buget. Metodele adaptive apoi aduc totul laolaltă, direcționând evaluările exact acolo unde funcția are nevoie de ele.

În Python, scipy.integrate.quad() este implicitul potrivit pentru probleme 1D. Restul din acest articol este pentru momentele când vrei să înțelegi ce se întâmplă „sub capotă” sau când implementezi tu însuți ceva. Dacă vrei să-ți dezvolți în continuare abilitățile numerice cu SciPy, urmează cursul nostru Introduction to Optimization in Python.