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:

Przedział całkowania <x_p, x_k> dzielimy na n + 1 równo odległych punktów x_o, x₁, x₂,..., x_n:

dla i = 0,1,2,...,n
obrazek

Dla każdych dwóch sąsiednich punktów wyznaczamy punkt środkowy t_i wg wzoru:

dla i = 1,2,...,n

Obliczamy odległość między dwoma sąsiednimi punktami.

Dla każdego wyznaczonego w ten sposób punktu obliczamy wartość funkcji f(x) w tym punkcie:

punkty podziałowe
dla i = 0,1,2,...,n f_i = f(x_i)
punkty środkowe
dla i = 1,2,...,n f_ti = 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,...,n

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

dla i = 1,2,...,n

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,...,n

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

dla i = 1,2,...,n

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

dla i = 1,2,...,n

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,...,n

Zwróćcie uwagę, iż wyrażenia w nawiasach są odpowiednio wartościami funkcji f_{i -1}, f_i oraz f_ti. Natomiast różnica

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

dla i = 1,2,...,n

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:

Jest to wzór wyliczania przybliżonej wartości całki oznaczonej za pomocą metody Simpsona. Ponieważ w obliczanych sumach wartości funkcji się powtarzają dwukrotnie (z wyjątkiem pierwszej i ostatniej), do obliczeń komputerowych stosujemy efektywniejszy wzór otrzymywania powyższej sumy:

do podrozdziału do strony

Specyfikacja problemu

Dane wejściowe

x_p	- początek przedziału całkowania, x_p ∈ R
x_k	- koniec przedziału całkowania, x_k ∈ R
n	- liczba punktów podziałowych, n ∈ 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 ∈ R

Zmienne pomocnicze

s_t	- suma wartości funkcji w punktach środkowych, s_t ∈ R
dx	- odległość między dwoma sąsiednimi punktami podziałowymi, dx ∈ R
x	- pozycja punktu podziałowego, x ∈ R
i	- licznik punktów podziałowych, i ∈ N

Lista kroków

K01:	s ← 0; s_t ← 0
K02:
K03:	Dla i = 1,2,...,n: wykonuj kroki K04...K06
K04:	x ← x_p + i ·dx
K05:
K06:	Jeśli i < n, to s ← s + f(x)
K07:
K08:	Zakończ

Schemat blokowy

Odczytujemy krańce przedziału całkowania <x_p,x_k>. Do obliczenia całki metodą Simpsona musimy zliczyć dwie sumy - wartości funkcji w punktach podziałowych x_i oraz wartości funkcji w punktach środkowych przedziałó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 punktami podziałowymi dx i rozpoczynamy pętlę, w której zmienna i pełni rolę numeru punktu podziałowego i środkowego. 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ść punktu podziałowego 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ę dx 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 punktów podziałowych od x₁ 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ę punktów podziałowych do 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 punktów podziałowych
  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 punktów podziałowych
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 punktów podziałowych
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 punktów podziałowych
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("Wartość całki wynosi : ",end="")
print("%.3f" % s)
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 punktów
  // podziałowych
  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

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