\documentclass[a4paper,11pt] {article}
\usepackage[utf8]{inputenc}
\usepackage[czech]{babel}
\usepackage[left=2cm,top=4cm,bottom=3cm]{geometry}
\usepackage{url}
\author{Radek Žemlička}
\title{Úloha č. 4: Výpočet integrálu -- dokumentace}

\begin{document}
\maketitle
	\section{Zadání}
		Napište program, který vypočítá integrál z nějaké funkce Simpsonovým pravidlem. Otestujte ho na takové nekonstantní funkci a takovém intervalu, kdy známý analytický výsledek bude a) nula, b) jedna. Diskutujte vliv počtu dělení intervalu na numerickou chybu.

	\section{Řešení}
		\subsection{Popis metody}
			\subsubsection{Numerické integrování obecně }
				Pro numerický výpočet určitého integrálu 
			
				\begin{equation}	
					\int_{a}^{b}{f(x) dx}
					\label{integral}
				\end{equation}

				použijeme možnost nahradit integrovanou funkci $f$ interpolačním polynomem.  \\
				Při definování určitého (Riemanova) integrálu se postupuje tak, že se interval $ \langle a,b \rangle $  rozdělí dělícími body $a = x_0 < x_1 < x_2 < ... < x_n = b$ .
				Normou dělení D rozumíme číslo \cite{nummet1} 
				\begin{equation}	
					 h = {max \atop i = 0,1,...,n-1} (x_{i+1} - x_i) 
					\label{norma}
				\end{equation}

				Pro výpočet tedy rozdělíme interval $ \langle a,b \rangle $ na dílčí intervaly (pro jednoduchost) stejné délky, a to uzlovými body $a = x_0 < x_1 < x_2 < ... < x_n = b$, kde $x_i = x_0 + i \cdot h, i = 0,1,...,n$ a kde $h = (b-a)/n$ je konstantní krok.

			
			\subsubsection{Simpsonovo pravidlo}
				Nechť počet stejně dlouhých dílčích intervalů, na které rozdělíme integrační interval $ \langle a,b \rangle $ integrálu $ \int_{a}^{b}{f(x) dx}$ je sudý.  Za těchto podmínek platí takzvané {\bf Simpsonovo pravidlo}: \cite{nummet1, nummet2}
				\begin{equation}
					\int_{a}^{b}{f(x)dx} \approx \frac{1}{3} h [f(x_0) + 4f(x_1) + 2f(x_2) + 4f(x_3)+...+ 4f(x_{n-1}) + f(x_n)] = S(f,h)
					\label{homer}
				\end{equation}

		\subsection{Vlastní program}
			Program na výpočet vlastního integrálu pomocí Simpsonova pravidla je psán v C++ (ačkoliv nevyužívá žádných funkcí, nebo vlastností, jejichž plnohodnotnou obdobu by C nemělo)  a je tak krátký, že je zde možné rozebrat celý jeho zdrojový kód.\\
			Hlavička programu obsahuje pouze volání knihovny {\tt iostream} (aby mohl být proveden výstup například na obrazovku) a knihovny {\tt cmath} (aby uživatel programu mohl použít funkce knihovny při zadávání  matematické funkce, kterou chce zintegrovat). \\
			Program obsahuje pouze tři funkce: {\tt double f}, {\tt double simpson} a {\tt int main}, a jednu pomocnou funkci {\tt double analytic}, která provede kontrolní výpočet pomocí přesného analytického tvaru integrálu z integrované funkce. 

			- Funkce {\tt main} nejprve deklaruje tři proměnné. Do těchto proměnných jsou uloženy velikosti horní a dolní meze Riemanova integrálu a počet dělících intervalů. (Uživatel může vložit hodnoty do těchto proměnných již ve zdrojovém kódu, nebo až prostřednictvím příkazové řádky po spuštění programu. Viz komentáře ve zdrojovém kódu k zjištění, jak tuto volbu provést.)  Funkce dále deklaruje promněnnou {\tt result}, do které vloží takovou hodnotu, jakou vrátí funkce {\tt simpson}. Parametry funkce {\tt simpson} jsou právě ony tři zmíněné proměnné. Funkce {\tt main} pak pošle hodnotu proměnné {\tt result} na standartní výstup. Pokud se uživatel rozhodne ke kontrole výsledku porovnáním s analytickým řešením, využije ještě následující dva řádky kódu. Pokud ne, musí řádky zakomentovat. První řádek deklaruje proměnnou {\tt double analyticresult}, do které vloží rozdíl, do analytického tvaru integrálu dosazených, mezních hodnot. Dosazení se docílí voláním funkce {\tt analytic}. Druhý řádek pak posílá analytický výsledek na standartní výstup. Nakonec funkce {\tt main}  program ukončí.  

			- Funkce {\tt simpson} prakticky provádí samotnou integraci Simpsonovou metodou, a proto ji zde vypíšeme celou:

			\begin{verbatim}
double simpson(double n, double a, double b)
{
        double h = (b - a)/n;
        double S = f(a) + 4 * f(a+h) +f(b);
        for(int i = 2; i != n; i += 2){
                S = S + ( 2 * f(a + i*h) );
                S = S + ( 4 * f(a + (i+1)*h) );
        }
        return (1.0/3.0)*h*S;
}
			\end{verbatim}

			Asi nepřekvapí, že proměnná {\tt a} obsahuje hodnotu dolní meze,   {\tt b} horní meze  a {\tt n} počet dělících intervalů. Funkce nejprve spočítá velikost dělícího intervalu {\tt h}, který potřebuje k dalším výpočtům a poté deklaruje proměnnou {\tt S}, do které uloží první, druhý a poslední člen řady ze vztahu $\ref{homer}$, aby následující for-cyklus mohl co nejefektivněji přičíst do {\tt S} členy zbylé. Funkce {\tt simpson} nakonec vrátí funkci {\tt main} vypočítanou aproximaci hodnoty integrálu.

			Funkce {\tt simpson} se samozřejmě během výpočtů funkčních hodnot funkce $f(x)$ odvolává na funkci  {\tt f}. 

			- Funkce {\tt f} je definovaná takto: 
			\begin{verbatim}
	double f(double x)
{
        double function = 3*x*x;    //vlastni matematicka funkce
        return function;
}
			\end{verbatim}

			Zřejmě {\tt 3*x*x} můžeme nahradit libovolnou integrovatelnou funkcí.\\
			- Funkce {\tt analytic} je prakticky stejná jako funkce {\tt f}, s tím rozdílem, že je v ní zapsán, místo matematické funkce, kterou chceme integrovat, analytický tvar integrálu dané matematické funkce. 


		\subsection{Chyba metody}

			Za předpokladu spojitosti čtvrté derivace $f^{(IV)}$ na intervalu $ \langle a,b \rangle $ platí pro velikost chyby metody tento vztah: \cite{nummet1, nummet2}

			\begin{equation}
				\int_{a}^{b}{f(x)dx} - S(f,h) = \frac{1}{180} (b-a) f^{(IV)} (\xi) \cdot h^4
			\label{chyba}
			\end{equation}
			kde $\xi \in (a,b)$.

			Všimněme si, že je-li funkce $P_3(x)$, kterou integrujeme, polynomem třetího stupně, pak je $\frac{d^4}{dx^4} P_3(x) = 0$. Použijeme-li tento výsledek ve vztahu $\ref{chyba}$, vidíme, že chyba aproximace integrálu Simpsonovým pravidlem je nulová pro všechny polynomy do třetího řádu. V ostatních případech je však vztah $\ref{chyba}$ pro určení chyby metody prakticky nepoužitelný, protože nemáme blíže určenou polohu bodů $\xi$, ve kterých počítáme derivace. Při omezenosti příslušných derivací na intervalu $ \langle a,b \rangle $ můžeme ovšem získat horní odhady chyb jednotlivých formulí. Dle skript, ze kterých jsem čerpal \cite{nummet1}, jsou však tyto odhady hodně pesimistické a skutečné chyby aproximací jsou většinou podstatně menší.



	\section{Závěr}
		Program byl dle zadání otestován na takových nekonstantních funkcích a takových intervalech, kdy známý analytický výsledek byl 1, resp. 0. V prvním případě program vypočítal skutečně  hodnotu 1, ve druhém jsme pak obdrželi výsledky například  řádu $10^{-17}$, což odpovídá velikosti strojového epsilon. Pokud byla během výpočtu v jednu chvíli průběžná hodnota Riemanova integrálu například jedna, načež k ní byla přičtena, resp. odečtena, hodnota takto malé velikosti, průběžný výsledek nebyl vzhledem k vlastnostem počítače, resp. charakteru strojového epsilon, změněn, což se v konečném výsledku projevilo takovouto odchylkou. V tomto případě se tedy nejedná přímo o chybu metody a výsledek můžeme považovat za správný, nebo přinejmenším dostačující. Ze vztahu $\ref{chyba}$ je zřejmé, že při zmenšování kroku h se chyba aproximace bude zmenšovat. Použitá literatura \cite{nummet1} mě však upozornila, že přílišné zmenšování nepřinese očekávaný efekt, protože poroste zaokrouhlovací chyba v důsledku narůstání množství výpočtů funkčních hodnot funkce $f$ a dále rostoucího počtu aritmetických operací. Teoreticky bychom se mohli s délkou kroku dostat pod "rozlišovací schopnost" počítače, a tam by mohly všechny integrály vyjít nulové.

\begin{thebibliography}{1}
  \bibitem[1]{nummet1} P. Zörnig:
    \emph{Numerické metody}. České vysoké učení technické v Praze
  \bibitem[2]{nummet2} P. Kocour:
    \emph{Numerické metody}..\\
    \url{http://www.kvd.zcu.cz/cz/materialy/numet/_numet.html}
\end{thebibliography}







\end{document}



