Vom calcula valoarea unei integrale definite

\( \int_a ^ b f (x) dx \)

Soluția descrisă aici a fost publicată într-una dintre disertațiile lui Thomas Simpson în 1743.

Formula lui Simpson

Fie \(n\) un număr natural oarecare. Împărțim segmentul de integrare \([a, b]\) în \(2n\) părți egale:

\(x_i = a + i h, ~~ i = 0 \ldots 2n,\)

\(h = \frac {b-a} {2n}.\)

Deci, să presupunem că luăm în considerare următorul segment \([x_ {2i-2}, x_ {2i}], i = 1 \ldots n\) . Înlocuim funcția \(f(x)\) pe el cu o parabolă \(P(x)\) care trece prin 3 puncte \((x_ {2i-2}, x_ {2i-1}, x_ {2i}, x_ {2i})\) . O astfel de parabolă există întotdeauna și este unică; ea poate fi găsită analitic. De exemplu, am putea să o construim folosind interpolarea polinomială Lagrange. Singurul lucru care rămâne de făcut este să integrăm acest polinom. Dacă faceți acest lucru pentru o funcție generală \(f\) , veți obține o expresie remarcabil de simplă:

\(\int_{x_ {2i-2}} ^ {x_ {2i}} f (x) ~dx \approx \int_{x_ {2i-2}} ^ {x_ {2i}} P (x) ~dx = \left(f(x_{2i-2}) + 4f(x_{2i-1})+(f(x_{2i})\right)\frac {h} {3}\)

Adăugând aceste valori pe toate segmentele, obținem formula finală a lui Simpson:

\(\int_a ^ b f (x) dx \approx \left(f (x_0) + 4 f (x_1) + 2 f (x_2) + 4f(x_3) + 2 f(x_4) + \ldots + 4 f(x_{2N-1}) + f(x_{2N}) \right)\frac {h} {3}\)

Eroare

Eroarea de aproximare a unei integrale prin formula lui Simpson este

\(-\tfrac{1}{90} \left(\tfrac{b-a}{2}\right)^5 f^{(4)}(\xi)\)

unde \(\xi\) este un număr oarecare între \(a\) și \(b\).

Eroarea este asimptotic proporțională cu \((b-a)^5\) . Cu toate acestea, derivatele de mai sus sugerează o eroare proporțională cu \((b-a)^4\) . Regula lui Simpson capătă un ordin în plus deoarece punctele în care se evaluează integrala sunt distribuite simetric în intervalul \([a, b]\)

Implementare

Aici, \(f(x)\) este o funcție definită de utilizator.

const int N = 1000 * 1000; // numărul de pași (deja înmulțit cu 2)
double simpson_integration(double a, double b){
    double h = (b - a) / N;
    double s = f(a) + f(b); // a = x_0 and b = x_2n
    for (int i = 1; i <= N - 1; ++i) { // Se referă la formula finală a lui Simpson
        double x = a + h * i;
        s += f(x) * ((i & 1) ? 4 : 2);
    }
    s *= h / 3;
    return s;
}