Całkowanie numeryczne - metoda prostokątów

powrót

jednym ze sposobów wyznaczenia pola obszaru ograniczonego wykresem funkcji, osią OX oraz dwiema prostymi równoległymi do osi OY jest podzielenie tego obszaru na prostokąty:

W naszym przykładzie podzielmy dany obszar na trzy prostokąty, które będą przecinać się z wykresem funkcji dokładnie po środku każdego boku. Oczywiście im większa liczba takich prostokątów tym wyznaczone pole jest dokładniejsze.

Wzór na pole prostokąta jest wszystkim znany:

$$P=x\cdot y$$

Długość jednego z boków wyznaczymy ze wzoru:

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

gdzie, $$n$$ to liczba prostokątów na jakie dzielimy szukane pole, natomiast $$b$$ i $$a$$ to proste, które ograniczają zadany obszar $$a<b$$.

Drugi bok wyznaczymy wyliczając wartość funkcji w środku boku prostokąta:

$$y=|f(srodek_{i})|$$

Środek pierwszego boku wyznaczymy ze wzoru:

$$srodek=a+\frac{b-a}{2n}$$

Każdy następny będziemy wyznaczać zwiększając poprzedni o długość boku $$x$$. 

Ostatecznie sumujemy pola wszystkich prostokątów otrzymując przybliżone pole zakreślanego obszaru:

$$P=x\cdot |f(srodek)|+x\cdot |f(srodek+x)|+x\cdot |f(srodek+2x)|+\cdots +x\cdot |f(srodek+(n-1)\cdot x)|=$$

$$x\cdot (|f(srodek|+|f(srodek+x|+|f(srodek+2x|+\cdots |f(srodek+(n-1)\cdot x)|$$

Rozwiazanie w C++:

//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 x = (b-a)/(double)n; //pierwszy bok - każdy prostokąt ma taki sam
  double S = 0.0; //zmienna będzie przechowywać sumę pól trapezów
  double srodek = a+(b-a)/(2.0*n); //środek pierwszego boku
 
  for(int i=0;i<n;i++)
  {
      S+=f(srodek); //obliczenie wysokości prostokąta
      srodek+=x; //przejście do następnego środka    
  }
  return S*x;
}
 
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;
}