Całki

Wstęp

Całkę można intuicyjnie wytłumaczyć jako obszar ograniczony przez osie odcięte oraz wykres funkcji f(x) na pewnym przedziale [a, b]. Istnieje wiele rodzajów całek. Jednak nie zawsze można wyznaczyć całkę dla funkcji f(x). W przypadku obliczeń komputerowych nie zawsze istnieje potrzeba wyznaczania wzoru całki. W aplikacjach możemy otrzymać przybliżone wyniki dzięki Całkowaniu numerycznemu - znanemu też jako kwadratura.

Całkowanie Numeryczne

Rozpatrujemy całkę nieujemną funkcji f(x) na przedziale [a, b]. Całkę traktujemy tutaj jako pole powierzchni ograniczone przez proste x = a, x = b, oś odciętą oraz od góry poprzez wykres funkcji.

W Całkowaniu numerycznym, aby obliczyć pole powierzchni zastępujemy funkcję podcałkową prostszą funkcją, którą potrafimy policzyć. Jednak takie uproszczenie powoduje, że wyniki są tylko przybliżone. Istnieją jednak różne warianty pozwalające obliczyć pole tą metodą. Dzięki temu przybliżenia są na tyle małe, że można je zaniedbać.

Metoda trapezów

Przedział całkowania dzielimy na n równych odcinków. Następnie obliczamy wartości funkcji f(x) w punktach a i b oraz w punktach, które są końcami dwóch odcinków. Wtedy całkę obliczamy następująco:

W programie MATLAB istnieje gotowa funkcja trapz, która pozwala obliczyć całkę oznaczoną funkcji f na przedziale [a, b]. Funkcja trapz przyjmuje dwa argumenty. Pierwszy to wektor, który zawiera argumenty węzłów . Drugi argument to też wektor, który zawiera wartości funkcji w węzłach . Należy pamiętać, że oba wektory muszą być równej długości.

Przypuśćmy, że mamy funkcję postaci:

i chcemy obliczyć całkę na przedziale .

1. >> E = 1:0.1:3;
2. >> ZW = E.^3-4*E.^2-E+10;
3. >> calka = trapz(E, ZW)
4. 
   calka =
   
       1.3400

Na początek należy ustalić dziedzinę. Im więcej argumentów w danym przdziale tym dokładniejszy uzyskamy wynik. Następnie obliczamy wartość funkcji dla każdego argumentu z E. Na koniec obliczamy całkę przy uzyciu funkcji trapz podając dziedzinę E oraz zbiór wartości ZW. W ten sposób uzyskujemy, że .

Metoda prostokątów

Metoda prostokątów jest bardzo podobna do metody trapezów. Dzielimy dany przedział [a, b] na n odcinków. W każdym odcinku wyliczamy wartości funkcji dla argumentów z początku i końca odcinka. Uśredniony wynik to wysokość prostokąta o podstawie równej długości odcinka.

Metodę tę można wykorzystać w programie MATLAB korzystając z kilku dodatokwych funkcji: diff - do obliczania różnic pomiędzamy elementami listy, abs - liczbę, każdą wartość wektora / macierzy zamianiemy na wartość absolutną oraz sum - pozwala zsumować wszystkie wartości z wektora lub macierzy

1. >> E = 1:0.01:3;
2. >> ZW = E.^3-4*E.^2-E+10;
3. >> calka = sum(abs((ZW(1:end-1)+ZW(2:end))./2).*diff(E));

Metoda parabol (Simpsona)

Jest to metoda wykorzystywana w programie MATLAB do obliczania całek. Do obliczania całek tą metodą można użyć funkcji quad. Funkcja ta przyjmuje trzy wymagane argumenty. Kolejno jest to wzór funkcji, końce przedziału całkowania (jako dwie liczby). Argumentami opcjonalnymi są kolejno dokładność obliczeń oraz czy program ma pokazywać raport z przebiegu całkowania.

Obliczanie całek przy pomocy funkcji quad jest bardzo wygodne:

1. >> quad('x.^3-4*x.^2-x+10', 1, 3)
2. 
   ans =
   
       1.3333

Jak widać funkcja quad obliczyła z maksymalną możliwą dokładnością. Przy pomocy funkcji trapz też moglibyśmy zwiększyć dokładność obliczeń zmniejsząc odcinki pomiędzy kolejnymi argumentami. W napisie, który jest pierwszym argumentem funkcji quad, można wpisać nazwę funkcji. W ten sposób można wywołać np. funkcje trygonometryczne lub własne funkcje z m plików.

1. >> quad('sin(x.^3)', 0, pi)
2. 
   ans =
   
       0.4158

W przypadku, gdy potrzebujemy tylko dokładności do dwóch miejsc do przecinku to mozemy dodać dodatkowy argument i dać mu wartość 0.01. Zmiana dokładności obliczeń może prowadzić do znacznego odbiegania wyników od prawidłowych:

1. >> quad('sin(x.^3)', 1, 3, 0.01)
2. 
   ans =
   
       0.9537

Ostatni argument, który jest opcjonalny pozwala na podejrzenie kroków obliczeń całki. Wystarczy ustawić jego wartość na inną niż 0:

1. >> quad('sin(x.^3)', 0, pi, 0.1, 1)
2.     9     0.0000000000     8.53193733e-01     0.1289597680
      11     0.8531937329     1.43520519e+00     0.7746322569
      13     2.2883989207     8.53193733e-01     0.0500587780
   
   ans =
   
       0.9537