\documentclass{article} \usepackage{amsmath} \usepackage{graphicx} \usepackage{color} \usepackage{listings} \usepackage{url} %\definecolor{darkgreen}{rgb}{0,0.5,0} \lstset{language=Python, basicstyle=\ttfamily\bfseries, commentstyle=\color{red}\itshape, stringstyle=\color{green}, showstringspaces=false, keywordstyle=\color{blue}\bfseries} \title{A Glimpse at Scipy} \author{FOSSEE} \date{June 2010} \begin{document} \maketitle \begin{abstract} This document shows a glimpse of the features of Scipy that will be explored during this course. \end{abstract} \section{Introduction} SciPy is open-source software for mathematics, science, and engineering. SciPy (pronounced ``Sigh Pie'') is a collection of mathematical algorithms and convenience functions built on the Numpy extension for Python. It adds significant power to the interactive Python session by exposing the user to high-level commands and classes for the manipulation and visualization of data. With SciPy, an interactive Python session becomes a data-processing and system-prototyping environment rivaling sytems such as \emph{Matlab, IDL, Octave, R-Lab, and Scilab}. \cite{scipy} %% \begin{quote} %% In 1998, ... I came across Python and its numerical extension %% (Numeric) while I was looking for ways to analyze large data sets %% ... using a high-level language. I quickly fell in love with %% Python programming which is a remarkable statement to make about a %% programming language. If I had not seen others with the same view, %% I might have seriously doubted my sanity. -- Travis Oliphant, Creator of Numpy %% \end{quote} \subsection{Sub-packages of Scipy} SciPy is organized into subpackages covering different scientific computing domains. These are summarized in the \underline{table \ref{subpkg}}. \begin{table} \caption{Sub-packages available in Scipy} \label{subpkg} \begin{tabular}{|l|l|} \hline \textbf{Subpackage} & \textbf{Description}\\ \hline \texttt{cluster} & Clustering algorithms\\ \hline \texttt{constants} & Physical and mathematical constants\\ \hline \texttt{fftpack} & Fast Fourier Transform routines\\ \hline \texttt{integrate} & Integration and ordinary differential equation solvers\\ \hline \texttt{interpolate} & Interpolation and smoothing splines\\ \hline \texttt{io} & Input and Output\\ \hline \texttt{linalg} & Linear algebra\\ \hline \texttt{maxentropy} & Maximum entropy methods\\ \hline \texttt{ndimage} & N-dimensional image processing\\ \hline \texttt{odr} & Orthogonal distance regression\\ \hline \texttt{optimize} & Optimization and root-finding routines\\ \hline \texttt{signal} & Signal processing\\ \hline \texttt{sparse} & Sparse matrices and associated routines\\ \hline \texttt{spatial} & Spatial data structures and algorithms\\ \hline \texttt{special} & Special functions\\ \hline \texttt{stats} & Statistical distributions and functions\\ \hline \texttt{weave} & C/C++ integration\\ \hline \end{tabular} \end{table} \subsection{Use of Scipy in this course} Following is a partial list of tasks we shall perform using Scipy, in this course. \begin{enumerate} \item Plotting \footnote{using \texttt{pylab} - see Appendix \ref{mpl}} \item Matrix Operations \begin{itemize} \item Inverse \item Determinant \end{itemize} \item Solving Equations \begin{itemize} \item System of Linear equations \item Polynomials \item Non-linear equations \end{itemize} \item Integration \begin{itemize} \item Quadrature \item ODEs \end{itemize} \end{enumerate} \section{A Glimpse of Scipy functions} This section gives a brief overview of the tasks that are going to be performed using Scipy, in future classes of this course. \subsection{Matrix Operations} Let $\mathbf{A}$ be the matrix \( \begin{bmatrix} 1 &3 &5\\ 2 &5 &1\\ 2 &3 &8 \end{bmatrix} \) To input $\mathbf{A}$ matrix into python, we do the following in \texttt{ipython}\footnote{\texttt{ipython} must be started with \texttt{-pylab} flag}\\ \begin{lstlisting} In []: A = array([[1,3,5],[2,5,1],[2,3,8]]) \end{lstlisting} \subsubsection{Inverse} The inverse of a matrix $\mathbf{A}$ is the matrix $\mathbf{B}$ such that $\mathbf{A}\mathbf{B} = \mathbf{I}$ where $\mathbf{I}$ is the identity matrix consisting of ones down the main diagonal. Usually $\mathbf{B}$ is denoted $\mathbf{B} = \mathbf{A}^{-1}$ . In SciPy, the matrix inverse of matrix $\mathbf{A}$ is obtained using \lstinline+inv(A)+. \begin{lstlisting} In []: inv(A) Out[]: array([[-1.48, 0.36, 0.88], [ 0.56, 0.08, -0.36], [ 0.16, -0.12, 0.04]]) \end{lstlisting} \subsubsection{Determinant} The determinant of a square matrix $\mathbf{A}$ is denoted $\left|\mathbf{A}\right|$. Suppose $a_{ij}$ are the elements of the matrix $\mathbf{A}$ and let $\mathbf{M}_{ij}=\left|\mathbf{A}_{ij}\right|$ be the determinant of the matrix left by removing the $i^{th}$ row and $j^{th}$ column from $\mathbf{A}$. Then for any row $i$ \[ \left|\mathbf{A}\right|=\sum_{j}\left(-1\right)^{i+j}a_{ij}\mathbf{M}_{ij} \] This is a recursive way to define the determinant where the base case is defined by accepting that the determinant of a $1\times1$ matrix is the only matrix element. In SciPy the determinant can be calculated with $det$ . For example, the determinant of \[ \mathbf{A}=\begin{bmatrix} 1 & 3 & 5\\ 2 & 5 & 1\\ 2 & 3 & 8\end{bmatrix}\] is \begin{eqnarray*} |\mathbf{A}| & = & 1\begin{vmatrix} 5 & 1\\ 3 & 8\end{vmatrix} -3\begin{vmatrix} 2 & 1\\ 2 & 8\end{vmatrix} +5\begin{vmatrix}2 & 5\\ 2 & 3\end{vmatrix}\\ & = & 1(5\cdot8-3\cdot1)-3(2\cdot8-2\cdot1)+5(2\cdot3-2\cdot5)=-25 \end{eqnarray*} In SciPy, this is computed as shown below \begin{lstlisting} In []: A = array([[1, 3, 5], [2, 5, 1], [2, 3, 8]]) In []: det(A) Out[]: -25.0 \end{lstlisting} \subsection{Solving Equations} \subsubsection{Linear Equations} Solving linear systems of equations is straightforward using the scipy command \lstinline+solve+. This command expects an input matrix and a right-hand-side vector. The solution vector is then computed. An option for entering a symmetrix matrix is offered which can speed up the processing when applicable. As an example, suppose it is desired to solve the following simultaneous equations: \begin{eqnarray} x+3y+5z & = & 10\\ 2x+5y+z & = & 8\\ 2x+3y+8z & = & 3\end{eqnarray} We could find the solution vector using a matrix inverse: \[ \left[\begin{array}{c} x\\ y\\ z\end{array}\right]=\left[\begin{array}{ccc} 1 & 3 & 5\\ 2 & 5 & 1\\ 2 & 3 & 8\end{array}\right]^{-1}\left[\begin{array}{c} 10\\ 8\\ 3\end{array}\right]=\frac{1}{25}\left[\begin{array}{c} -232\\ 129\\ 19\end{array}\right]=\left[\begin{array}{c} -9.28\\ 5.16\\ 0.76\end{array}\right] \] However, it is better to use the solve command which can be faster and more numerically stable. In this case it however gives the same answer. \begin{lstlisting} In []: A = array([[1, 3, 5], [2, 5, 1], [2, 3, 8]]) In []: b = array([[10], [8], [3]]) In []: dot(inv(A), b) Out[]: array([[-9.28], [ 5.16], [ 0.76]]) In []: solve(A,b) Out[]: array([[-9.28], [ 5.16], [ 0.76]]) \end{lstlisting} \subsubsection{Polynomials} Solving a polynomial is straightforward in scipy using the \lstinline+roots+ command. It expects the coefficients of the polynomial in their decreasing order. For example, let's find the roots of $x^3 - 2x^2 - \frac{1}{2}x + 1$ are $2$, $\sqrt{2}$ and $-\sqrt{2}$. This is easy to see. \begin{align*} x^3 - 2x^2 - \frac{1}{2}x + 1 = 0\\ x^2(x-2) - \frac{1}{2}(x-2) = 0\\ (x-2)(x^2 - \frac{1}{2}) = 0\\ (x-2)(x - \frac{1}{\sqrt{2}})(x + \frac{1}{\sqrt{2}}) = 0 \end{align*} We do it in scipy as shown below: \begin{lstlisting} In []: coeff = array([1, -2, -2, 4]) In []: roots(coeff) \end{lstlisting} \subsubsection{Non-linear Equations} To find a root of a set of non-linear equations, the command \lstinline+fsolve+ is needed. For example, the following example finds the roots of the single-variable transcendental equation \[ x+2\cos\left(x\right)=0,\] and the set of non-linear equations \begin{align} x_{0}\cos\left(x_{1}\right) &= 4,\\ x_{0}x_{1}-x_{1} &= 5 \end{align} The results are $x=-1.0299$ and $x_{0}=6.5041,\, x_{1}=0.9084$ . \begin{lstlisting} In []: def func(x): ...: return x + 2*cos(x) In []: def func2(x): ...: out = [x[0]*cos(x[1]) - 4] ...: out.append(x[1]*x[0] - x[1] - 5) ...: return out In []: from scipy.optimize import fsolve In []: x0 = fsolve(func, 0.3) In []: print x0 -1.02986652932 In []: x02 = fsolve(func2, [1, 1]) In []: print x02 [ 6.50409711 0.90841421] \end{lstlisting} \subsection{Integration} % To be done in the lab. \subsubsection{Quadrature} The function \texttt{quad} is provided to integrate a function of one variable between two points. The points can be $\pm\infty$ ($\pm$ \texttt{inf}) to indicate infinite limits. For example, suppose you wish to integrate the expression $e^{\sin(x)}$ in the interval $[0,2\pi]$, i.e. $\int_0^{2\pi}e^{\sin(x)}dx$, it could be computed using \begin{lstlisting} In []: def func(x): ...: return exp(sin(x)) In []: from scipy.integrate import quad In []: result = quad(func, 0, 2*pi) In []: print result (7.9549265210128457, 4.0521874164521979e-10) \end{lstlisting} \subsubsection{ODE} We wish to solve an (a system of) Ordinary Differential Equation. For this purpose, we shall use \lstinline{odeint}. As an illustration, let us solve the ODE \begin{align} \frac{dy}{dt} = ky(L-y)\\ L = 25000,\,k = 0.00003,\,y(0) = 250 \nonumber \end{align} We solve it in scipy as shown below. \begin{lstlisting} In []: from scipy.integrate import odeint In []: def f(y, t): ...: k, L = 0.00003, 25000 ...: return k*y*(L-y) ...: In []: t = linspace(0, 12, 60) In []: y0 = 250 In []: y = odeint(f, y0, t) \end{lstlisting} Note: To solve a system of ODEs, we need to change the function to return the right hand side of all the equations and the system and the pass the required number of initial conditions to the \lstinline{odeint} function. \appendix \section{Plotting using Pylab}\label{mpl} The following piece of code, produces the plot in Figure \ref{fig:sin} using \texttt{pylab}\cite{pylab} in \texttt{ipython}\footnote{start \texttt{ipython} with \texttt{-pylab} flag}\cite{ipy} \begin{lstlisting} In []: x = linspace(0, 2*pi, 50) In []: plot(x, sin(x)) In []: title('Sine Curve between 0 and $\pi$') In []: legend(['sin(x)']) \end{lstlisting} \begin{figure}[h!] \begin{center} \includegraphics[scale=0.4]{sine} \end{center} \caption{Sine curve} \label{fig:sin} \end{figure} \begin{thebibliography}{9} \bibitem{scipy} Eric Jones and Travis Oliphant and Pearu Peterson and others, \emph{SciPy: Open source scientific tools for Python}, 2001 -- , \url{http://www.scipy.org/} \bibitem{pylab} John D. Hunter, ``Matplotlib: A 2D Graphics Environment,'' \emph{Computing in Science and Engineering}, vol. 9, no. 3, pp. 90-95, May/June 2007, doi:10.1109/MCSE.2007.55 \bibitem{ipy} Fernando Perez, Brian E. Granger, ``IPython: A System for Interactive Scientific Computing,'' \emph{Computing in Science and Engineering}, vol.~9, no.~3, pp.~21-29, May/June 2007, doi:10.1109/MCSE.2007.53. \end{thebibliography} \end{document}