Lab. 6 Całkowanie numeryczne

Całka oznaczona

$$I=\int_{a}^{b} f(x) d x$$

pole powierzchni pod funkcją na odcinku $[a,b]$

Metody

• kwadratura - przybliżenie całkowanej funkcji za pomocą $Q(f)=\sum_{i=0}^{n} A_{i} f\left(x_{i}\right)$
• kwadratury Newtona-Cotesa, o stałych odległościach węzłów
• kwadratury Gaussa
• całkowanie Monte Carlo

Reguła prostokątów

$$\int_{a}^{b} f(x) dx \approx (b-a) f\left(\frac{a+b}{2}\right)$$

Reguła trapezów

$$\int_{a}^{b} f(x) d x \approx(b-a)\left(\frac{f(a)+f(b)}{2}\right)$$

Reguła parabol (Simpsona)

$$\int_{a}^{b} f(x) d x \approx \frac{b-a}{6}\left[f(a)+4 f\left(\frac{a+b}{2}\right)+f(b)\right]$$

Złożone kwadratury

Dzielimy przedział $[a,b]$ na $n$ odcinków z krokiem $h = \frac{b-a}{n}$

$$ x_0 =a, \qquad x_i = a + ih, \qquad x_n = b$$

wartość całki będzie sumą kwadratur prostych z każdego odcinka $[x_i, x_{i+1}]$

• złożona reguła prostokątów (midpoint)
  $$I = h \sum_{i=0}^{n-1} f\left(\frac{x_i + x_{i+1}}{2}\right) \qquad \text{błąd} \quad \frac{b-a}{24} h^{2} f^{(2)}(\xi)$$
• złożona reguła trapezów
  $$I = h \left(\frac{f(x_0)+ f(x_n)}{2}+\sum_{i=1}^{n-1} f(x_i) \right) \qquad \text{błąd} \quad \frac{b-a}{12} h^{2} f^{(2)}(\xi)$$
• złożona reguła Simpsona
  $$I = \frac{h}{3} \sum_{j=1}^{n / 2}\left ( f(x_{2 j-2})+4 f(x_{2 j-1})+f(x_{2 j})\right) \qquad \text{błąd} \quad \frac{b-a}{180} h^{4} f^{(4)}(\xi)$$

Przykład

$$\int_{0}^{2} f(x) d x$$

Dla kwadratur prostych ($n=1$)

midpoint.c

#include <stdio.h>
#include <math.h>

double f1(double x)
{
   return x * x; 
}

double midpoint_rule(double a, double b, int n, double (*func)(double))
{
   double s, fa, fb, x, h;
   int j;

   h = (b-a)/((double) n);
   s=0.0;

   for (j=0; j < n; j++)
   {
      x = a + (j+0.5) * h;
      s += (*func)(x);
   }
   s *= h;

   return s;
} 

int main()
{
   int n;
   double a=0.0, b=2.0;
   
   printf("Number of partitions = ");
   scanf("%d", &n) ;
   
   printf("f1  = %lf\n", midpoint_rule(a, b, n, f1)) ;
   
   return 0;
}

Zadanie

Znając wartość poniższej całki $$\int_{0}^{1} \frac{1}{1+x^{2}} dx = \frac{\pi}{4}$$ napisz program, który wyznaczy przybliżoną wartość $\pi$ za pomocą całkowania numerycznego. Wynik wyznacz za pomocą 3 metod: złożonej metody prostokątów, trapezów i Simpsona. Ilość węzłów dla każdej metody dobierz automatycznie poprzez iteracyjne zwiększanie liczby przedziałów $n =2, 4, 8, 16, \ldots $ aż do sytuacji, gdy wartość przybliżona całki dla kolejnych wartości $n$ będzie mniejsza od założonej dokładności $|I_{n-1}-I_n|<\epsilon$. Dla każdej z metod program wypisuje przybliżoną wartość $\pi$ oraz ilość podziałów $n$ potrzebnych do osiągnięcia żądanej precyzji. Obliczenia wykonaj dla $\epsilon=10^{-6}$.

Dodatkowe ćwiczenia

Ćw. 1

Wyznacz całki z przykładu z zajęć za pomocą metody złożonej trapezów oraz Simpsona i porównaj dokładność tych metod.

Ćw. 2

Napisz program wyznaczający przybliżoną wartość całki $\int_0^{2.5} f(x) dx$ funkcji o wartościach zawartych w tabeli