%% -*- mode: latex; coding: utf-8 -*- \documentclass[a4paper, 10]{article} % This document is (C) and ©, Mahindra École Centrale, 2014,15 % Written by Lilian Besson (https://bitbucket.org/lbesson/) for the MA 101 course (2014). %% Well understand input, and nice output \usepackage[english]{babel} \usepackage[utf8]{inputenc} \usepackage{lmodern} \usepackage[T1]{fontenc} \usepackage{cmap,enumerate,multicol} % {fixltx2e} \usepackage{amsmath,amsfonts,amssymb,amsthm,bbm,mathabx} \usepackage{mathrsfs,graphicx,xcolor,color,float,framed,scrextend} \everymath{\displaystyle} \usepackage{enumitem} % use {enumerate} if enumitem is not available \usepackage{caption,subcaption} \usepackage{listings} % Include Python code! (and any other language) \usepackage{textcomp} % Thanks to http://tex.stackexchange.com/a/50538 (\textquotesingle{} can be used to input ') \usepackage[scale=0.78]{geometry} \newcommand{\inv}[1]{\frac{1}{#1}} \newcommand{\mytitle}[0]{CS101 project : numerical integration} \usepackage{fancyhdr} \pagestyle{fancyplain} \renewcommand{\headrulewidth}{0.2pt} \renewcommand{\footrulewidth}{0.2pt} %% If necessary, change here the header and footer \lhead{\it{\mytitle{}}} %% Left header \chead{\fancyplain{}{\href{http://www.MahindraEcoleCentrale.edu.in/portal/course/view.php?id=27}{\textbf{CS101}}}} %% Center header \rhead{\it{\today}} %% Right header \lfoot{\small{Please, report any issue to \href{mailto:CS101_at_crans_._org}{\texttt{{CS101}{@}{crans.org}}}}} %% Left footer % \cfoot{\thepage/\pageref{LastPage}} %% Center footer \rfoot{\tt{\href{http://www.MahindraEcoleCentrale.edu.in/}{Mahindra École Centrale}}, 2015} %% Right footer \usepackage[plainpages=false,pdfcenterwindow=true, pdftoolbar=false,pdfmenubar=false,backref, pdftitle={\mytitle{}}, %% If necessary, change here the default PDFAuthor value pdfauthor={Mahindra École Centrale, Hyderabad, Telangana}, linkcolor=black,citecolor=black,filecolor=black,urlcolor=black]{hyperref} %% Allow to add URL hyperlink with the command \href. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Code embedding. %% Custom colors for code embedding. \definecolor{dkgreen}{rgb}{0,0.4,0} \definecolor{gray}{rgb}{0.5,0.5,0.5} \definecolor{purple}{rgb}{0.58,0,0.82} %% From https://en.wikibooks.org/wiki/LaTeX/Source_Code_Listings#Customizing_captions \lstset{ % backgroundcolor=\color{white}, % choose the background color; you must add \usepackage{color} or \usepackage{xcolor} basicstyle=\ttfamily, % \texttt\small, % the size of the fonts that are used for the code, FIXME \ttfamily breakatwhitespace=false, % sets if automatic breaks should only happen at whitespace breaklines=true, % sets automatic line breaking captionpos=b, % sets the caption-position to bottom commentstyle=\small\color{dkgreen}\emph, % comment style % deletekeywords={...}, % if you want to delete keywords from the given language % escapeinside={\%*}{*)}, % if you want to add LaTeX within your code % extendedchar=true, % lets you use non-ASCII characters; for 8-bits encodings only, does not work with UTF-8 frame=single, % adds a frame around the code keywordstyle=\small\color{blue}, % keyword style language=python, % the (default) language of the code fontadjust=false, % if you want to add more keywords to the set % morekeywords={define,domain,objects,init,goal,problem,action,parameters,precondition,effect,types,requirements,strips,typing}, numbers=left, % where to put the line-numbers; possible values are (none, left, right) numbersep=10pt, % how far the line-numbers are from the code numberstyle=\small\color{gray}, % the style that is used for the line-numbers rulecolor=\color{black}, % if not set, the frame-color may be changed on line-breaks within not-black text (e.g. comments (green here)) showspaces=false, % show spaces everywhere adding particular underscores; it overrides 'showstringspaces' showstringspaces=false, % underline spaces within strings only showtabs=false, % show tabs within strings adding particular underscores stepnumber=1, % the step between two line-numbers. If it's 1, each line will be numbered stringstyle=\small\color{purple}, % string literal style tabsize=4, % sets default tabsize to 2 spaces prebreak = \raisebox{0ex}[0ex][0ex]{\ensuremath{\hookleftarrow}}, % pour la fin des lignes. aboveskip={0.5\baselineskip}, belowskip={0.5\baselineskip}, %title=\lstname % show the filename of files included with \lstinputlisting; also try caption instead of title % title=\tiny{File \textcolor{blue}{\url{\lstname}}} % show the filename of files included with \lstinputlisting; also try caption instead of title %% FIXME title ! upquote=true % Thanks to http://tex.stackexchange.com/a/50535 } \begin{document} \numberwithin{equation}{section} %% Option to number equations \DeclareGraphicsExtensions{.pdf,.png,.jpg,.svg} %% Order to load the images \title{\mytitle{}} %% Define Title of the document %% -*- coding:utf8; mode:latex -*- %% LaTeX file automatically generated with [strapdown2pdf](https://bitbucket.org/lbesson/bin/src/master/strapdown2pdf) %% from /home/lilian/cs101/projects/index.html, the jeudi 12 mars 2015, 15:36:29 (UTC+0530). \section*{Instructions for the programming project} This document is explaining to you what you will have to do (and how you will do it) for the CS101 programming project in $2014/15$ at \emph{Mahindra Ecole Centrale}. % \hspace{\fill}\rule{.6\linewidth}{0.4pt}\hspace{\fill} \paragraph{Topic:} You and your team-mates will work on some algorithms used to numerically compute integrals (1D or more). % This is not something new, as we studied some methods in MA101 (October, November), but you will also discover some other approaches. % The goal of this project is to implement correctly some algorithms that will allow us to compute numerical approximations of 1D integrals (area under the curve), and possibly 2D integrals (surface) or 3D integrals (volume). % At the end, you will be able to compute the same integral $\int_{x_{min}}^{x_{max}} f(x) \mathrm{d}x$ with half-a-dozen of different approaches, and compare their efficiency and accuracy. % The user interface we want you to provide will have to be similar to this: \begin{lstlisting}[language=python] import integrals # importing the file integrals.py def f1(x): # here we define an example function return x**2 + x - 4 # Using the function riemann_left that you wrote in integrals.py area = integrals.riemann_left(f1, xmin=0, xmax=10, n=1000) print "For f: x -> x**2 + x - 4, the area under the curve on [0, 10] is approximatively", area \end{lstlisting} \paragraph{General instructions:} Please refer to the main document we gave you, it contains general advices, details about the deadline, and tells you how to submit your work. \hspace{\fill}\rule{.6\linewidth}{0.4pt}\hspace{\fill} \section{What do you have to do?} \subsection{In a file called \texttt{integrals.py}} \begin{enumerate} \item Start by implementing the \textsc{Riemann} sum algorithm, as seen in MA101\footnote{ For this entire project, a good reference is the \textsc{Apostol} Vol.II book, from page 602, paragraph 15.19 to 15.23.}. You will first define a function \texttt{riemann\_left}, that implements the left \textsc{Riemann} rectangle method. Such function can be defined with \texttt{def riemann\_left(f, xmin, xmax, n=1e5)}, with \texttt{f} being the function\footnote{ It should be callable, with one argument returning a numerical value, exactly like the example given above.} for which we want to compute $\int_{xmin}^{xmax} f(x) \mathrm{d}x$, so \texttt{xmin} (or \texttt{x\_min}) and \texttt{xmax} (or \texttt{x\_max}) are respectively left and right bounds of the integral. The last parameter \texttt{n} is the number of rectangles to be used. Remember that the left \textsc{Riemann} rule is stating that $\int_{x_{min}}^{x_{max}} f(x) \mathrm{d}x \simeq \left( \frac{x_{max} - x_{min}}{n} \right) \sum_{k=0}^{n-1} f(x_{k}) $, where $x_k = x_{min} + h \times k $ for $0 \leq k \leq n-1$ (for $n$ big enough). Similarly, implement \texttt{riemann\_right} and \texttt{riemann\_center} for right and centred \textsc{Riemann} sums. \item Similarly, use the \emph{trapezoidal rule} to define a function called \texttt{trapez}. Which of these first 4 functions is the more accurate? Try to compare them on the same function, while increasing the number of rectangles. Give their complexity (as functions of \texttt{n} the number of rectangles). \item The \href{http://code.activestate.com/recipes/577263-numerical-integration-using-monte-carlo-method/}{Monte-Carlo integration} is an completely different algorithm, based on numbers $(x, y)$ randomly picked in a rectangle $[ x_{min}, x_{max} ] \times [ y_{min}, y_{max} ]$ (with $y_{min} \leq \min_{x_{min} \leq x \leq x_{max}} f(x) \leq y_{max}$). Please refer to \href{https://en.wikipedia.org/wiki/Monte_Carlo_integration}{en.wikipedia.org/wiki/Monte\_Carlo\_integration} for more details, examples and illustrations. % The basic concept is to randomly pick a huge number of points in this rectangle, and count how many of them are below the curve $y = f(x)$ (easy by comparing the picked value $y_i$ and the computed value $f(x_i)$). % So if you picked $n$ points, and $k$ are below the curve, the area under the curve is approximatively equal to the area of the rectangle ($(x_{max} - x_{min}) \times (y_{max} - y_{min})$) times the observed probability of getting a point under the curve (the rate $k / n$). Is this method more or less efficient than the \textsc{Riemann} method? If you use $n = 1000$ rectangles and if you pick $n = 1000$ random points, which of the two methods give a more accurate approximation of the integral? \item Implement the \textsc{Simpson} rule, which is similar to the trapezoidal rule (reference is given below). Is it more accurate? Is the $O$ time complexity worse than for the previous methods? \item Finally, use the \emph{bases of plotting with Python}\footnote{ You can go through the first examples of \href{http://www.labri.fr/perso/nrougier/teaching/matplotlib/}{that tutorial} (\href{http://www.labri.fr/perso/nrougier/teaching/matplotlib/}{www.labri.fr/perso/nrougier/teaching/matplotlib/}).} that will be taught soon in lectures, in order to illustrate the different methods. You must illustrate \emph{at least} the 3 \textsc{Riemann} rectangle methods (left, right, centred) on different examples. This part can look like this (you can use that code, and adapt it for right and centred rectangles): \begin{lstlisting}[language=python] from pylab import * # we load the PyLab environment def plot_riemann_left(f, xmin, xmax, name_f, n=10): """ Plot the function f from xmin to xmax, and n left Riemann rectangles.""" figure() # new figure title("Riemman left rectangles for " + name_f + " with " + str(n) + " rectangles.") xvalues = linspace(xmin, xmax, 1000) # values for the x axis yvalues = [ f(x) for x in xvalues ] # and y axis plot(xvalues, yvalues, 'r-') # plots the function as a red line h = (xmax - xmin) / float(n-1) xi = [ xmin + i*h for i in xrange(0, n-1) ] yi = [ f(x) for x in xi ] # left rectangles! # Now plot the rectangles for x, y in zip(xi, yi): # Rectangle (x,0), (x,y), (x+h,y), (x+h,0) plot([x, x, x+h, x+h, x], [0, y, y, 0, 0], 'b-') # One example (to moved to the tests.py file after) print plot_riemann_left(cos, -4, 4, "cos", n=15) \end{lstlisting} Below is included three illustration of the \textsc{Riemann} sum methods (left, centered and right), cf. Fig.\ref{fig:riemann}. \begin{figure} % https://en.wikibooks.org/wiki/LaTeX/Floats,_Figures_and_Captions#Subfloats \centering \begin{subfigure}[b]{0.30\textwidth} \includegraphics[width=\textwidth]{riemannleft} % \caption{Left rectangles} \label{fig:left} \end{subfigure}% ~ %add desired spacing between images, e. g. ~, \quad, \qquad, \hfill etc. %(or a blank line to force the subfigure onto a new line) \begin{subfigure}[b]{0.30\textwidth} \includegraphics[width=\textwidth]{riemanncenter} % \caption{Center rectangles} \label{fig:center} \end{subfigure} ~ %add desired spacing between images, e. g. ~, \quad, \qquad, \hfill etc. %(or a blank line to force the subfigure onto a new line) \begin{subfigure}[b]{0.30\textwidth} \includegraphics[width=\textwidth]{riemannright} % \caption{Right rectangles} \label{fig:right} \end{subfigure} \caption{Three rectangles methods being illustrated with PyLab (MatPlotLib)} \label{fig:riemann} \end{figure} \end{enumerate} \subsection{On-line references, and references in \textsc{Apostol}, Vol.II} \begin{itemize} \item Main references: \href{https://en.wikipedia.org/wiki/Numerical_integration}{en.wikipedia.org/wiki/Numerical\_integration} explains well the problem and give links to different approaches (see also \href{http://mathworld.wolfram.com/NumericalIntegration.html}{mathworld.wolfram.com/NumericalIntegration.html}). \href{https://en.wikipedia.org/wiki/Newton\%E2\%80\%93Cotes_formulas}{en.wikipedia.org/wiki/Newton-Cotes\_formulas} give details about \textsc{Newton-Cotes} formula, the generalization of \textsc{Riemann} and trapezoidal formulas. \item For \textsc{Riemann} rectangles rules, \href{https://en.wikipedia.org/wiki/Riemann_sum\#Methods}{en.wikipedia.org/wiki/Riemann\_sum\#Methods} give details, formulas and illustration of the different methods (left, right, middle). \item % \href{http://mathworld.wolfram.com/TrapezoidalRule.html}{mathworld.wolfram.com/TrapezoidalRule.html} and \href{https://en.wikipedia.org/wiki/Trapezoidal_rule}{en.wikipedia.org/wiki/Trapezoidal\_rule} explains well the trapezoidal rule. \\ You can also refer to the Theorem 15.13, in \textsc{Apostol} Vol.II (page 604). \item Simpson's rule is also well explained on-line\footnote{ Also here \href{https://en.wikipedia.org/wiki/Simpson\%27s_rule}{en.wikipedia.org/wiki/Simpsons\_rule}.} at \href{http://mathworld.wolfram.com/SimpsonsRule.html}{mathworld.wolfram.com/SimpsonsRule.html}. \\ It is also proved and explained on Theorem 15.15, in \textsc{Apostol} Vol.II (page 609). \end{itemize} \subsection{In a file called \texttt{tests.py}} This other file is mandatory, and will show many examples of use of all the algorithms you implemented. \begin{enumerate} \item Write \emph{at least three examples} (of different kinds, like one polynomial, one trigonometric and another function) for \emph{every} algorithm you implemented in the \texttt{integrals.py} file. \item Do not use random or meaningless examples, try to take examples that can be computed mathematically (in order to check the correctness and precision of your results). \item For example, this \texttt{tests.py} file will look like this: \begin{lstlisting}[language=python] import math import integrals area = integrals.riemann(math.sin, 0, math.pi, n=1000) # area should be almost 2.0 print "For sinus, the area under its curve on [0, pi] is approximatively", area # many more to do ! def finv(x): return 1.0 / x area = integrals.montecarlo(finv, 1, 4, n=10000) # area should be almost ln(4) which is approximatively 1.34 print "For the function inverse, the area under its curve on [0, pi] is approximatively", area \end{lstlisting} \item Feel free to get examples from the MA102 course, like the examples seen in lectures, or the problems in either the exams or the tutorial sheets. Any interesting examples that you want to illustrate or compute in your project is most surely welcome. % will help you having a better mark. \end{enumerate} \subsection{In your project report} \begin{itemize} \item Explain any choice or special hypotheses you chose to follow for all the algorithms you implemented. E.g. if one function expect its argument \texttt{f} to be continuous or differentiable (or more), specify it. \item If you choose to implement any extra features to the \texttt{integrals.py} file (it can be an extra function or method), please explain it in details in your report. \item For \emph{each} function and method you write, give its time complexity ($O(n), O(n^2)$ etc) in its \emph{docstring}\footnote{ A good example of such practice is \href{https://bitbucket.org/lbesson/python-demos/src/master/CS101_Second_Mid-Term_Exam__Problem_II.py\#cl-74}{the solution for Problem III for CS second mid term exam} (given on Moodle).}. Here, a difficulty can be to know what is the parameter to count as size of the problem (it can be the precision threshold \texttt{h} or the number of rectangles or random points \texttt{n}). \end{itemize} \hspace{\fill}\rule{.6\linewidth}{0.4pt}\hspace{\fill} \section{Possible extra job?} Any of the following task will give you extra marks\footnote{ It is not possible to aim the best mark without trying any of these$\dots$}. \begin{enumerate} \item Do you have any idea\footnote{ \underline{Hint:} Wikipedia and Internet can help you, but do not copy-paste anything without trying to understand.} of a non-random method that can be used to compute surfaces, volumes or higher-dimension integrals ($2D$, $3D$ or more) ? \item For all your functions and methods, you can take the time to check that the given arguments are valid. For example, \texttt{n} should be non-negative, and \texttt{x\_{max} $\geq$ x\_{min}}.\\ \underline{Hint:} it is easy to perform such checks, by using the \texttt{assert} keyword. Example: \texttt{assert n > 0 and x\_max >= x\_min} is a very easy way to check this. \item Implement any of these additional methods: \begin{itemize} \item These two methods are a little bit harder to program but more efficient : \begin{itemize} \item \textsc{Boole}'s rule (as explained here \href{http://mathworld.wolfram.com/BoolesRule.html}{mathworld.wolfram.com/BoolesRule.html}), \item \textsc{Romberg}'s method (as explained here \href{https://en.wikipedia.org/wiki/Romberg\%27s_method}{en.wikipedia.org/wiki/Romberg's\_method}). \end{itemize} \item The general Gaussian quadrature method is really harder, but not out of your reach, if you have time: \href{http://mathworld.wolfram.com/GaussianQuadrature.html}{mathworld.wolfram.com/GaussianQuadrature.html}. \item Finally, advanced methods such as adaptive quadrature or adaptive \textsc{Simpson} rule can also be considered (\href{https://en.wikipedia.org/wiki/Adaptive_quadrature}{en.wikipedia.org/wiki/Adaptive\_quadrature}). \end{itemize} \end{enumerate} % \subsection{Bonus: exact integration \textbf{(Harder!)}} % Without having to program anything\footnote{ I only expect a small paragraph explaining what you would have tried, in English. No Python is required.}, if your task would have been to compute exact integrals (mathematically correct, not just numerical approximations), what would have been your approach? % What do we need to implement? Any idea can be given on your project report. % \underline{Hint:} this bonus question is quite hard, do not spend too much time on it. \hspace{\fill}\rule{.6\linewidth}{0.4pt}\hspace{\fill} \section{Extra advices} \begin{itemize} \item Each function or method should be correctly implemented. The \emph{correctness} of everything you write \emph{is crucial} (the required examples are in fact here to help you detect any issue). \item Try to be careful about the special cases (empty integrals, division by zero, integer division etc). \item Do not cheat from your friends or from Internet. %%, I will detect it. \item Apply the rule ``the more efficient, the better'' for each function/method you write. Do not try to optimize each line, what is important is the \emph{overall complexity} of the algorithms (e.g. $O(n^2$) is better than $O(n^3)$). You can take some ideas from the exams or the labs. \item But keep in mind that \emph{readability and simplicity of your programs are important}! \item Please contact me if needed (\texttt{{CS101}{@}{crans}{.}{org}} or \texttt{{Lilian}{.}{Besson}{@}{ens}{-}{cachan}{.}{fr}}). \end{itemize} \hspace{\fill}\rule{.6\linewidth}{0.4pt}\hspace{\fill} \newpage{} \section{Symbolic and numerical integration in Python} Of course, all this has already been written in Python. The goal is not to use these packages, but \underline{to write your own \texttt{integrals.py} file, and all these algorithms by yourself.} \emph{If you are curious}, you can compare your results with the ones obtained with one of the following package. These two modules are already included in Anaconda, and should be available directly in Spyder. \subsection{With \textbf{SciPy}/\textbf{NumPy}} The main reference is the \href{http://docs.scipy.org/doc/scipy/}{scipy.integrate} module (Scientific Python, cf. \\ \href{http://docs.scipy.org/doc/scipy/reference/tutorial/integrate.html#integration-scipy-integrate}{docs.scipy.org/doc/scipy/reference/tutorial/integrate.html\#integration-scipy-integrate}). \subsection{With \textbf{SymPy}} Another good reference is the \href{http://docs.sympy.org/latest/}{\texttt{sympy}} (Symbolical Python) module, which provides a general purpose \href{http://docs.sympy.org/latest/modules/integrals/integrals.html#sympy.integrals.integrate}{sympy.integrate} function, designed to compute formal integrals (ie. exact value and not numerical approximations). For more details, the \texttt{sympy.mpmath} package provide implementation of the Gaussian quadrature and other numerical methods (\href{http://docs.sympy.org/latest/modules/mpmath/calculus/integration.html}{docs.sympy.org/latest/modules/mpmath/calculus/integration.html}) . Last link is \href{http://docs.sympy.org/latest/modules/integrals/integrals.html#numeric-integrals}{docs.sympy.org/latest/modules/integrals/integrals.html\#numeric-integrals} which explains how to use and compute \textsc{Gaussian} quadratures. % \paragraph{Warning:} % Do not use these modules to cheat, we will see it: % \textbf{you have to implement everything by yourself.} \end{document}