%% Całkowanie numeryczne
% Opracował dr inż. Robert Jakubowski, Politechnika Rzeszowska, Katedra Samolotów i Silników Lotniczych
% Table of Contents
% Wprowadzenie
% Metoda prostokątów
% Metoda trapezów
%     Metoda trapezów dla danych punktowych
% Całka Simpsona
%     Adaptacja metody Simpsona do obliczeń dla danych punktowych
% Całkowanie z wykorzystaniem funkcji wewnętrznych Matlab
% Funkcje przygotowane do wykładu
%% Wprowadzenie
% Całkowanie numeryczne polega na przybliżonym wyznaczeniu wartości całki oznaczonej
% S=\int_{\textrm{x1}}^{\textrm{x2}} f\left(x\right)\textrm{dx}
% Przedstawione wcześniej metody realizowały obliczenia w oparciu o Symbolic Toolbox. Tym razem przyjrzymy się metodom, które wykorzystuje się w obliczeniach numerycznych.
% Zakładając, że mamy doczynienia z całką właściwą, czyli ciągłą w przedziale całkowania i o skończonych granicach, to obliczenia wykonuje się z wykorzystaniem kwadratur w ostaci:
% S=\int_{\textrm{x1}}^{\textrm{x2}} f\left(x\right)\textrm{dx}\approx \sum_{i=1}^n F_i
% Gdzie F_i jest przyblizoną wartością pola pod funkcją podcałkową dla elementarnego przedziału \Delta x_i na jaki został podzielony przedział całkowania.
%% Metoda prostokątów
%
% Rys 1. Graficzne zobrazowanie metody prostokątów
% Metoda prostokątów przedstawiona graficznie na rys. 1 polega na podzieleniu całkowanego przedziału na n odcinków i wyliczeniu sumy pól utworzonych w ten sposób prostokątów. Prostokąty można utworzyć na trzy sposoby w zależńości od metody wyboru węzłów. Może to być metoda (rys. 2)
% lewych prostokątów gdy węzeł do utworzenia prostokąta jest położony z lewej strony
% prawych prostokątów, gdy wężeł leży z prawej strony
% środkowych prostokątów gdy węzeł leży pośrodku; x_m =\frac{1}{2}\left(x_1 +x_2 \right)
%
% Rys 2. Metody tworzenia prostokątów a) metoda lewych prostokątów, b) metoda prawych prostokątów, c) metoda środkowych prostokątów
% Dla pierwszej metody funkcja całkowania metodą prostokątów została przedstawiona w funkcji calka_prostok_l
% Spróbujmy wykorzystać przygotowaną funkcję i policzyć wartość całki. Policzmy całkę z funkcji sin(x) od 0 do \piczyli:
% S=\int_0^{\pi } \sin \left(x\right)\textrm{dx}

f=@(x) sin(x); %definicja funkcji
p_c=[0,pi]; %przedział całkowanie
n=10; % ilość przedziałów na które podzielony zostanie przedział całkowania
S1=calka_prostok_l(f,p_c,n)

%
% Sprawdzmy wynik obliczeń realizując wykorzystując wcześniej poznaną metodę obliczeń symbolicznych

syms z
fz=sin(z)
S=double(int(fz,z,0,pi))

% Wynik obliczenia całki dla metody dokładnej dał wartość 2, która jest nieco wyzsza, niż otrzymana metodą prostokątów
%
% Spróbujmy obliczyć wartość całki gdy weźmiemy pod uwagę metodę całkowania z wykorzystaniem środkowych prostokątnych. Metodę taką przedstawiono w funkcji calka_prostok_srod. Weżmy pod uwagę te same parametry całkowania, co poprzednio

S2=calka_prostok_srod(f,p_c,n)

% Tym razem wyznaczona wartość całki jest większa niż dla rozwiązania dokładnego
%
% W ramach ćwiczeń proszę przygotować funkcję całka_prostok_p, która będzie wyliczać całkę metodą prawych prostokątów w polu do którego odsyła link. Następnie poniżej należy wyznaczyć wartość całki tą metodą



%
% Sprawdzmy jak zwiększenie ilości węzłów całkowania wpłynie na dokładność obliczeń. Zwiększmy n do 50

n=50;
S1_1=calka_prostok_l(f,p_c,n)
S2_1=calka_prostok_srod(f,p_c,n)

% Otrzymane wyniki pokazują, że jesteśmy coraz bliżej rozwiązania dokładnego. Gdyby dalej zagęszczać przedział całkowania to wyliczana wartość będzie dokładniejsza
%
%% Metoda trapezów
% Kolejną metodą, która daje wynik obarczony mniejszym błędem jest metoda trapezów. Polega ona na przybliżeniu obszaru całkowania w trapezami jak przedstawiono to na rysunku 3.
%
%
% Rys 3 Przykład całkowania metodą trapezów
% W przedstawionym przypadku pole pojedynczego przedziału bedzie obliczane zgodnie z definicją pola trapezu:
% F_1 =0,5\left(f\left(x_{12} \right)+f\left(x_{11} \right)\right)\left(x_{12} -x_{11} \right)
% Funkcję liczącą całkę metodą trapezów przedstawiono w funkcja calka_trapez

n=10;
St=calka_trapez(f,p_c,n)
St_1=calka_trapez(f,p_c,50)

% W przypadku obliczeń całki z funkcji sinus dostajemy tak samo dokłdny wynik jak metodą lewych prostokątów. Oczywiście dla wiekszości innych funkcji metoda ta da znacznie lepsze wyniki obliczeń
%
%% Metoda trapezów dla danych punktowych
% Metodę trapezów można wykorzystać do obliczeń całki z wartości dyskretnych. Na przykład, gdy dla zadanych wartości x mamy podane wartości y i oczywiście dane są posortowane wg rosnących wartości x. Dla takiego przypadku funkcja do obliczeń całki metodą trapezów przedstawiono w funkcji calka_trapez_xy
% Przygotujmy dane w postaci punktów x i y

x=linspace(-1,1.5,21);
y=x.^3+2*x.^2+2;
%punkty przedstawiono graficznie 
plot(x,y,'p')
axis([x(1),x(end),0,10])
S_t_xy=calka_trapez_xy(x,y)

% Ponieważ do otrzymania punktów wykorzystaliśmy funkcję, więc zweryfikujmy wynik obliczeń wykonując obliczenia na zmiennych symbolicznych

syms z
fz=z.^3+2*z.^2+2;
S_d=double(int(fz,z,-1,1.5))

% Widzimy, że otrzymany wynik dla obliczeń dokładnych nieznacznie rózni się od wyniku uzyskanego metodą całkowania trapezowego
%
%% Całka Simpsona
% W tym przypadku pole elementarnej figury F jest obliczone jako obszar pod krzywą drugiego stopnia. Do tego, aby znaleźć funkcję drugiego stopnia niezbędne są trzy punkty, dlatego też w wybranym przedzialeod x11 do x12 dla wyznaczenia elementarnego pola wprowadza się dodatkowo x1m=0.5*(x11+x12). Na bazie tych trzech punktów można opisać wielomian drugiego stopnia, przybliżający całkowaną funkcję (rys 4). Teorię dotyczącą tworzenia takiej funkcji zawarto w pracy [1]. Finalnie pole po scałkowaniu funkcji pod funkcją kwadratową jest wyrażane
%
% F_1 =\frac{x_{12} -x_{11} }{6}\left(f\left(x_{11} \right)+4f\left(x_{1m} \right)+f\left(x_{12} \right)\right)
% gdzie x_{1m} =\frac{1}{2}\left(x_{11} +x_{12} \right)
% Oczywiście punkt końcowy tego przedziału x12 jest punktem początkowym kolejnego przedziały czyli x21=x12 itd. Stąd funkcja całkowania metodą Simpsona może być przedstawiona jako
% S=\sum_{i=1}^n \frac{x_{i2} -x_{i1} }{6}\left(f\left(x_{i1} \right)+4f\left(0\ldotp 5\left(x_{i1} +x_{i2} \right)\right)+f\left(x_{i2} \right)\right)
% gdzie n - jest liczbą przedziałów na który został podzielony obszar całkowania
% Metodę tę wykorzystano w funkcji calka_simpsona, Wykorzystując tą funkcję spróbujmy wyznaczyć wartość całki funkcji sin jak robiliśmy to dla wcześniejszych metod i porównajmy uzyskane wyniki

n=10;
Ss=calka_simpsona(f,p_c,n)
n=50;
Ss_1=calka_simpsona(f,p_c,n)
tab1=table([S; S], [10; 50], [S1; S1_1], [S2; S2_1], [St; St_1], [Ss; Ss_1]);
tab1.Properties.VariableNames={'Wynik dokładny' 'liczba przedziałów' 'm. lewych prost.'...
    'm. środkowych prost.' 'm. trapezów' 'm. Simpsona'};
tab1.Properties.Description=' Porównanie wyników całkowania różnymi metodami';
tab1

% Tab. 1. Porównanie wyników całkowania różnymi metodami
% Wyniki przedstawione w tab.1 pokazują, że zwiększając ilość przedziałów całkowania podnosi się dokładność obliczeń. W przedztsawionych przykładach już przy przyjętej liczbie przedziałów 10 metoda Simpsona daje rozwiązanie bliskie rozwiązaniu dokładnemu. W kolejności otrzymanych wyników można stwierdzić, że metoda Simpsona pozwala wyznaczyć wartość całki z największą dokładnością, później jest metoda trapezów i najmniej dokładną jest metoda prostokątów.
%
%% Adaptacja metody Simpsona do obliczeń dla danych punktowych
% Metode Simpsona podobnie jak metodę trapezów można dostosować do całkowania na bazie punktów. Przyjrzyjmy się rysunkowi 4. Do określenia krzywej drugiego stopnia dla pojedynczego segmentu potrzebne jest trzy punkty, zatem dla pierwszego segmentu musimy mieć taką ilość punktów. Każdy kolejny segment wymaga dołożenia kolejnych dwóch punktów, bowiem pierwszy punkt w tym segmencie jest ostatnim w segmencie poprzedzającym. Wynika stąd, że dla przyjętego toku rozumowania niezbedne jest aby policzyć całkę Simpsona minimum trzy punkty, oraz drubim wymaganiem jest aby liczba punktów była nieparzysta. Wtedy można wydzielić (n-1)/2 segmentów, na bazie których zostanie wykonane całkowanie metodą Simpsona.
%
% Zależność na całkę Simpsona dla zbioru punktów można zapisać
% S=\sum_{i=1}^{\frac{\left(n-1\right)}{2}} \frac{\left(x\left(2i+1\right)-x\left(2i-1\right)\right)}{6}\left(f\left(x\left(2i-1\right)\right)+4f\left(x\left(2i\right)\right)+f\left(x\left(2i+1\right)\right)\right)
% Zależność tę przedstawiono w funkcji calka_simpsona_xy
% Spróbujmy wykorzystać tą metodę do obliczenia całki na podstawie wcześniej wyznaczonych punktów x i y  (przy metodzie trapezów). Celow wprowadziłem tam nieparzystą liczbę punktów. Porównajmy otrzymane wyniki

S_s_xy=calka_simpsona_xy(x,y)
tab2=table({'Metoda dokładna'; 'Metoda trapezów'; 'Metoda Sipsona'},[S_d;S_t_xy;S_s_xy], [0;abs(S_t_xy-S_d);abs(S_s_xy-S_d)]);
tab2.Properties.VariableNames={'Metoda całkowania' 'Wynik całkowania' 'błąd względny'};
tab2

% Tab. 2. Porównanie wyników całkowania dla matody trapezów i metody Simpsona na zbioerze punktów
% Zestawienie wyników pokazuje, że metoda Simpsona w tym wypadku też daje rozwiązanie bardzo bliskie rozwiązania dokładnego, Metoda trapezów w tym wględzie daje nieznaczny błąd.
%
%% Całkowanie z wykorzystaniem funkcji wewnętrznych Matlab
% Można realizować całkowanie wykorzystując standardowo występujące w matlabie funkcje do całkowania integral, integral2, integral3. Funkcje te pozwalają na wyznaczanie wartości całki oznaczonej jednej dwóch i trzech zmiennych. Wywołanie funkcji integral1 jest następujące S=integral(fun,pocz_prz,kon_prz),
% gdzie: fun -  jest funkcją, pocz_prz - początek przedziału całkowania, kon_prz - koniec przedziału całkowania.
% Spróbujmy wykonać całkowanie na funkcji wcześniej rozpatrywanej F(x)=x.^3+2*x.^2+2 w przedziale jak wcześniej od -1 do 1.5.

f=@(x) x.^3+2*x.^2+2
S=integral(f,-1,1.5)

% Pozostałe funkcje mają podobną strukturę wywołania funkcji.
% funkcja integer2 jest wywoływana S=integral2(f(x,y),x_min,x_max,y_min,y_max).
% Zróbmy obliczenia dla wyznaczania połowy pola koła, które zrobiono dla całkowania z wykorzystaniem symbolic toolbox. Wykorzystamy funkcję f\left(r,\alpha \right)=\sqrt{r^2 {\cos \left(\alpha \right)}^2 +r^2 {\sin \left(\alpha \right)}^2 }

fk=@(r,alpha) ((r.*cos(alpha)).^2+(r.*sin(alpha)).^2).^0.5

%  Policzmy wartość całki dla r1=0 do r2=5 i kąta alpha1=0 i alpha2=pi, co jest równowarzne wyliczeniu pola koła o promieniu R=5. Ponieważ pierwszą zmienną w fynkcji jest r zatem jako pierwsze będą granice całkowania po r, a następnie po alpha

S_k05=integral2(fk,0,5,0,2*pi)
%Sprawdzenie wyniku z obliczeń dla koła
S_k05_spr=pi*5^2

%
% W ramach ćwiczeń można policzyć objętość kuli wykorzystując całkę potrójną przedstawioną integral3, wywołanej jako V=integral2(f(x,y,z),x_min,x_max,y_min,y_max,z_min,z_max). Do zależności na objetość kuli można wykorzystać funkcję przedstawioną w materiałach Działania w Matlabie na wyrażeniach symbolicznych cz. 2




%
%
%% Funkcje przygotowane do wykładu
% Funkcja calka_prostok_l
% Zmiennymi wejścowymi będę funkcja f, przedział całkowania [x1, x2] i ilość elemntów na których dzielony jest przedział. Zmienną wyjściową jest wartość całki S

function S=calka_prostok_l(f,x,n)
if length(x)~=2 
   %Komunikat o niewłaściwym przedziale zmiennej x
    Disp "Przedział całkowania powinien być wektorem dwuelementowym [x1, x2]"
    % Przypisanie wyniku jako NotaNumber
    S=NaN;
else
    xn=linspace(x(1),x(2),n+1);
    S=0;
    % Iteracyjne wyliczenie wartości całki poprzez sumowanie pól
    for i=1:n
        S=S+f(xn(i))*(xn(i+1)-xn(i));
    end
end    
end

% Funkcja calka_prostok_srod
% Zmiennymi wejścowymi będę funkcja f, przedział całkowania [x1, x2] i ilość elemntów na których dzielony jest przedział. Zmienną wyjściową jest wartość całki S

function S=calka_prostok_srod(f,x,n)
if length(x)~=2 
   %Komunikat o niewłaściwym przedziale zmiennej x
    Disp "Przedział całkowania powinien być wektorem dwuelementowym [x1, x2]"
    % Przypisanie wyniku jako NotaNumber
    S=NaN;
else
    xn=linspace(x(1),x(2),n+1);
    S=0;
    % Iteracyjne wyliczenie wartości całki poprzez sumowanie pól przy
    % założeniu że węzeł będzie punkem śtodkowym
    for i=1:n
        S=S+f(0.5*(xn(i)+xn(i+1)))*(xn(i+1)-xn(i));
    end
end    
end

% Funkcja całka_prostok_p
% W podobny sposób jak wyżej należy przygotować funkcję do całkowania metodą prawych prostokątów. Można wykorzystać te same zmienne wejściowe i wyjściowe, jak powyżej



% Funkcja calka_trapez
% Zmiennymi wejścowymi będę funkcja f, przedział całkowania [x1, x2] i ilość elemntów na których dzielony jest przedział. Zmienną wyjściową jest wartość całki S

function S=calka_trapez(f,x,n)
if length(x)~=2 
   %Komunikat o niewłaściwym przedziale zmiennej x
    Disp "Przedział całkowania powinien być wektorem dwuelementowym [x1, x2]"
    % Przypisanie wyniku jako NotaNumber
    S=NaN;
else
    xn=linspace(x(1),x(2),n+1);
    S=0;
    % Iteracyjne wyliczenie wartości całki poprzez sumowanie pól trapezów 
    for i=1:n
        S=S+(xn(i+1)-xn(i))*0.5*(f(xn(i))+f(xn(i+1)));
    end
end    
end

% Funkcja calka_trapez_xy
% Funkcja oblicza wartość całki metodą trapezów dla dyskretnego zbioru danych posortowanego wg roznących wartości x. Zmiennymi wejściowymi są wartości x i odpowiadające im wartości y

function S=calka_trapez_xy(x,y)
n=length(x);
S=0;
for i=1:n-1
    S=S+0.5*(y(i+1)+y(i))*(x(i+1)-x(i));
end
end

% Funkcja calka_simpsona
% Zmiennymi wejścowymi będę funkcja f, przedział całkowania [x1, x2] i ilość elemntów na których dzielony jest przedział. Zmienną wyjściową jest wartość całki S

function S=calka_simpsona(f,x,n)
if length(x)~=2 
   %Komunikat o niewłaściwym przedziale zmiennej x
    Disp "Przedział całkowania powinien być wektorem dwuelementowym [x1, x2]"
    % Przypisanie wyniku jako NotaNumber
    S=NaN;
else
    xn=linspace(x(1),x(2),n+1);
    S=0;
    % Iteracyjne wyliczenie wartości całki poprzez sumowanie pól trapezów 
    for i=1:n
        xm=0.5*(xn(i+1)+xn(i));
        S=S+(xn(i+1)-xn(i))/6*(f(xn(i))+4*f(xm)+f(xn(i+1)));
    end
end    
end

% Funkcja calkasimpsona_xy
% Funkcja oblicza wartość całki metodą Simpsona dla dyskretnego zbioru danych posortowanego wg roznących wartości x. Zmiennymi wejściowymi są wartości x i odpowiadające im wartości y

function S=calka_simpsona_xy(x,y)
n=length(x);
n1=length(y);
if n1~=n
    % zabezpieczenie dla źle wprowadzonych danych wektor x i y mają różne
    % długości
    S=NaN;
    disp 'wektor wejściowy x i y muszą mieć tę samą długość'
elseif n<3
    %Zabezpieczenie gdy liczba elementów jest zbyt mała
    S=NaN;
    disp 'Nie można policzyć całki Simpsona dla mniej niż trzech punktów'
elseif round(n/2)==n/2
    %Zabezpieczenie gdy liczba elementów jest parzysta
    S=NaN;
    disp 'Wymagane jest do obliczeń, aby liczba punktów była nieparzysta'
else
    % Obliczenie całki
    S=0;
    for i=1:(n-1)/2
        S=S+(x(2*i+1)-x(2*i-1))/6*(y(2*i-1)+4*y(2*i)+y(2*i+1));
    end
end
end

% Literatura
% 1) Pańczyk B. i in.; Metody Numeryczne w przykładach, Politechnika Lubelska, Lublin 2012
% 2) Fortuna Z., Macukow B., Wąsowski J.; Metody Numeryczne; Wydawnictwa Naukowo-Techniczne, Warszawa 2001