|
1 \documentclass{article} |
|
2 |
|
3 \usepackage{amsmath} |
|
4 \usepackage{graphicx} |
|
5 \usepackage{color} |
|
6 \usepackage{listings} |
|
7 \usepackage{url} |
|
8 %\definecolor{darkgreen}{rgb}{0,0.5,0} |
|
9 |
|
10 \lstset{language=Python, |
|
11 basicstyle=\ttfamily\bfseries, |
|
12 commentstyle=\color{red}\itshape, |
|
13 stringstyle=\color{green}, |
|
14 showstringspaces=false, |
|
15 keywordstyle=\color{blue}\bfseries} |
|
16 |
|
17 \title{A Glimpse at Scipy} |
|
18 \author{FOSSEE} |
|
19 \date{June 2010} |
|
20 |
|
21 \begin{document} |
|
22 \maketitle |
|
23 |
|
24 \begin{abstract} |
|
25 This document shows a glimpse of the features of Scipy that will be |
|
26 explored during this course. |
|
27 \end{abstract} |
|
28 |
|
29 \section{Introduction} |
|
30 |
|
31 SciPy is open-source software for mathematics, science, and |
|
32 engineering. |
|
33 |
|
34 SciPy (pronounced ``Sigh Pie'') is a collection of mathematical |
|
35 algorithms and convenience functions built on the Numpy extension for |
|
36 Python. It adds significant power to the interactive Python session by |
|
37 exposing the user to high-level commands and classes for the |
|
38 manipulation and visualization of data. With SciPy, an interactive |
|
39 Python session becomes a data-processing and system-prototyping |
|
40 environment rivaling sytems such as \emph{Matlab, IDL, Octave, R-Lab, |
|
41 and Scilab}. \cite{scipy} |
|
42 |
|
43 %% \begin{quote} |
|
44 %% In 1998, ... I came across Python and its numerical extension |
|
45 %% (Numeric) while I was looking for ways to analyze large data sets |
|
46 %% ... using a high-level language. I quickly fell in love with |
|
47 %% Python programming which is a remarkable statement to make about a |
|
48 %% programming language. If I had not seen others with the same view, |
|
49 %% I might have seriously doubted my sanity. -- Travis Oliphant, Creator of Numpy |
|
50 %% \end{quote} |
|
51 |
|
52 \subsection{Sub-packages of Scipy} |
|
53 |
|
54 SciPy is organized into subpackages covering different scientific |
|
55 computing domains. These are summarized in the \underline{table |
|
56 \ref{subpkg}}. |
|
57 |
|
58 \begin{table} |
|
59 \caption{Sub-packages available in Scipy} |
|
60 \label{subpkg} |
|
61 \begin{tabular}{|l|l|} |
|
62 \hline |
|
63 |
|
64 \textbf{Subpackage} & \textbf{Description}\\ |
|
65 |
|
66 \hline |
|
67 |
|
68 \texttt{cluster} & Clustering algorithms\\ |
|
69 |
|
70 \hline |
|
71 |
|
72 \texttt{constants} & Physical and mathematical constants\\ |
|
73 |
|
74 \hline |
|
75 |
|
76 \texttt{fftpack} & Fast Fourier Transform routines\\ |
|
77 |
|
78 \hline |
|
79 |
|
80 \texttt{integrate} & Integration and ordinary differential equation |
|
81 solvers\\ |
|
82 |
|
83 \hline |
|
84 |
|
85 \texttt{interpolate} & Interpolation and smoothing splines\\ |
|
86 |
|
87 \hline |
|
88 |
|
89 \texttt{io} & Input and Output\\ |
|
90 |
|
91 \hline |
|
92 |
|
93 \texttt{linalg} & Linear algebra\\ |
|
94 |
|
95 \hline |
|
96 |
|
97 \texttt{maxentropy} & Maximum entropy methods\\ |
|
98 |
|
99 \hline |
|
100 |
|
101 \texttt{ndimage} & N-dimensional image processing\\ |
|
102 |
|
103 \hline |
|
104 |
|
105 \texttt{odr} & Orthogonal distance regression\\ |
|
106 |
|
107 \hline |
|
108 |
|
109 \texttt{optimize} & Optimization and root-finding routines\\ |
|
110 |
|
111 \hline |
|
112 |
|
113 \texttt{signal} & Signal processing\\ |
|
114 |
|
115 \hline |
|
116 |
|
117 \texttt{sparse} & Sparse matrices and associated routines\\ |
|
118 |
|
119 \hline |
|
120 |
|
121 \texttt{spatial} & Spatial data structures and algorithms\\ |
|
122 |
|
123 \hline |
|
124 |
|
125 \texttt{special} & Special functions\\ |
|
126 |
|
127 \hline |
|
128 |
|
129 \texttt{stats} & Statistical distributions and functions\\ |
|
130 |
|
131 \hline |
|
132 |
|
133 \texttt{weave} & C/C++ integration\\ |
|
134 |
|
135 \hline |
|
136 \end{tabular} |
|
137 \end{table} |
|
138 |
|
139 \subsection{Use of Scipy in this course} |
|
140 Following is a partial list of tasks we shall perform using Scipy, in |
|
141 this course. |
|
142 |
|
143 \begin{enumerate} |
|
144 \item Plotting \footnote{using \texttt{pylab} - see Appendix |
|
145 \ref{mpl}} |
|
146 \item Matrix Operations |
|
147 \begin{itemize} |
|
148 \item Inverse |
|
149 \item Determinant |
|
150 \end{itemize} |
|
151 \item Solving Equations |
|
152 \begin{itemize} |
|
153 \item System of Linear equations |
|
154 \item Polynomials |
|
155 \item Non-linear equations |
|
156 \end{itemize} |
|
157 \item Integration |
|
158 \begin{itemize} |
|
159 \item Quadrature |
|
160 \item ODEs |
|
161 \end{itemize} |
|
162 \end{enumerate} |
|
163 \section{A Glimpse of Scipy functions} |
|
164 |
|
165 This section gives a brief overview of the tasks that are going to be |
|
166 performed using Scipy, in future classes of this course. |
|
167 |
|
168 \subsection{Matrix Operations} |
|
169 |
|
170 Let $\mathbf{A}$ be the matrix |
|
171 \( |
|
172 \begin{bmatrix} |
|
173 1 &3 &5\\ |
|
174 2 &5 &1\\ |
|
175 2 &3 &8 |
|
176 \end{bmatrix} |
|
177 \) |
|
178 |
|
179 To input $\mathbf{A}$ matrix into python, we do the following in |
|
180 \texttt{ipython}\footnote{\texttt{ipython} must be started with |
|
181 \texttt{-pylab} flag}\\ |
|
182 |
|
183 \begin{lstlisting} |
|
184 In []: A = array([[1,3,5],[2,5,1],[2,3,8]]) |
|
185 \end{lstlisting} |
|
186 |
|
187 \subsubsection{Inverse} |
|
188 |
|
189 The inverse of a matrix $\mathbf{A}$ is the matrix $\mathbf{B}$ such |
|
190 that $\mathbf{A}\mathbf{B} = \mathbf{I}$ where $\mathbf{I}$ is the |
|
191 identity matrix consisting of ones down the main diagonal. Usually |
|
192 $\mathbf{B}$ is denoted $\mathbf{B} = \mathbf{A}^{-1}$ . In SciPy, the |
|
193 matrix inverse of matrix $\mathbf{A}$ is obtained using |
|
194 |
|
195 \lstinline+inv(A)+. |
|
196 \begin{lstlisting} |
|
197 In []: inv(A) |
|
198 Out[]: |
|
199 array([[-1.48, 0.36, 0.88], |
|
200 [ 0.56, 0.08, -0.36], |
|
201 [ 0.16, -0.12, 0.04]]) |
|
202 \end{lstlisting} |
|
203 |
|
204 \subsubsection{Determinant} |
|
205 |
|
206 The determinant of a square matrix $\mathbf{A}$ is denoted |
|
207 $\left|\mathbf{A}\right|$. Suppose $a_{ij}$ are the elements of the |
|
208 matrix $\mathbf{A}$ and let |
|
209 $\mathbf{M}_{ij}=\left|\mathbf{A}_{ij}\right|$ be the determinant of |
|
210 the matrix left by removing the $i^{th}$ row and $j^{th}$ column from |
|
211 $\mathbf{A}$. Then for any row $i$ |
|
212 |
|
213 \[ \left|\mathbf{A}\right|=\sum_{j}\left(-1\right)^{i+j}a_{ij}\mathbf{M}_{ij} \] |
|
214 |
|
215 This is a recursive way to define the determinant where the base case |
|
216 is defined by accepting that the determinant of a $1\times1$ matrix is |
|
217 the only matrix element. In SciPy the determinant can be calculated |
|
218 with $det$ . For example, the determinant of |
|
219 |
|
220 \[ \mathbf{A}=\begin{bmatrix} 1 & 3 & 5\\ 2 & 5 & 1\\ 2 & 3 & 8\end{bmatrix}\] |
|
221 |
|
222 is |
|
223 |
|
224 \begin{eqnarray*} |
|
225 |\mathbf{A}| & = & 1\begin{vmatrix} 5 & 1\\ 3 & 8\end{vmatrix} |
|
226 -3\begin{vmatrix} 2 & 1\\ 2 & 8\end{vmatrix} |
|
227 +5\begin{vmatrix}2 & 5\\ 2 & 3\end{vmatrix}\\ |
|
228 & = & 1(5\cdot8-3\cdot1)-3(2\cdot8-2\cdot1)+5(2\cdot3-2\cdot5)=-25 |
|
229 \end{eqnarray*} |
|
230 |
|
231 In SciPy, this is computed as shown below |
|
232 |
|
233 \begin{lstlisting} |
|
234 In []: A = array([[1, 3, 5], [2, 5, 1], [2, 3, 8]]) |
|
235 In []: det(A) |
|
236 Out[]: -25.0 |
|
237 \end{lstlisting} |
|
238 |
|
239 \subsection{Solving Equations} |
|
240 |
|
241 \subsubsection{Linear Equations} |
|
242 |
|
243 Solving linear systems of equations is straightforward using the scipy |
|
244 command \lstinline+solve+. This command expects an input matrix and a |
|
245 right-hand-side vector. The solution vector is then computed. An |
|
246 option for entering a symmetrix matrix is offered which can speed up |
|
247 the processing when applicable. As an example, suppose it is desired |
|
248 to solve the following simultaneous equations: |
|
249 |
|
250 \begin{eqnarray} x+3y+5z & = & 10\\ 2x+5y+z & = & 8\\ 2x+3y+8z & = & 3\end{eqnarray} |
|
251 |
|
252 We could find the solution vector using a matrix inverse: |
|
253 |
|
254 \[ \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] \] |
|
255 |
|
256 However, it is better to use the solve command which can be faster and |
|
257 more numerically stable. In this case it however gives the same |
|
258 answer. |
|
259 |
|
260 \begin{lstlisting} |
|
261 In []: A = array([[1, 3, 5], [2, 5, 1], [2, 3, 8]]) |
|
262 In []: b = array([[10], [8], [3]]) |
|
263 In []: dot(inv(A), b) |
|
264 Out[]: |
|
265 array([[-9.28], |
|
266 [ 5.16], |
|
267 [ 0.76]]) |
|
268 |
|
269 In []: solve(A,b) |
|
270 Out[]: |
|
271 array([[-9.28], |
|
272 [ 5.16], |
|
273 [ 0.76]]) |
|
274 \end{lstlisting} |
|
275 |
|
276 \subsubsection{Polynomials} |
|
277 |
|
278 Solving a polynomial is straightforward in scipy using the |
|
279 \lstinline+roots+ command. It expects the coefficients of the |
|
280 polynomial in their decreasing order. For example, let's find the |
|
281 roots of $x^3 - 2x^2 - \frac{1}{2}x + 1$ are $2$, $\sqrt{2}$ and |
|
282 $-\sqrt{2}$. This is easy to see. |
|
283 |
|
284 \begin{align*} |
|
285 x^3 - 2x^2 - \frac{1}{2}x + 1 = 0\\ |
|
286 x^2(x-2) - \frac{1}{2}(x-2) = 0\\ |
|
287 (x-2)(x^2 - \frac{1}{2}) = 0\\ |
|
288 (x-2)(x - \frac{1}{\sqrt{2}})(x + \frac{1}{\sqrt{2}}) = 0 |
|
289 \end{align*} |
|
290 |
|
291 We do it in scipy as shown below: |
|
292 \begin{lstlisting} |
|
293 In []: coeff = array([1, -2, -2, 4]) |
|
294 In []: roots(coeff) |
|
295 \end{lstlisting} |
|
296 |
|
297 \subsubsection{Non-linear Equations} |
|
298 |
|
299 To find a root of a set of non-linear equations, the command |
|
300 \lstinline+fsolve+ is needed. For example, the following example finds |
|
301 the roots of the single-variable transcendental equation |
|
302 |
|
303 \[ x+2\cos\left(x\right)=0,\] |
|
304 |
|
305 and the set of non-linear equations |
|
306 |
|
307 \begin{align} |
|
308 x_{0}\cos\left(x_{1}\right) &= 4,\\ |
|
309 x_{0}x_{1}-x_{1} &= 5 |
|
310 \end{align} |
|
311 |
|
312 The results are $x=-1.0299$ and $x_{0}=6.5041,\, x_{1}=0.9084$ . |
|
313 |
|
314 \begin{lstlisting} |
|
315 In []: def func(x): |
|
316 ...: return x + 2*cos(x) |
|
317 |
|
318 In []: def func2(x): |
|
319 ...: out = [x[0]*cos(x[1]) - 4] |
|
320 ...: out.append(x[1]*x[0] - x[1] - 5) |
|
321 ...: return out |
|
322 |
|
323 In []: from scipy.optimize import fsolve |
|
324 In []: x0 = fsolve(func, 0.3) |
|
325 In []: print x0 |
|
326 -1.02986652932 |
|
327 |
|
328 In []: x02 = fsolve(func2, [1, 1]) |
|
329 In []: print x02 |
|
330 [ 6.50409711 0.90841421] |
|
331 \end{lstlisting} |
|
332 |
|
333 \subsection{Integration} |
|
334 % To be done in the lab. |
|
335 \subsubsection{Quadrature} |
|
336 |
|
337 The function \texttt{quad} is provided to integrate a function of one |
|
338 variable between two points. The points can be $\pm\infty$ ($\pm$ |
|
339 \texttt{inf}) to indicate infinite limits. For example, suppose you |
|
340 wish to integrate the expression $e^{\sin(x)}$ in the interval |
|
341 $[0,2\pi]$, i.e. $\int_0^{2\pi}e^{\sin(x)}dx$, it could be computed |
|
342 using |
|
343 |
|
344 \begin{lstlisting} |
|
345 In []: def func(x): |
|
346 ...: return exp(sin(x)) |
|
347 |
|
348 In []: from scipy.integrate import quad |
|
349 In []: result = quad(func, 0, 2*pi) |
|
350 In []: print result |
|
351 (7.9549265210128457, 4.0521874164521979e-10) |
|
352 \end{lstlisting} |
|
353 |
|
354 \subsubsection{ODE} |
|
355 |
|
356 We wish to solve an (a system of) Ordinary Differential Equation. For |
|
357 this purpose, we shall use \lstinline{odeint}. As an illustration, let |
|
358 us solve the ODE |
|
359 |
|
360 \begin{align} |
|
361 \frac{dy}{dt} = ky(L-y)\\ |
|
362 L = 25000,\,k = 0.00003,\,y(0) = 250 \nonumber |
|
363 \end{align} |
|
364 |
|
365 We solve it in scipy as shown below. |
|
366 |
|
367 \begin{lstlisting} |
|
368 In []: from scipy.integrate import odeint |
|
369 In []: def f(y, t): |
|
370 ...: k, L = 0.00003, 25000 |
|
371 ...: return k*y*(L-y) |
|
372 ...: |
|
373 In []: t = linspace(0, 12, 60) |
|
374 In []: y0 = 250 |
|
375 In []: y = odeint(f, y0, t) |
|
376 \end{lstlisting} |
|
377 |
|
378 Note: To solve a system of ODEs, we need to change the function to |
|
379 return the right hand side of all the equations and the system and the |
|
380 pass the required number of initial conditions to the |
|
381 \lstinline{odeint} function. |
|
382 |
|
383 \appendix |
|
384 |
|
385 \section{Plotting using Pylab}\label{mpl} |
|
386 |
|
387 The following piece of code, produces the plot in Figure \ref{fig:sin} |
|
388 using \texttt{pylab}\cite{pylab} in \texttt{ipython}\footnote{start |
|
389 \texttt{ipython} with \texttt{-pylab} flag}\cite{ipy} |
|
390 |
|
391 \begin{lstlisting} |
|
392 In []: x = linspace(0, 2*pi, 50) |
|
393 In []: plot(x, sin(x)) |
|
394 In []: title('Sine Curve between 0 and $\pi$') |
|
395 In []: legend(['sin(x)']) |
|
396 \end{lstlisting} |
|
397 |
|
398 \begin{figure}[h!] |
|
399 \begin{center} |
|
400 \includegraphics[scale=0.4]{sine} |
|
401 \end{center} |
|
402 \caption{Sine curve} |
|
403 \label{fig:sin} |
|
404 \end{figure} |
|
405 |
|
406 |
|
407 |
|
408 \begin{thebibliography}{9} |
|
409 \bibitem{scipy} |
|
410 Eric Jones and Travis Oliphant and Pearu Peterson and others, |
|
411 \emph{SciPy: Open source scientific tools for Python}, 2001 -- , |
|
412 \url{http://www.scipy.org/} |
|
413 |
|
414 \bibitem{pylab} |
|
415 John D. Hunter, ``Matplotlib: A 2D Graphics Environment,'' |
|
416 \emph{Computing in Science and Engineering}, vol. 9, no. 3, |
|
417 pp. 90-95, May/June 2007, doi:10.1109/MCSE.2007.55 |
|
418 |
|
419 \bibitem{ipy} |
|
420 Fernando Perez, Brian E. Granger, ``IPython: A System for |
|
421 Interactive Scientific Computing,'' \emph{Computing in Science and |
|
422 Engineering}, vol.~9, no.~3, pp.~21-29, May/June 2007, |
|
423 doi:10.1109/MCSE.2007.53. |
|
424 |
|
425 \end{thebibliography} |
|
426 \end{document} |