First introduction to SciPy =========================== The second very important python package we look at is SciPy. The SciPy documentation says: SciPy 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. You can read a longer informal introduction at: http://docs.scipy.org/doc/scipy/reference/tutorial/general.html The home of the SciPy python package is: http://www.scipy.org In the following we will make some examples for a few selected parts of SciPy. But SciPy is a huge library and we can not make an example or point you to every useful function from it. There also is no good reason to duplicate its documentation here. If you encounter a (scientific programming) problem you should go to the library reference: http://docs.scipy.org/doc/scipy/reference/#reference and search for the tools you need. The documentation is very good and well structured. Linear Algebra -------------- At the end of the last chapter we only looked at the very basics of linear algebra. In this section we show some more features. First we have to import the ``linalg`` submodule of ``scipy`` .. sourcecode:: ipython In [1]: from numpy import * In [2]: from scipy.linalg import * Let's do a first computation: .. sourcecode:: ipython In [5]: A = array([[1,2,3],[4,5,6],[7,8,8]]) In [6]: det(A) Out[6]: 3.0 next we can do a LU decomposition of ``A``: .. sourcecode:: ipython In [11]: P, L, U = lu(A) In [12]: L Out[12]: array([[ 1. , 0. , 0. ], [ 0.14285714, 1. , 0. ], [ 0.57142857, 0.5 , 1. ]]) In [13]: U Out[13]: array([[ 7. , 8. , 8. ], [ 0. , 0.85714286, 1.85714286], [ 0. , 0. , 0.5 ]]) In [14]: dot(L,U) Out[14]: array([[ 7., 8., 8.], [ 1., 2., 3.], [ 4., 5., 6.]]) In [15]: dot(P,dot(L,U)) Out[15]: array([[ 1., 2., 3.], [ 4., 5., 6.], [ 7., 8., 8.]]) We can find out the eigenvalues and eigenvectors of this matrix: .. sourcecode:: ipython In [17]: EW, EV = eig(A) In [18]: EW Out[18]: array([ 15.55528261+0.j, -1.41940876+0.j, -0.13587385+0.j]) In [19]: EV Out[19]: array([[-0.24043423, -0.67468642, 0.51853459], [-0.54694322, -0.23391616, -0.78895962], [-0.80190056, 0.70005819, 0.32964312]]) In [20]: dot(A, EV[:,0]) Out[20]: array([ -3.74002233, -8.50785632, -12.47378976]) In [21]: n = dot(A, EV[:,0]) - EW[0]*EV[:,0] In [22]: n Out[22]: array([ -3.55271368e-15+0.j, -1.24344979e-14+0.j, -8.88178420e-15+0.j]) In [23]: norm(n) Out[23]: 1.5688358818848954e-14 Indeed, we found the eigenvalues and eigenvectors. Solving systems of linear equations is not more difficult: .. sourcecode:: ipython In [26]: v = array([[2,3,5]]) In [27]: s = solve(A, transpose(v)) In [28]: s Out[28]: array([[-2.33333333], [ 3.66666667], [-1. ]]) We can check the solution: .. sourcecode:: ipython In [29]: dot(A, s) - v Out[29]: array([[ 4.44089210e-16, -1.00000000e+00, -3.00000000e+00], [ 1.00000000e+00, 8.88178420e-16, -2.00000000e+00], [ 3.00000000e+00, 2.00000000e+00, 3.55271368e-15]]) Oups, broadcasting happened! We have to transpose the vector ``v`` (or squeeze the result of ``dot``): .. sourcecode:: ipython In [30]: dot(A, s) - transpose(v) Out[30]: array([[ 4.44089210e-16], [ 8.88178420e-16], [ 3.55271368e-15]]) If you want you can read on here http://docs.scipy.org/doc/scipy/reference/tutorial/linalg.html to see some more examples on how to do linear algebra with SciPy. For a complete overview of the (dense) linear algebra functionalities of SciPy please see: http://docs.scipy.org/doc/scipy/reference/linalg.html Sparse linear algebra --------------------- SciPy has some routines for computing with sparse and potentially very large matrices. The necessary tools are in the submodule ``scipy.sparse``. We make one example on how to construct a large matrix: .. sourcecode:: ipython In [1]: from scipy.sparse import * In [2]: from scipy import rand In [3]: A = lil_matrix((1000, 1000)) In [4]: A[0,:100] = rand(100) In [5]: A[1,100:200] = A[0,:100] In [6]: A.setdiag(rand(1000)) In [7]: A Out[7]: <1000x1000 sparse matrix of type '' with 1199 stored elements in LInked List format> Next we solve a linear equation: .. sourcecode:: ipython In [8]: from scipy.sparse.linalg import * In [9]: A = A.tocsr() In [10]: A Out[10]: <1000x1000 sparse matrix of type '' with 1199 stored elements in Compressed Sparse Row format> In [11]: b = rand(1000) In [12]: x = spsolve(A, b) The two largest and smallest eigenvalues of ``A``: .. sourcecode:: ipython In [26]: EWs, EVs = eigs(A, 2, which="SM") In [27]: EWs Out[27]: array([ 0.00107174+0.j, 0.00211868+0.j]) In [28]: EWl, EVl = eigs(A, 2, which="LM") In [29]: EWl Out[29]: array([ 0.99854009+0.j, 0.99786107+0.j]) For further information please see the reference for sparse matrices: http://docs.scipy.org/doc/scipy/reference/sparse.html and sparse linear algebra: http://docs.scipy.org/doc/scipy/reference/sparse.linalg.html Fourier Transforms ------------------ There are plenty of different Fourier and related transformations available. Some function are available in the ``numpy.fft`` submodue: http://docs.scipy.org/doc/numpy/reference/routines.fft.html others in the ``scipy.fftpack`` submodule: http://docs.scipy.org/doc/scipy/reference/fftpack.html It's a bit unfortunate that some of these functions live in NumPy rather than SciPy. Optimization and root finding ----------------------------- SciPy has functions for univariate and multivariate as well as unconstrained and constrained optimization techniques. As an example we compute the first 4 maxima of the Bessel function :math:`J_{\frac{1}{2}}(x)` given an initial guess for each one and a bounded interval to search in: .. sourcecode:: ipython In [1]: from scipy.special import jv # Import the Bessel J_n special function In [2]: from scipy.optimize import fminbound In [3]: approx = [1, 6, 12, 18] In [4]: for guess in approx: ...: f = lambda x: -jv(0.5, x) ...: xmax = fminbound(f, guess-4, guess+4) ...: print(xmax) ...: 1.16556123601 7.78988391603 14.1017257582 20.3958424877 This specific Bessel function looks like: .. image:: images/bessel_maxima.png :scale: 100 % Consult the documentation for more details on the various different routines of the ``scipy.optimize`` submodule: http://docs.scipy.org/doc/scipy/reference/optimize.html Other important submodules of SciPy ----------------------------------- - http://docs.scipy.org/doc/scipy/reference/special.html All the ususal higher special functions from mathematics and physics. Almost all of them are implemented as universal functions operating on NumPy arrays. - http://docs.scipy.org/doc/scipy/reference/stats.html Large collection statistical functions and random distributions. - http://docs.scipy.org/doc/scipy/reference/constants.html Physical and mathematical constants and units. The values are taken from the "CODATA Recommended Values of the Fundamental Physical Constants 2010". - http://docs.scipy.org/doc/scipy/reference/integrate.html Some functions for numerical integration and methods for solving ODEs. Further information on SciPy ---------------------------- There are many other submodules of SciPy and we have no place to introduce all. (Additionally not all are relevant for our purpose.) Anyway, if you are interested or need more information, look into the reference guide: - `Scipy Reference Guide `_