Całkowanie numeryczne - metoda trapezów

powrót

Jedną z metod wyznaczania pola ograniczonego wykresem funkcji, osią OX oraz dwiema prostymi równoległymi do osi OY jest metoda trapezów. Nazwa tej metody wzięła się stąd, że dany obszar dzielimy na wiele trapezów prostokątnych. 

Załóżmy, że mamy funkcję $$y=f(x)$$, dwie proste $$x=a$$ i $$x=b$$, gdzie $$a<b$$:

Będziemy wyznaczać zacieniowane pole dzieląc ten obszar na mniejsze części. W omawianym przykładzie, dla czytelności rysunku podzielę ten obszar na cztery trapezy(w praktyce, aby otrzymać dokładne pole, należy go podzielić na kilka tysięcy takich czworokątów):

Następnie punkty przecięcia się tych podobszarów z wykresem funkcji połączymy ze sobą otrzymując cztery trapezy prostokątne: T1, T2, T3 oraz T4.

Suma pól tych trapezów będzie przybliżonym polem ograniczonego obszaru. Im większa liczba trapezów tyn dokładniejsze pole.

Pole pojedynczego trapezu liczymy oczywiście ze wzoru:

$$P=\frac{a+b}{2}\cdot h$$

Zacznijmy od wyznaczenia wysokości trapezu. Będzie ona równa:

$$h=\frac{b-a}{n}$$

gdzie $$n$$ to liczba trapezów, na które dzielmy omawiany obszar. Pozostaje nam jeszcze znalezienie długości podstaw, dla pierwszego trapezu będą one równe:

dłuższa_podstawa = $$|f(a)|$$

krótsza_podstawa = $$|f(a+h)|$$

Wartość funkcji może być ujemna więc musimy wyznaczyć jej wartość bezwzględną.

Pole i-tego trapezu liczymy ze wzoru:

$$P_i=\frac{|f(a+(i-1)\cdot h)|+|f(a+i\cdot h)|}{2}\cdot h$$

Sumując poszczególne pola trapezów otrzymujemy przybliżone szukane pole:

Reasumujmy:

Pole danego obszaru złożonego z $$n$$ trapezów będzie równe:

$$P=\frac{|f(a)|+|f(a+h)|}{2}\cdot h+\frac{|f(a+h)|+|f(a+2\cdot h)|}{2}\cdot h+\cdots +\frac{|f(a+(n-1)\cdot h)|+|f(a+n\cdot h)|}{2}\cdot h=$$

$$=h\cdot 0,5\cdot (|f(a)|+|f(a+h)|)+|f(a+h)|+|f(a+2\cdot h)|+\cdots +|f(a+(n-1)\cdot h)|+|f(a+n\cdot h)|)$$

Rozwiązanie w C++ dla następującej funkcji:

$$f(x)=x^2+x+2$$

//algorytm.edu.pl
#include<iostream>
#include<iomanip>
using namespace std;

double f(double x)
{
	//funkcja zawsze przyjmuje wartosci dodatnie
	//więc można pominąć wartosć bezwzględną
	return x*x+x+2;
}

double Pole(int a, int b, int n)
{
	double h = (b-a)/(double)n; //wysokosć trapezów
	double S = 0.0; //zmienna będzie przechowywać sumę pól trapezów
	double podstawa_a = f(a), podstawa_b;
		
	for(int i=1;i<=n;i++)
	{
		podstawa_b = f(a+h*i);
		S += (podstawa_a+podstawa_b);
		podstawa_a = podstawa_b;
	}
	return S*0.5*h;
}

int main()
{
	int a, b, n;
	cout<<"Podaj przedział [a, b]\na = ";
	cin>>a;
	cout<<"b = ";
	cin>>b;
	cout<<"Podaj liczbę trapezów: ";
	cin>>n;
	
	if(!(a<b))
		cout<<"To nie jest przedział!";
	else
		cout<<"Pole figury wynosi: "<<fixed<<setprecision(2)<<Pole(a, b, n);	

	return 0;
}