Serwis Edukacyjny
w I-LO w Tarnowie
obrazek

Materiały dla uczniów liceum

  Wyjście       Spis treści       Wstecz       Dalej  

obrazek

Autor artykułu: mgr Jerzy Wałaszek

©2026 mgr Jerzy Wałaszek

obrazek

Całkowanie

Metoda Simpsona

SPIS TREŚCI REMANENT
Podrozdziały
Wstęp
Algorytm
Przykładowa implementacja
 

Wstęp

Metoda Simpsona (ang. Simpson's Rule) jest klasyczną i jedną z najpopularniejszych numerycznych metod przybliżonego obliczania całek oznaczonych (kwadratur numerycznych). Chociaż nosi nazwisko angielskiego matematyka Thomasa Simpsona (1710 – 1761), jej historia jest doskonałym przykładem tego, jak ten sam koncept matematyczny był niezależnie odkrywany w różnych epokach i szerokościach geograficznych.

Antyczne korzenie (wkład Archimedesa i Keplera w stereometrię)

Archimedes położył geometryczny fundament pod metodę Simpsona, rozwiązując ten sam problem 2000 lat wcześniej za pomocą geometrii.

Związkiem między nimi jest parabola i obliczanie pola powierzchni pod jej łukiem. Metoda Simpsona opiera się na przybliżaniu dowolnej funkcji za pomocą łuków paraboli, a Archimedes jako pierwszy w historii precyzyjnie obliczył pole takiego obszaru.

Archimedes w swoim słynnym dziele Kwadratura paraboli udowodnił niezwykłe twierdzenie geometryczne.

Wykazał, że pole odcinka paraboli (obszaru zamkniętego przez parabolę i prostą) wynosi dokładnie \(\frac{4}{3}\) pola trójkąta wpisanego w ten odcinek, mającego tę samą podstawę i wysokość.

Do udowodnienia tego użył tzw. metody wyczerpywania (antycznego odpowiednika rachunku całkowego), wpisując w parabolę nieskończony ciąg coraz mniejszych trójkątów i sumując ich pola jako szereg geometryczny.

Jeśli weźmiemy geometryczne twierdzenie Archimedesa i zapiszemy je we współczesnym języku algebry oraz układu współrzędnych kartezjańskich, otrzymamy dokładnie ten sam wzór, który Thomas Simpson opublikował w 1743 roku dla pojedynczego przedziału.

W metodzie Simpsona pole pod parabolą rozpiętą na trzech punktach \(x_0\), \(x_1\), \(x_2\) oblicza się ze wzoru:

\[\int\limits_{x_{0}}^{x_{2}}f(x)dx\approx \frac{h}{3}[f(x_{0})+4f(x_{1})+f(x_{2})]\]

Jeśli podstawisz do tego wzoru czystą funkcję kwadratową, wynik będzie w 100% identyczny z proporcją \(\frac{4}{3}\) odkrytą przez Archimedesa. Simpson nie musiał wyprowadzać nowej geometrii pola paraboli – ona była znana od starożytności dzięki Archimedesowi. Simpson przeniósł ją na poziom analityczny i połączył w ciąg wielu przedziałów.

Łącznikiem między Archimedesem a Simpsonem był Johannes Kepler (1571–1630). W 1615 roku Kepler szukał sposobu na szybkie liczenie objętości beczek. Znał dzieła Archimedesa, wziął jego geometrię paraboli i uprościł ją do praktycznego wzoru (reguły beczki), który potrzebował tylko trzech pomiarów: dna beczki, jej góry i najszerszego punktu w środku. Simpson ponad wiek później ubrał tę regułę Keplera i geometrię Archimedesa w nowoczesny calculus.

Bez matematycznego geniuszu Archimedesa, który jako pierwszy ujarzmił pole powierzchni krzywej parabolicznej, metoda Simpsona nie miałaby swojego geometrycznego rdzenia.

Prawdziwy odkrywca: Bonaventura Cavalieri (1635 r.)

Pierwszym matematykiem, który formalnie sformułował i opublikował geometryczny odpowiednik tej metody, był włoski geometra Bonaventura Cavalieri (1598–1647) – uczeń Galileusza. W swoim dziele o geometrii niepodzielnych (1635 r.) opisał regułę, która odpowiada dzisiejszemu całkowaniu numerycznemu za pomocą wielomianów drugiego stopnia.

James Gregory i pierwsza formalna notacja (1668 r.)

Większość współczesnych historyków matematyki zgadza się, że autorem analitycznej wersji tego wzoru jest szkocki matematyk i astronom James Gregory (1638–1675). W 1668 roku opublikował on pracę Exercitationes Geometricae, w której podał ogólne wzory na przybliżanie pól pod krzywą za pomocą łuków paraboli. Simpson urodził się dopiero 42 lata po publikacji Gregory'ego.

Thomas Simpson i popularyzacja metody (1743 r.)

Angielski matematyk-samouk Thomas Simpson (1710 – 1761) opublikował ten wzór w swoim podręczniku Mathematical Dissertations on a Variety of Physical and Analytical Subjects w 1743 roku.

Co właściwie zrobił Simpson? Nie odkrył tej metody od zera. Wziął znane już w kręgach akademickich założenia (współpracował m.in. z pracami Newtona i Rogera Cotesa nad wzorami Newtona-Cotesa) i opisał je w sposób niezwykle klarowny, przystępny i łatwy do zastosowania w ówczesnej inżynierii oraz nawigacji morskiej. Podręcznik Simpsona stał się bestsellerem, przez co nazwisko autora na stałe przylgnęło do tej kwadratury.

Niemiecki ślad: Wzór Keplera-Simpsona

W niemieckiej literaturze matematycznej i technicznej metoda ta do dziś bywa nazywana regułą Simpsona-Keplera lub Keplersche Fassregel, co oddaje sprawiedliwość historyczną pierwotnym obserwacjom Keplera.

Matematyczna kwintesencja algorytmu

Zamiast przybliżać pole pod krzywą za pomocą prostych odcinków (jak w prostszej metodzie trapezów), metoda Simpsona dzieli przedział całkowania na parzystą liczbę podprzedziałów i nad każdą parą rozpina łuk paraboli (wielomian stopnia 2).

Dzięki temu, przy stosunkowo niewielkim nakładzie obliczeniowym (zliczanie wartości funkcji w punktach węzłowych i odpowiednie ich ważenie współczynnikami 1, 4, 2, 4..., 1), metoda ta cechuje się bardzo wysokim rzędem dokładności – jest dokładna dla wszystkich wielomianów stopnia do trzeciego włącznie.


do podrozdziału  do strony 

Algorytm

Metoda Simpsona (ang. Simpson's Rule) jest jedną z dokładniejszych metod przybliżonego całkowania. W metodzie prostokątów całka oznaczona przybliżana była funkcjami stałymi – liczyliśmy sumę pól prostokątów. W metodzie trapezów całkę przybliżaliśmy za pomocą funkcji liniowych - obliczaliśmy sumy pól trapezów. W metodzie Simpsona stosujemy jako przybliżenie parabolę – będziemy obliczali sumy wycinków obszarów pod parabolą. Zasada jest następująca:

W przedziale całkowania \({[x_p ; x_k]}\) wyznaczamy \({n + 1}\) równoodległych punktów \({x_0, x_1, x_2 ,\dots, x_n}\), które obejmują cały przedział całkowania. Najpierw liczymy odległość pomiędzy dwoma kolejnymi punktami (ponieważ punkty są równoodległe, odległość ta jest stała):

\[
h = \frac{x_k - x_p}{n}
\]

Teraz łatwo policzymy położenie każdego punktu \(x_i\):

dla \(i = 0,1,2,\dots,n\):
\[
x_i = x_p + i \cdot h
\]

Następnie wyliczamy punkt leżący w środku podprzedziału \({[x_{i - 1},x_i]}\)

dla \(i = 1,2,\dots,n\):
\[
x_{t_i} = x_{i-1} + \frac{h}{2} = x_i - \frac{h}{2}
\]

Dla każdego z wyznaczonych w ten sposób punktów obliczamy wartość funkcji \(f(x)\):

współrzędne \(y\)
punktów podziałowych
dla i = 0,1,2,...,n
\(y_i = f(x_i)\)
współrzędne \(y\)
punktów środkowych
dla i = 1,2,...,n
\(y_{t_i} = f(x_{t_i})\)

W każdym podprzedziale \({[x_{i - 1} ; x_i]}\) przybliżamy funkcję za pomocą paraboli \(g(x)\) o następującej postaci:

dla \(i = 1,2,\dots,n\):
\[
g_i(x) = a_i \cdot x^2 + b_i \cdot x + c_i; \; x \in [x_{i-1},x_i], \; \text{dla } i = 1,2,\dots,n
\]

Parabola \(g_i(x)\) musi przechodzić przez punkty: \({(x_{i-1}, y_{i - 1})}\), \({(x_{t_i}, y_{t_i})}\), \({(x_i, y_i)}\). Współczynniki \(a_i\), \(b_i\) i \(c_i\) wyznaczymy zatem z układu trzech równań:

dla \(i = 1,2,\dots,n\):
\[
\begin{cases}
a_i x_{i-1}^2 + b_i x_{i-1} + c_i = y_{i-1} \\
a_i x_{t_i}^2 + b_i x_{t_i} + c_i = y_{t_i} & \text{dla } i = 1,2,\dots,n \\
a_i x_i^2 + b_i x_i + c_i = y_i
\end{cases}
\]

Uwaga:

W metodzie Simpsona chodzi o wyznaczenia pola pod parabolą w danym podprzedziale, a nie jej współczynników dla tego podprzedziału. Możemy zatem pójść inną drogą. Załóżmy, iż powyższe współczynniki są znane (ostatecznie możemy je przecież wyliczyć).

Pole pod parabolą w przedziale \({[x_{i-1} ; x_i]}\) będzie równe całce oznaczonej:

dla \(i = 1,2,\dots,n\):
\[
P_i = \int\limits_{x_{i-1}}^{x_i} g_i(x) \, dx = \int\limits_{x_{i-1}}^{x_i} \left(a_i x^2 + b_i x + c_i\right) dx; \quad x \in [x_{i-1}, x_i]
\]

Funkcja pierwotna jest bardzo prosta w tym przypadku i ma wzór następujący:

dla \(i = 1,2,\dots,n\):
\[
G_i(x) = \int g_i(x) \, dx = \frac{a_i}{3}x^3 + \frac{b_i}{2}x^2 + c_i x + C
\]

Wartość całki obliczymy zgodnie z definicją Newtona-Leibniza:

dla \(i = 1,2,\dots,n\):
\[
\begin{aligned}
\int\limits_{x_{i-1}}^{x_i} g_i(x) \, dx &= G_i(x_i) - G_i(x_{i-1}) \\
\int\limits_{x_{i-1}}^{x_i} g_i(x) \, dx &= \frac{a_i}{3}x_i^3 + \frac{b_i}{2}x_i^2 + c_i x_i - \frac{a_i}{3}x_{i-1}^3 - \frac{b_i}{2}x_{i-1}^2 - c_i x_{i-1} \\
\int\limits_{x_{i-1}}^{x_i} g_i(x) \, dx &= \frac{a_i}{3}x_i^3 - \frac{a_i}{3}x_{i-1}^3 + \frac{b_i}{2}x_i^2 - \frac{b_i}{2}x_{i-1}^2 + c_i x_i - c_i x_{i-1} \\
\int\limits_{x_{i-1}}^{x_i} g_i(x) \, dx &= \frac{a_i}{3}\left(x_i^3 - x_{i-1}^3\right) + \frac{b_i}{2}\left(x_i^2 - x_{i-1}^2\right) + c_i\left(x_i - x_{i-1}\right)
\end{aligned}
\]

Teraz postaramy się uprościć maksymalnie otrzymane wyrażenie. W tym celu wyciągamy przed nawias wspólny czynnik i całość dzielimy przez 6:

dla \(i = 1,2,\dots,n\):
\[
\begin{aligned}
\int\limits_{x_{i-1}}^{x_i} g_i(x) \, dx &= \frac{(x_i - x_{i-1})}{6} \left\{ 2a_i\left(x_i^2 + x_i x_{i-1} + x_{i-1}^2\right) + 3b_i\left(x_i + x_{i-1}\right) + 6c_i \right\} \\
\int\limits_{x_{i-1}}^{x_i} g_i(x) \, dx &= \frac{(x_i - x_{i-1})}{6} \left( 2a_i x_i^2 + 2a_i x_i x_{i-1} + 2a_i x_{i-1}^2 + 3b_i x_i + 3b_i x_{i-1} + 6c_i \right) \\
\int\limits_{x_{i-1}}^{x_i} g_i(x) \, dx &= \frac{(x_i - x_{i-1})}{6} \left\{ \left(a_i x_{i-1}^2 + b_i x_{i-1} + c_i\right) + \left(a_i x_i^2 + b_i x_i + c_i\right) + a_i\left(x_{i-1} + x_i\right)^2 + 2b_i\left(x_{i-1} + x_i\right) + 4c_i \right\} \\
\int\limits_{x_{i-1}}^{x_i} g_i(x) \, dx &= \frac{(x_i - x_{i-1})}{6} \left\{ \left(a_i x_{i-1}^2 + b_i x_{i-1} + c_i\right) + \left(a_i x_i^2 + b_i x_i + c_i\right) + 4 \left( a_i \left( \frac{x_{i-1} + x_i}{2} \right)^2 + b_i \frac{x_{i-1} + x_i}{2} + c_i \right) \right\} \\
\int\limits_{x_{i-1}}^{x_i} g_i(x) \, dx &= \frac{(x_i - x_{i-1})}{6} \left\{ \left(a_i x_{i-1}^2 + b_i x_{i-1} + c_i\right) + \left(a_i x_i^2 + b_i x_i + c_i\right) + 4\left(a_i x_{t_i}^2 + b_i x_{t_i} + c_i\right) \right\}
\end{aligned}
\]

Zwróćcie uwagę, iż wyrażenia w nawiasach są odpowiednio wartościami funkcji \({y_{i -1}}\), \(y_i\) oraz \(y_{t_i}\). Natomiast różnica

dla \(i = 1,2,\dots,n\):
\[
x_i - x_{i - 1}
\]

jest odległością \(h\) pomiędzy dwoma sąsiednimi punktami podziałowymi. Zatem po uproszczeniu otrzymujemy ostateczny wzór:

dla \(i = 1,2,\dots,n\):
\[
\int\limits_{x_{i-1}}^{x_i} g_i(x) \, dx = \frac{h}{6} \left(y_{i-1} + y_i + 4y_{t_i}\right)
\]

Wzór ten pozwala wyliczyć pole obszaru pod parabolą aproksymującą funkcję \(f(x)\) w przedziale \({[x_{i-1} ; x_i]}\). Wartość całej całki otrzymamy sumując te pola, czyli:

\[
\int\limits_{x_p}^{x_k} f(x) \, dx \approx \frac{h}{6} \left( \sum_{i=1}^{n} y_{i-1} + \sum_{i=1}^{n} y_i + 4 \sum_{i=1}^{n} y_{t_i} \right)
\]

Jest to wzór wyliczania przybliżonej wartości całki oznaczonej za pomocą metody Simpsona. Ponieważ w obliczanych sumach wartości funkcji w węzłach się powtarzają dwukrotnie na krańcach stykających się przedziałów (z wyjątkiem pierwszej i ostatniej wartości, ponieważ te węzły nie nie są wspólne dla przedziałów pierwszego i ostatniego n), do obliczeń komputerowych stosujemy efektywniejszy wzór otrzymywania powyższej sumy:

\[
\begin{aligned}
&\sum_{i=1}^{n} y_{i-1} + \sum_{i=1}^{n} y_i = y_0 + y_n + 2 \sum_{i=1}^{n-1} y_i \\
&\int\limits_{x_p}^{x_k} f(x) \, dx \approx \frac{h}{6} \left( y_0 + y_n + 2 \sum_{i=1}^{n-1} y_i + 4 \sum_{i=1}^{n} y_{t_i} \right)
\end{aligned}
\]

do podrozdziału  do strony 

Przykładowa implementacja

Program liczy całkę oznaczoną:

\[
\int\limits_0^{\pi} \left(\sin(x) + cos(x)\right) \, dx
\]

Całka ta ma wartość dokładną (analityczną) 2. Program liczy całkę dla 100 segmentów (101 punktów) metodą Simpsona następnie wyświetla wynik. Wyjaśnienia są w komentarzach

C++
// Całkowanie metodą Simpsona
// (C)2026 mgr Jerzy Wałaszek
// Metody numeryczne 0083
//---------------------------

#define _USE_MATH_DEFINES
#include <iostream>
#include <windows.h>
#include <iomanip>
#include <cmath>

using namespace std;

// Definicje stałych
//------------------

// Liczba segmentów
const int N = 100;
// Przedział całkowania
const double Xp = 0;
const double Xk = M_PI;
// Szerokość segmentów
const double H = (Xk - Xp) / N;

// Całkowana funkcja
double f(double x)
{
  return sin(x) + cos(x);
}

// Metoda Simpsona
double s_r()
{
  int i;
  double x, s, s1, s2;

  // Punkt startowy środków
  x = Xp + H / 2;
  // Wartość całki
  s = f(Xp) + f(Xk);
  s1 = 0;
  s2 = f(x);
  for(i = 1; i < N; i++)
  {
    s1 += f(Xp + i * H);
    s2 += f(x + i * H);
  }
  s1 *= 2;
  s2 *= 4;
  s += s1 + s2;
  s *= H / 6;
  return s;
}

// Program główny
//---------------
int main()
{
  SetConsoleOutputCP(CP_UTF8);
  SetConsoleCP(CP_UTF8);

  cout << setprecision(7)
       << fixed;

  cout
  << "Obliczanie całki oznaczonej "
     "metodą Simpsona\n\n"
  << "wynik : " << setw(9)
  << s_r() << endl << endl;

  system("pause");
  return 0;
}
Wynik: 
Obliczanie całki oznaczonej metodą Simpsona

wynik : 2.0000000

Zwróć uwagę na dokładność wyniku – specjalnie zwiększyłem precyzję do 7 cyfr po przecinku. Poeksperymentuj z tym programem, np. zmniejszaj liczbę segmentów i porównuj wyniki, wyciągając wnioski.

Zastanów się nad usprawnieniem pętli, np. poprzez zastąpienie mnożenia \({i * H}\) dodawaniem stałej \(H\) w każdym obiegu. Kiedy mnożenie ma przewagę nad dodawaniem (rozważ kumulacje błędów przy dodawaniu)?

Python (dodatek) 
# Całkowanie metodą Simpsona
# (C)2026 mgr Jerzy Wałaszek
# Metody numeryczne 0083
#---------------------------

from math import sin, cos, pi

# Definicje stałych
#------------------

# Liczba segmentów
n = 100
# Przedział całkowania
xp = 0
xk = pi
# Szerokość segmentów
h = (xk - xp) / n

# Całkowana funkcja
def f(x):
    return sin(x) + cos(x)

# Metoda Simpsona
def s_r():
    # Punkt startowy środków
    x = xp + h / 2
    # Wartość całki
    s = f(xp) + f(xk)
    s1 = 0
    s2 = f(x)
    for i in range(1,n):
        s1 += f(xp + i * h)
        s2 += f(x  + i * h)
    s1 *= 2
    s2 *= 4
    s += s1 + s2
    s *= h / 6
    return s

# Program główny
#---------------

print(f"Obliczanie całki oznaczonej "
      f"metodą Simpsona\n\n"
      f"wynik : {s_r():9.7f}\n")

input("Naciśnij Enter...")

do podrozdziału  do strony 
Zespół Przedmiotowy
Chemii-Fizyki-Informatyki
w I Liceum Ogólnokształcącym
im. Kazimierza Brodzińskiego
w Tarnowie
ul. Piłsudskiego 4
©2026 mgr Jerzy Wałaszek

Materiały tylko do użytku dydaktycznego. Ich kopiowanie i powielanie jest dozwolone pod warunkiem podania źródła oraz niepobierania za to pieniędzy.
Pytania proszę przesyłać na adres email: i-lo@eduinf.waw.pl
Serwis wykorzystuje pliki cookies. Jeśli nie chcesz ich otrzymywać, zablokuj je w swojej przeglądarce.

Informacje dodatkowe.