Całkowanie numeryczne: Metoda trapezów

Opisana w poprzednim rozdziale metoda prostokątów nie jest zbyt dokładna, ponieważ pola użytych w niej prostokątów źle odwzorowują powierzchnię pola pod krzywą (dokładność odwzorowania rośnie wraz ze wzrostem liczby prostokątów). Dużo lepszym rozwiązaniem jest zastosowanie zamiast nich trapezów o wysokości \(h\) i podstawach równych odpowiednio wartości funkcji w punktach krańcowych. Sama zasada nie zmienia się.

W przedziale całkowania \({[x_p ; x_k]}\) wyznaczamy \({n+1}\) równoodległych punktów \({x_0,x_1,\dots,x_n}\). Punkty te wyznaczamy w prosty sposób wg wzoru:

\[
x_i = x_p + \frac{i}{n} \cdot (x_k - x_p) \quad \text{dla } i = 0,1,2,\dots,n
\]

Obliczamy odległość między dwoma sąsiednimi punktami - będzie to wysokość każdego trapezu:

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

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

\[
y_i = f(x_i), \; \text{dla} \; i = 0,1,2,\dots,n
\]

Pole pod wykresem funkcji przybliżane jest polami \(n\) trapezów. Pole i-tego trapezu obliczamy wg wzoru:

\[
P_i = \frac{y_{i-1} + y_i}{2} \cdot h, \; \text{dla} \; i = 1,2,\dots,n
\]

Przybliżona wartość całki jest sumą pól wszystkich otrzymanych w ten sposób trapezów:

\[
s = P_1 + P_2 + \dots + P_n
\]

czyli

\[
\begin{aligned}
s &= \frac{y_0 + y_1}{2} \cdot h + \frac{y_1 + y_2}{2} \cdot h + \frac{y_2 + y_3}{2} \cdot h + \dots + \frac{y_{n-2} + y_{n-1}}{2} \cdot h + \frac{y_{n-1} + y_n}{2} \cdot h \\
s &= \frac{h}{2} \cdot (y_0 + y_1 + y_1 + y_2 + y_2 + y_3 + \dots + y_{n-2} + y_{n-1} + y_{n-1} + y_n) \\
s &= \frac{h}{2} \cdot (y_0 + 2y_1 + 2y_2 + \dots + 2y_{n-1} + y_n) \\
s &= h \cdot \left(y_1 + y_2 + \dots + y_{n-1} + \frac{y_0 + y_n}{2}\right)
\end{aligned}
\]

Wyprowadzony na końcu wzór jest podstawą przybliżonego wyliczania całki w metodzie trapezów.

\(
\int\limits_{x_p}^{x_k} f(x) \, dx \approx h \cdot \left( \sum_{i=1}^{n-1} f(x_p + i \cdot h) + \dfrac{f(x_p) + f(x_k)}{2} \right)
\)

Przykład:

Obliczymy ręcznie przybliżoną wartość całki oznaczonej z funkcji \({f(x) = \sin(x)}\)w przedziale \({[0 ; \pi]}\). W przedziale wyznaczamy \({n + 1 = 5}\) punktów, \({n = 4}\) – dzieje się tak dlatego, iż zaczynamy od punktu \(x_0\), a kończymy na \(x_n\), co daje \({n + 1}\) punktów, które podzielą przedział na dokładnie \(n\) równych segmentów:

\[
\begin{aligned}
x_0 &= 0 \\
x_1 &= \frac{1}{4}\pi \\
x_2 &= \frac{2}{4}\pi \\
x_3 &= \frac{3}{4}\pi \\
x_4 &= \pi
\end{aligned}
\]

Odległość między dwoma sąsiednimi punktami (czyli długość segmentu) wynosi:

\[
h = \frac{x_k - x_p}{n} = \frac{\pi}{4} \approx 0{,}7854\dots
\]

Dla każdego z wyznaczonych punktów \(x_i\) obliczamy wartość funkcji \({y_i = f(x_i) = \sin(x_i)}\):

\[
\begin{array}{lcl}
y_0 = f(x_0) = \sin 0 & = & 0{,}0000 \\
y_1 = f(x_1) = \sin\frac{1}{4}\pi & \approx & 0{,}7071 \\
y_2 = f(x_2) = \sin\frac{2}{4}\pi & = & 1{,}0000 \\
y_3 = f(x_3) = \sin\frac{3}{4}\pi & \approx & 0{,}7071 \\
y_4 = f(x_4) = \sin\pi & = & 0{,}0000
\end{array}
\]

Obliczamy sumę pól trapezów:

\[
\begin{aligned}
s &= h \cdot \left( y_1 + y_2 + y_3 + \frac{y_0 + y_4}{2} \right) \\
s &\approx 0{,}7854 \cdot (0{,}7071 + 1 + 0{,}7071) = 1{,}8961
\end{aligned}
\]

Otrzymaliśmy wynik identyczny jak w przypadku metody prostokątów - dlaczego?

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; \({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

h	- długość segmentu między dwoma sąsiednimi punktami podziałowymi; \({h \in \mathbb{R}}\)
\(i\)	- licznik segmentów; \({i \in \mathbb{N}}\)

Lista kroków

K01:	\(s \leftarrow 0\)
K02:	\(h \leftarrow \dfrac{x_k - x_p}{n}\)
K03:	Dla \(i = 1,2,\dots,n-1\) : wykonuj: \(s \leftarrow s + f \left(x_p + i \cdot h \right)\)
K04:	\(s \leftarrow \left(s + \dfrac{f(x_p) + f(x_k)}{2} \right) \cdot h\)
K05:	Zakończ

Schemat blokowy

Po odczytaniu informacji o krańcach \(x_p\) i \(x_k\) oraz liczbie segmentów \(n\) przedziału całkowania ustawiamy sumę \(s\) na 0, obliczamy odległość \(h\) pomiędzy sąsiednimi punktami podziałowymi i ustawiamy licznik segmentów na 1.

Rozpoczynamy pętlę iteracyjną, która wykona się \({(n - 1)}\) razy. Wewnątrz pętli obliczamy i-ty punkt podziałowy oraz wartość funkcji w tym punkcie, którą dodajemy do sumy \(s\).

Po zakończeniu pętli do sumy \(s\) dodajemy średnią wartość funkcji na krańcach przedziału i całość przemnażamy przez \(h\). Po tej operacji \(s\) zawiera wartość przybliżoną całki. Zwracamy ją użytkownikowi i kończymy algorytm.

do podrozdziału do strony

Prezentowane poniżej programy wyliczają całkę oznaczoną funkcji \({f(x) = x^2 + 2x}\). Przedział całkowania jest dzielony na \({n = 1000}\) segmentów. W przedziale \({[0 ; 1]}\) całka oznaczona ma wartość dokładną równą \({\frac{4}{3} = 1{,}3333\dots}\).

C++

// Obliczanie całki oznaczonej
// metodą trapezów
// ---------------------------
// (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 trapezów podziałowych
  const int N = 1000;

  double xp,xk,s,h;
  int i;

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

  cout << "Obliczanie calki oznaczonej\n"
          " za pomoca metody trapezow\n"
          "---------------------------\n"
          "(C)2004 mgr J.Walaszek I LO\n\n"
          "f(x) = x * x + 2 * x\n\n"
          "Podaj poczatek przedzialu "
          "calkowania\n\n"
          "xp = ";
  cin >> xp;
  cout << "\nPodaj koniec przedzialu "
          "calkowania\n\n"
          "xk = ";
  cin >> xk;
  cout << endl;
  s  = 0;
  h = (xk - xp) / N;
  for(i = 1; i < N; i++)
    s += f(xp + i * h);
  s = (s + (f(xp) + f(xk)) / 2) * h;
  cout << "Wartosc calki wynosi : "
       << s
       << endl << endl;
  system("pause");
  return 0;
}

Pascal

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

program int_trapez;

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

// Program główny

// liczba trapezów podziałowych
const N = 1000;

var
  xp,xk,s,h : double;
  i : integer;

begin
  writeln('Obliczanie calki oznaczonej');
  writeln(' za pomoca metody trapezow');
  writeln('--------------------------');
  writeln('(C)2004 mgr J.Walaszek I LO');
  writeln;
  writeln('f(x) = x * x + 2 * x');
  writeln;
  writeln('Podaj poczatek ' +
          'przedzialu calkowania');
  writeln;
  write('xp = '); readln(xp);
  writeln;
  writeln('Podaj koniec ' +
          'przedzialu calkowania');
  writeln;
  write('xk = '); readln(xk);
  writeln;
  s  := 0;
  h := (xk - xp) / N;
  for i := 1 to N - 1 do
    s := s + f(xp + i * h);
  s := (s + (f(xp) + f(xk)) / 2)* h;
  writeln('Wartosc calki wynosi : ',
          s:0:3);
  writeln;
  writeln('Nacisnij klawisz Enter...');
  readln;
end.

Basic

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

Declare Function f(x As Double) _
                     As Double

' Program główny

' liczba trapezów podziałowych
const N = 1000

Dim As Double xp,xk,s,h
Dim As Integer i

Print "Obliczanie  calki oznaczonej"
Print " za pomoca  metody trapezow"
Print "----------------------------"
Print "(C)2004 mgr J.Walaszek  I LO"
Print
Print "f(x) = x * x + 2 * x"
Print
Print "Podaj poczatek " & _
      "przedzialu calkowania"
Print
Input "xp = ", xp
Print
Print "Podaj koniec " & _
      "przedzialu calkowania"
Print
Input "xk = ", xk
Print
s  = 0
h = (xk - xp) / N
for i = 1 to N - 1
  s += f(xp + i * h)
Next
s = (s + (f(xp) + f(xk)) / 2) * h
Print Using "Wartosc calki " & _
      "wynosi : ####.###"; s

' Gotowe
Print
Print "KONIEC. Nacisnij " & _
      "dowolny klawisz..."
Sleep
End

' Tutaj definiujemy funkcję
Public Function f(ByVal x As Double) _
                  As Double
    Return x * x + 2 * x
End Function

Python (dodatek)

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

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

# Program główny

# liczba trapezów podziałowych
n = 1000

print("Obliczanie całki oznaczonej")
print(" za pomocą metody trapezów")
print("---------------------------")
print("(C)2026 mgr J.Wałaszek I LO")
print()
print("f(x) = x * x + 2 * x")
print()
print("Podaj początek "
      "przedziału całkowania")
print()
xp = float(input("xp = "))
print()
print("Podaj koniec "
      "przedziału całkowania")
print()
xk = float(input("xk = "))
print()
s  = 0
h = (xk - xp) / n
for i in range(1,n):
    s += f(xp + i * h)
s = (s + (f(xp) + f(xk)) / 2) * h
print(f"Wartość całki "
      f"wynosi : {s:.3f}")

# Gotowe
print()
input("Naciśnij Enter...")

Wynik:

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

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

Podaj początek przedziału całkowania

xp = 0

Podaj koniec przedziału całkowania

xk = 1

Wartość całki wynosi : 1.333

Naciśnij Enter...

JavaScript

<html>
<head>
  <title>Całkowanie numeryczne
  metodą trapezów</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 trapezów
          &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_trapez();"
          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ą trapezów
// ---------------------------
// (C)2004 mgr Jerzy Wałaszek

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

function js_trapez()
{
  // liczba trapezów
  // podziałowych
  var N = 1000;
  var xp,xk,s,h,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;
    h = (xk - xp) / N;
    for(i = 1; i < N; i++)
      s += f(xp + i * h);
    s = (s + (f(xp) + f(xk)) /
         2) * h;
    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 trapezów

Algorytm całkowania metodą trapezów

Opis algorytmu

Specyfikacja problemu

Dane wejściowe

Dane wyjściowe

Zmienne pomocnicze

Lista kroków

Schemat blokowy

Przykładowe programy