Całkowanie numeryczne: Metoda Simpsona

Metoda Simpsona jest najdokładniejszą z opisanych tutaj 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:

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\):
\[
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)\):

punkty podziałowe
dla i = 0,1,2,...,n \(y_i = f(x_i)\)
punkty środkowe
dla i = 1,2,...,n \(y_{t_i} = f(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})}\), \({(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 t_i^2 + b_i 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. 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 t_i^2 + b_i 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

Specyfikacja problemu

Dane wejściowe

\(x_p\)	- początek przedziału całkowania; \({x_p \in \mathbb{R}}\)
\(x_k\)	- koniec przedziału całkowania; \({x_k \in \mathbb{R}}\)
\(n\)	- liczba segmentów podziałowych; \({n \in \mathbb{N}}\)
\(f(x)\)	- funkcja rzeczywista, której całkę liczymy

Dane wyjściowe

\(s\)

- przybliżona wartość całki oznaczonej funkcji \(f(x)\) w przedziale \({[x_p ; x_k]}\); \({s \in \mathbb{R}}\)

Zmienne pomocnicze

\(s_t\)	- suma wartości funkcji w punktach środkowych; \({s_t \in \mathbb{R}}\)
\(h\)	- odległość na osi OX między dwoma sąsiednimi węzłami; \({h \in \mathbb{R}}\)
\(x\)	- współrzędna \(x\) węzła; \({x \in \mathbb{R}}\)
\(i\)	- licznik segmentów podziałowych; \({i \in \mathbb{N}}\)

Lista kroków

K01:	\(s \leftarrow 0; \; s_t \leftarrow 0\)
K02:	\(h = \frac{x_k - x_p}{n}\)
K03:	Dla \(i = 1,2,\dots,n\): wykonuj kroki K04...K06
K04:	\(x \leftarrow x_p + i \cdot h\)
K05:	\(s_t \leftarrow s_t + f(x - \frac{h}{2})\)
K06:	Jeśli \(i < n\), to \(s \leftarrow s + f(x)\)
K07:	\(s \leftarrow \frac{h}{6} \left(f(x_p) + f(x_k) + 2s + 4s_t \right)\)
K08:	Zakończ

Schemat blokowy

Odczytujemy krańce przedziału całkowania \({[x_p ; x_k]}\) oraz liczbę segmentów \(n\). Do obliczenia całki metodą Simpsona musimy zliczyć dwie sumy – wartości funkcji w węzłach \(x_i\) oraz wartości funkcji w punktach środkowych segmentów \(t_i\). Pierwszą sumę będziemy obliczać w zmiennej \(s\), drugą w \(s_t\). Obie na początku przyjmują wartość 0. Wyznaczamy dalej odległość pomiędzy dwoma sąsiednimi węzłami \(h\) i rozpoczynamy pętlę, w której zmienna \(i\) pełni rolę numeru segmentu. Pętla ta wykonuje się \(n\) razy od \({i = 1}\) do \({i = n}\) włącznie, czyli jest to zwykła pętla iteracyjna typu FOR.

W pętli wyznaczamy wartość węzła \(x_i\) i umieszczamy ją w zmiennej \(x\). Następnie obliczamy wartość funkcji w punkcie środkowym \(t_i\), który jest odległy o połowę \(h\) od wyznaczonego wcześniej punktu \(x_i\). Wartość tę dodajemy do sumy \(s_t\).

Drugą sumę tworzymy w zmiennej \(s\). Jednakże powinna ona zawierać jedynie wartości funkcji dla węzłów od \(x_1\) do \({x_{n - 1}}\). Dlatego przed sumowaniem sprawdzamy, czy indeks \(i\) jest w odpowiednim zakresie.

Po zakończeniu pętli wyznaczamy wartość całki w zmiennej \(s\) zgodnie z podanym wzorem, wyprowadzamy ten wynik dla użytkownika i kończymy wykonywanie algorytmu.

do podrozdziału do strony

W celu zobrazowania dokładności metody Simpsona zmniejszyliśmy w naszych przykładach liczb segmentów \({n = 10}\). Wynik obliczenia całki jest niezmieniony w stosunku do metod prostokątów i trapezów. Zachęcamy do eksperymentów z liczbą \(N\).

C++

// Obliczanie całki oznaczonej
// metodą Simpsona
// ---------------------------
//(C)2004 mgr Jerzy Wałaszek

#include <iomanip>
#include <iostream>
#include <cstdlib>

using namespace std;

// Tutaj definiujemy funkcję
double f(double x)
{
  return(x * x + 2 * x);
}

// Program główny
int main()
{
  //liczba segmentów
  const int N = 10;
  double xp,xk,s,st,dx,x;
  int i;

  // 3 cyfry po przecinku
  cout << setprecision(3)
  // format stałoprzecinkowy
       << fixed;

  cout << "Obliczanie calki oznaczonej\n"
          " za pomoca metody Simpsona\n"
          "---------------------------\n"
          "(C)2004 mgr J.Walaszek  I LO\n\n"
          "f(x) = x * x + 2 * x\n\n"
          "Poczatek przedzialu calkowania\n\n"
          "xp = ";
  cin >> xp;
  cout << "\nKoniec przedzialu calkowania\n\n"
          "xk = ";
  cin >> xk;
  cout << endl;
  s  = 0; st = 0;
  dx = (xk - xp) / N;
  for(i = 1; i <= N; i++)
  {
    x = xp + i * dx;
    st += f(x - dx / 2);
    if(i < N) s += f(x);
  }
  s = dx / 6 * (f(xp) + f(xk) + 2 * s + 4 * st);
  cout << "Wartosc calki wynosi : "
       << s
       << endl << endl;
  system("pause");
  return 0;
}

Pascal

// Obliczanie całki oznaczonej
// metodą Simpsona
// ---------------------------
// (C)2004 mgr Jerzy Wałaszek

program int_simpson;

// Tutaj definiujemy funkcję
function f(x : double) : double;
begin
  f := x * x + 2 * x;
end;

// Program główny

// liczba segmentów
const N = 10;

var
  xp,xk,s,st,dx,x : double;
  i : integer;

begin
  writeln('Obliczanie calki oznaczonej');
  writeln(' za pomoca metody Simpsona');
  writeln('---------------------------');
  writeln('(C)2004 mgr J.Walaszek I LO');
  writeln;
  writeln('f(x) = x * x + 2 * x');
  writeln;
  writeln('Poczatek przedzialu calkowania');
  writeln;
  write('xp = '); readln(xp);
  writeln;
  writeln('Koniec przedzialu calkowania');
  writeln;
  write('xk = '); readln(xk);
  writeln;
  s  := 0; st := 0;
  dx := (xk - xp) / N;
  for i := 1 to N do
  begin
    x  := xp + i * dx;
    st := st + f(x - dx / 2);
    if i < N then s := s + f(x);
  end;
  s := dx / 6 * (f(xp) + f(xk) + 2 * s + 4 * st);
  writeln('Wartosc calki wynosi : ',s:0:3);
  writeln;
  writeln('Nacisnij klawisz Enter...');
  readln;
end.

Basic

' Obliczanie całki oznaczonej
' metodą Simpsona 
' ---------------------------
' (C)2004 mgr Jerzy Wałaszek

Declare Function f(x As Double) _
                     As Double

' Program główny

' liczba segmentów
const N = 10

Dim As Double xp,xk,s,st,x,dx
Dim As Integer i

Print "Obliczanie  calki oznaczonej"
Print " za pomoca  metody Simpsona"
Print "----------------------------"
Print "(C)2004 mgr J.Walaszek  I LO"
Print
Print "f(x) = x * x + 2 * x"
Print
Print "Poczatek przedzialu calkowania"
Print
input "xp = ", xp
Print
Print "Koniec przedzialu calkowania"
Print
Input "xk = ", xk
Print
s  = 0: st = 0
dx = (xk - xp) / N
for i = 1 to N
  x  = xp + i * dx
  st += f(x - dx / 2)
  if i < N then s += f(x)
Next
s = dx / 6 * _
    (f(xp) + f(xk) + 2 * s + 4 * st)
print Using _
      "Wartosc calki wynosi : ####.###";s
Print
Print "Nacisnij klawisz Enter..."
Sleep
End

' Tutaj definiujemy funkcję
function f(x as Double) As Double
  f = x * x + 2 * x
End Function

Python (dodatek)

# Obliczanie całki oznaczonej
# metodą Simpsona 
# ---------------------------
# (C)2026 mgr Jerzy Wałaszek

# Tutaj definiujemy funkcję
def f(x):
    return x * x + 2 * x

# Program główny

# liczba segmentów
n = 10

print("Obliczanie  całki oznaczonej")
print(" za pomocą  metody Simpsona")
print("----------------------------")
print("(C)2026 mgr J.Wałaszek  I LO")
print()
print("f(x) = x * x + 2 * x")
print()
print("Początek przedziału całkowania")
print()
xp = float(input("xp = "))
print()
print("Koniec przedziału całkowania")
print()
xk = float(input("xk = "))
print()
s  = 0
st = 0
dx = (xk - xp) / n
for i in range(1, n + 1):
    x  = xp + i * dx
    st += f(x - dx / 2)
    if i < n: s += f(x)
s = (dx / 6 * 
    (f(xp) + f(xk) + 2 * s + 4 * st))
print(f"Wartość całki wynosi : {s:.3f}")
print()
input("Naciśnij Enter...")

Wynik:

Obliczanie  całki oznaczonej
 za pomocą  metody Simpsona
----------------------------
(C)2026 mgr J.Wałaszek  I LO

f(x) = x * x + 2 * x

Początek przedziału całkowania

xp = 0

Koniec przedziału całkowania

xk = 1

Wartość całki wynosi : 1.333

Naciśnij Enter...

JavaScript

<html>
<head>
  <title>Całkowanie numeryczne
  metodą Simpsona</title>
</head>
<body>

<div style="overflow-x: auto;"
     align="center">
  <table
  border="0"
  cellpadding="4"
  style="border-collapse:
         collapse">
    <tr>
      <td nowrap>
        <form
        name="frm"
        style="text-align: center;
               background-color:
               #E7E7DA">
          <b>
          Obliczanie całki
          oznaczonej
          <br>&nbsp;&nbsp;
          za pomocą
          metody Simpsona
          &nbsp;&nbsp;</b><br>
          (C)2026 mgr
          Jerzy Wałaszek
          <hr>
          Całkowana funkcja:<br>
          f(x) =
          x<sup>2</sup> + 2x
          <hr>
          Przedział całkowania<br>
          x<sub>p</sub> =
          <input
          name="xp_inp"
          size="15"
          value="0"
          type="text"
          style="text-align:
                 right">
          <br>
          x<sub>k</sub> =
          <input
          name="xk_inp"
          size="15"
          value="1"
          type="text"
          style="text-align:
                 right">
          <hr>
          <input
          onclick="js_simpson();"
          value="Oblicz całkę"
          name="B1"
          type="button">
          <hr>
          Wartość całki wynosi:
          <div id="out">.</div>
        </form>
      </td>
    </tr>
  </table>
</div>

<script language="JavaScript">

// Obliczanie całki oznaczonej
// metodą Simpsona
// ---------------------------
// (C)2004 mgr Jerzy Wałaszek

// Tutaj definiujemy funkcję
function f(x)
{
  return(x * x + 2 * x);
}

function js_simpson()
{
  // liczba segmentów
  var N = 10;
  var xp,xk,s,st,dx,x,i,t;

  xp = parseFloat(document.frm
       .xp_inp.value);
  xk = parseFloat(document.frm
       .xk_inp.value);
  if(isNaN(xp) || isNaN(xk))
    t = "<font color=red><b>" +
        "Popraw dane wejściowe!" +
        "</b></font>";
  else
  {
    s  = 0;
    st = 0;
    dx = (xk - xp) / N;
    for(i = 1; i <= N; i++)
    {
      x = xp + i * dx;
      st += f(x - dx / 2);
      if(i < N) s += f(x);
    };
    s = dx / 6 * (f(xp) + f(xk) +
        2 * s + 4 * st);
    t = Math.floor(s * 1000) /
        1000;
  };
  document.getElementById("out")
  .innerHTML = t;
}

</script>

</body>
</html>

do podrozdziału do strony

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.

Metoda Simpsona

Algorytm całkowania metodą Simpsona

Uwaga:

Opis algorytmu

Specyfikacja problemu

Dane wejściowe

Dane wyjściowe

Zmienne pomocnicze

Lista kroków

Schemat blokowy

Przykładowe programy