{ "metadata": { "name": "", "signature": "sha256:e289a3bd3e479dcd0b5cdf328d307db0c70fde2493de2c6cced29d9d19873686" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "[![](https://bytebucket.org/davis68/resources/raw/f7c98d2b95e961fae257707e22a58fa1a2c36bec/logos/baseline_cse_wdmk.png?token=be4cc41d4b2afe594f5b1570a3c5aad96a65f0d6)](http://cse.illinois.edu/)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Scientific Python for Engineers\n", "\n", "\n", "## Scientific Computing\n", "\n", "\n", "In its purest sense, science describes both the simple hypothesize\u2013test\u2013revise method at its heart and the accrued body of knowledge obtained thereby: sharp at the center, fuzzy at the edges. Computation and modeling fits into this cycle by allowing models to be generated and tested against empirical natural results as well as in cases where no physical results are attainable.\n", "\n", "Computational science and engineering is the related discipline, lying at the intersection of applied mathematics, computer science, and traditional science and engineering. As a discipline, CSE uses high-performance computing (HPC), numerical algorithms, and physical insight to further engineering and scientific ends.\n", "\n", "Although most scientific code is written in Fortran or C/C++, Python has rapidly become a serious and viable contender as a complete scientific computing language, particular in competition with MATLAB. The packages discussed in this notebook, NumPy and SciPy, are the current enablers of that prowess. As much of their underlying features take advantage of compiled C libraries, their performance is comparable to traditional compiled languages, yet they retain the dynamic flexibility of Python.\n", "\n", "## Contents\n", "- [Scientific Computing](#intro)\n", "- [NumPy/SciPy](#numpyscipy)\n", "- [Motivation 1](#motiv1)\n", "- [Motivation 2](#motiv2)\n", "- [Numerical Python](#numpy)\n", " - [Arrays & Data](#array)\n", " - [Views, Copies, & Slices](#views)\n", " - [File I/O with `ndarray`s](#fileio)\n", " - [Vectorization](#vector)\n", " - [Linear Algebra](#linalg)\n", " - [Polynomials](#poly)\n", " - [Random Sampling](#random)\n", "- [Scientific Python](#scipy)\n", " - [Special Constants](#const)\n", " - [Integrals & Differential Equations](#intdiff)\n", " - [Data Interpolation](#interp)\n", " - [Optimization](#opt)\n", " - [Special Functions](#sepcfun)\n", " - [Statistics](#stat)\n", "- [Integrated Examples](#int-ex)\n", " - [Finite Difference Heat Equation](#fdhe)\n", " - [Mandelbrot Set (Fractals & Complex Plane)](#mandel)\n", "- [References](#refs)\n", "- [Credits](#credits)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "\n", "## NumPy/SciPy\n", "\n", "This notebook discusses the main features of [Numerical Python](http://www.numpy.org/) (`numpy`), [Scientific Python](http://www.scipy.org/) (`scipy`), and their relatives in an engineering context.\n", "\n", "Although any installation of [IPython](http://ipython.org/) will work with a version of this notebook, we recommend that you download and install either the [Enthought Canopy Distribution](https://www.enthought.com/products/canopy/) (free for academic users) or [Anaconda](https://continuum.io/downloads). To launch the notebook, open a command terminal, type `ipython notebook tutorial.ipynb`, and press Return.\n", "\n", "A few notes on this tutorial:\n", "\n", "- The code in this tutorial is written for Python 3, so if you have Python 2 then you will need to make some sensible modifications.\n", "\n", "- Code blocks starting with `$` are intended to be run on the command line, not executed as Python code." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Standard Header\n", "\n", "As we will be utilizing a number of packages with reasonably long names, we will adopt the _de facto_ standard module abbreviations in the following header. We also ensure that our [division behavior is sensible](http://www.python.org/dev/peps/pep-0238/) by importing from `__future__`: _i.e._, promotion to `double` will occur from `int` or `long` data types involving division: `1/2 == 0.5`. Although this is the default in Python 3, it is a trivial way to help this notebook work in Python 2 if that's what you are using." ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Getting Help" ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "\n", "## Motivation 1: Working with Arrays\n", "\n", "\n", "As a partial motivation for our discussion of NumPy, let us briefly consider the problem of an object moving under the influence of a gravitational force field.\n", "\n", "This problem will obey the standard law for an object subject to gravitational acceleration,\n", "$$y(t) = g t^{2} + v_{0}(t) t + x_0 \\,\\text{.}$$\n", "\n", "We take the object to be initially stationary and falling from an initial position $x_{0} = 1000 \\,\\text{m}$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### For a single point" ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Python has a simple syntax which allows you to express mathematical and numerical expressions succinctly. Its syntax was in part inspired by MATLAB's, although making some distinct changes and improvements.\n", "\n", "- For MATLAB users: what is different in the line `y = a*t**2 + v*t + x0` from what you would typically expect?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### For a list\n", "\n", "Often, however, we don't just want to evaluate a function like this at one point, but at an array of points: for instance, time steps. The intuitive way to do this in Python is to try using a `list`. So take the code snippet above and try to use it on a `list` of numbers:" ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### For an array\n", "\n", "Basic pure Python cannot go beyond this point trivially: you have to introduce loops or some more advanced logic to get this to work properly with lists. A better way is NumPy.\n", "\n", "The NumPy solution to this problem introduces the `array` object (more pedantically, the `ndarray` object). This object has common mathematical operators overloaded to work properly with it out-of-the-box; thus," ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### As a function\n", "\n", "You can generalize this calculation a little bit more by abstracting it out into a function.\n", "\n", "Many features of `array`-based calculations will work as you would intuitively expect in your functions, although sometimes you need to tweak things to make them work correctly. A set of powerful methods in the SciPy library (including NumPy, SciPy, MatPlotLib, and others) allows you to _vectorize_ this function and make it apply to sets of numbers." ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we can trivially utilize NumPy functions, like `arange`, to calculate our result vector." ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Further improvements could be made in the function and the physical model. However, this suffices to make our basic case: NumPy and SciPy greatly extend the power of Python to investigate numerical and physical models. This notebook will examine some of the features and their applications." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "\n", "## Motivation 2: Speed\n", "\n", "Consider a code solving the Laplace equation $\\nabla^2 u = 0$ over a square grid ([source](http://technicaldiscovery.blogspot.com/2011/06/speeding-up-python-numpy-cython-and.html)). The pure Python version follows:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "from numpy import zeros\n", "\n", "dx = 0.1\n", "dy = 0.1\n", "dx2 = dx*dx\n", "dy2 = dy*dy\n", "\n", "def py_update(u):\n", " nx, ny = u.shape\n", " for i in range(1,nx-1):\n", " for j in range(1, ny-1):\n", " u[i,j] = ((u[i+1, j] + u[i-1, j]) * dy2 +\n", " (u[i, j+1] + u[i, j-1]) * dx2) / (2*(dx2+dy2))\n", "\n", "def calc(N, Niter=100, func=py_update, args=()):\n", " u = zeros([N, N])\n", " u[0] = 1\n", " for i in range(Niter):\n", " func(u,*args)\n", " return u" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This code takes a _long_ time (if you don't want to wait, go to Kernel->Interrupt in the menu above). About *** on my machine." ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's compare a NumPy version." ] }, { "cell_type": "code", "collapsed": false, "input": [ "def num_update(u):\n", " u[1:-1,1:-1] = ((u[2:,1:-1]+u[:-2,1:-1])*dy2 + \n", " (u[1:-1,2:] + u[1:-1,:-2])*dx2) / (2*(dx2+dy2))" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So (on my machine) that's (6.18 s)/(21.9 ms) = 282$\\times$ speedup! (You can do even better moving to more optimized code, but we'll ignore that for now.) The point is that for any serious numerical work, NumPy is the package to build on.\n", "\n", "(By the by, you can also do better if you include C code directly, for instance using the `scipy.weave` package which we won't discuss further. In this particular case, it's not much better. We will also discuss the incorporation of C and Fortran code and libraries into Python in a few lessons.)" ] }, { "cell_type": "code", "collapsed": false, "input": [ "from scipy import weave\n", "\n", "def weave_update(u):\n", " code = \"\"\"\n", " int i, j;\n", " for (i=1; i\n", "## Numerical Python\n", "\n", "NumPy is the foundational numerical module in modern Python (both 2 and 3). NumPy provides the basic array manipulation and mathematical routines that are drawn upon by scientific classes. The common classes (which will discuss here) are either built directly on NumPy (SciPy, Pandas, MatPlotLib) or can readily interface with it (`mpmath`)." ] }, { "cell_type": "code", "collapsed": false, "input": [ "from IPython.display import Image\n", "Image(filename='numpy.png', width=800, embed=True)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "NumPy is implemented in combination of Python and C so it executes [much faster](http://wiki.scipy.org/PerformancePython) than the corresponding pure Python code would. In particular, for array-based problems, NumPy provides vastly improved performance, comparable to [MATLAB](http://www.mathworks.com/). NumPy provides the following key features which we will touch on in this notebook:\n", "\n", "- a powerful generic $N$-dimensional array object\n", "\n", "- vectorized functions\n", "\n", "- useful linear algebra, Fourier transform, and random number capabilities (these are greatly extended by SciPy)\n", "\n", "NumPy also contains tools for integrating C/C++ and Fortran code with Python for numerical applications, although this falls beyond the scope of this tutorial." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "## Arrays & Data\n", "\n", "#### Arrays\n", "\n", "The core data type of `numpy` is the [`ndarray`](http://docs.scipy.org/doc/numpy/reference/arrays.ndarray.html), or colloquially the _array_. Arrays are row-major multidimensional containers of items of the same size and type. Or, more intuitively, think of them as lists which can act as mathematical vectors and matrices without any trouble." ] }, { "cell_type": "code", "collapsed": true, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You can make them directly from a list, as above, or you can use one of the built-in helper functions to do it.\n", "\n", "For instance, `numpy.ones` takes the _dimensions_ of the desired array, as either a list or a tuple." ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "These arrays now behave \"properly\"\u2014that is, it's a little harder to run into a mathematical operation that doesn't just work as you may expect." ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Naturally, additions or multiplications for arrays of different shapes do not work." ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As you can see, values promote to the expected data type, although the default is double `float64`. If you need to, you may also specify the desired variable type according to [those available in NumPy](http://docs.scipy.org/doc/numpy/reference/arrays.scalars.html#built-in-scalar-types) (notable examples include double `float64`, complex double `complex128`, quad `float128`, and long `int64`). These reflect the data types available in the underlying C implementation of the numerical code." ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Ranges\n", "\n", "Numerical calculations often require sequential arrays for coordinates. These can be nontrivial to create without a loop in pure Python, but naturally NumPy supports a way\u2014actually two\u2014to create coordinate arrays or ranges trivially.\n", "\n", "##### `arange`\n", "\n", "[`arange`](http://docs.scipy.org/doc/numpy/reference/generated/numpy.arange.html) is the floating-point counterpart to `range`. `range` only returns integers and fails on floating-point arguments. `arange` returns an array of evenly-spaced floating-point values.\n", "\n", "**Python 2 v. 3**: In Python 2, `range` generates a `list` of values. In Python 3, `range` returns an iterable instead (equivalent to Python 2's `xrange`). If you're not sure of the difference yet, the behavior should still be as you generally expect." ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### `linspace`\n", "\n", "Another option you have to create ranges is [`linspace`](http://docs.scipy.org/doc/numpy/reference/generated/numpy.linspace.html), which is familiar to MATLAB users." ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Careful!** `linspace` is actually preferred for most cases, since numerical drift (discussed in a moment) leads `arange` to accrue small numerical truncation errors when using non-integer steps. `np.linspace` is often a better choice, if less conveniently expressed.\n", "\n", "Some other functions you should investigate for creating useful ranges are [`r_`](http://wiki.scipy.org/Numpy_Example_List#r_), [`c_`](http://wiki.scipy.org/Numpy_Example_List#c_), and [`item`](http://wiki.scipy.org/Numpy_Example_List#item)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Exercises\n", "\n", "Use these cells to familiarize yourself with both IPython notebooks and NumPy/SciPy code." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- Create a $5\\times5$ identity matrix." ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- Create a range from -1.7 to 3.4 with 100 intervals using `linspace`. Then do it using `arange`." ] }, { "cell_type": "code", "collapsed": true, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* Create an array `x` with the numbers $0, 0.5, 1.0$.\n", "* Create an array `y` with the numbers $-5.3, -1.8, 1.5$.\n", "* Place `x` and `y` into an array as two subsequent rows of data. (We haven't talked about this yet\u2014try a few things to see what works.)\n", "* Output the result using `print`. _Hint: remember to enclose the variable names in a list._" ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- Remember our earlier example with a falling object? Now we want you to use the documentation to find out how the [`max`](http://docs.scipy.org/doc/numpy/reference/generated/numpy.amax.html#numpy.amax) and [`argmax`](http://docs.scipy.org/doc/numpy/reference/generated/numpy.argmax.html#numpy.argmax) methods work. Modify the last line of the following code to find the maximum height attained by the object and the corresponding time. (Hint: you need to get the index.)" ] }, { "cell_type": "code", "collapsed": false, "input": [ "def y_fall(t, x0, v0):\n", " a = -9.8\n", " y = a*t**2 + v0*t + x0\n", " return y\n", "\n", "t = np.arange(0,300,0.1)\n", "print(y_fall(t, x0, v))" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "### Views, Copies, & Slices\n", "\n", "\n", "When working with an array, we can distinguish between when we want to refer to the same array (or parts of it) by _different names_, or when we want to actually work with copies of data.\n", "\n", "With base data types, simple assignment makes a copy ($y$ doesn't change when $x$ changes)." ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "However, with arrays (particularly large ones), we don't want to incautiously make copies which could rapidly consume our available memory. Thus the behavior for assignment in NumPy is to _not copy an array, but refer to it_." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Simple assignment (no copy)\n", "\n", "\n", "The variable names simply refer to the same object and the same data. Any change made to one happens identically to the other. This is the default behavior for NumPy assignment." ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### View (shallow copy)\n", "\n", "\n", "[_Views_](http://docs.scipy.org/doc/numpy/reference/generated/numpy.ndarray.view.html) are ways of conveniently referring to the same data (perhaps a given row or column) by different variable names. They are constructed using _slicing_ or the `view` method.\n", "\n", "Views particularly find a good usage in referring to slices, such as ghost cells, which require regular access but are inconveniently addressed by a complex slicing expression." ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Copy (deep copy)\n", "\n", "\n", "The `copy` method completely replicates the object and data of the array, creating an independent instance of all data. This operation can be relatively slow and memory-intensive and should sometimes be avoided if possible." ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Slicing\n", "\n", "NumPy slicing simply extends the traditional Python concept to $n$ dimensions. Slices are views of the original array, and so are useful for manipulating or iterating over subsets. The basic format for each dimension is _start_:_stop_:_step_. Any or all of these may be omitted." ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Exercises\n", "- Create a $4 \\times 4$ matrix $A$ consisting of all ones. Take a slice of the first column in a variable called `A_col`. Change the contents of this slice to twos.\n", "- Now make a slice of the second row called `A_row` and change its contents to threes." ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Sorting\n", "\n", "Sorting by a nested element in a list of lists is rather complicated, requiring the definition of a sorting function for `sorted`, for instance. NumPy provides a trivial solution:" ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "\n", "### File I/O with `ndarray`s\n", "\n", "NumPy arrays can be written to disk in either a plain-text or binary format for future use. This is especially convenient for crash recovery (interim matrices can be recovered) and result storage.\n", "\n", "The [`save`](http://docs.scipy.org/doc/numpy/reference/generated/numpy.save.html) function creates a binary `.npy` file containing the array. This file may be loaded again using [`load`](http://docs.scipy.org/doc/numpy/reference/generated/numpy.load.html)." ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The `savetxt` function creates a text file (of user-specified extension) containing the array. This file may be loaded using `loadtxt` or the more robust `genfromtxt` (which handles missing values)." ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Binary formats are preferred over plain text when large data sets are stored." ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "\n", "### Vectorization\n", "\n", "In many cases, existing functions (from outside the NumPy/SciPy family) require modification to work properly with `ndarray`s. This is particularly the case with `if` statements over arrays (which return boolean arrays rather than single values).\n", "\n", "For these cases, NumPy provides the convenience function [`vectorize`](http://docs.scipy.org/doc/numpy/reference/generated/numpy.vectorize.html), which creates an alternate version of the function which can operate on arrays. (The resulting `vectorize`d function is not optimized, however, and should be avoided in performance-critical cases in favor of hand-tuned code.)" ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Why does this fail? Let's evaluate the conditional test in the `if` statement." ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So we get an array out, which the `if` statement doesn't know what to do with. NumPy knows that things like this will occasionally be a problem, so a convenience function, `vectorize`, is provided in order to make things work." ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You can also hand-code a function, for instance using the useful [`where`](http://docs.scipy.org/doc/numpy/reference/generated/numpy.where.html) function (try it):" ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's actually pull `where` out so you can see what it does." ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "heading", "level": 4, "metadata": {}, "source": [ "Exercises" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- What else can go inadvertently wrong? Try the following code. Can it be properly vectorized?" ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "\n", "### Linear Algebra\n", "\n", "[`numpy.linalg`](docs.scipy.org/doc/numpy/reference/routines.linalg.html) will be your workhorse module if you are interested in classical numerics: eigenvalues, matrix solution, decomposition, etc. It has been designed to reflect a MATLAB-like syntax as well for ease of reading and code conversion.\n", "\n", "Many basic features of linear algebra have been promoted from the `numpy.linalg` module to the main `numpy` module and are available in both namespaces." ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Matrix multiplication\n", "\n", "The `*` operator provides elementwise multiplication between two arrays. `numpy.linalg` additionally provides the dot product, inner product, outer product, tensor dot product along specified axes, and support for Einstein summation notation." ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "heading", "level": 4, "metadata": {}, "source": [ "Eigenvalues and eigenvectors" ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "heading", "level": 4, "metadata": {}, "source": [ "Matrix inversion and solution" ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Matrix inversion, while convenient formally, is a tedious memory-hogging process in practice and should be avoided in almost all cases. Use the more intelligent `linalg.solve` command or execute [LU decomposition](http://mathworld.wolfram.com/LUDecomposition.html) (`scipy.linalg.lu`) (where inversion can be done by backsubstituion), etc." ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "heading", "level": 4, "metadata": {}, "source": [ "Other matrix functions" ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "NumPy, by default, uses `ndarray`s, and so `numpy.linalg` is built around these objects. Additionally, `numpy.matlib` provides support for `matrix`, an alternative to `ndarray` which is rarely used but provides an alternative convenience syntax for matrix\u2013matrix multiplication (`*`) and matrix exponentiation (`**`), for instance." ] }, { "cell_type": "heading", "level": 4, "metadata": {}, "source": [ "Truss forces example" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### \u00bb Truss forces\n", "\n", "A common problem in mechanics is the solution of forces in a truss. Tis is solved statically by the _method of joints_, in which an equation is writen for each node of the truss and resulting set of linear equations is solved.\n", "\n", "\n", "\n", "The system of linear equations on the right can be solved in matrix form. Let us write $T x = f$.\n", "\n", "You could write the solution to this as $x = T^{-1}f$. While atractive formally, it is ofen far too expensive to calculate and store an inverse matrix in memory for large problems. Matrix inversion is a brute-force solution to a linear algebra problem. There are a number of clever ways to solve matrices built into NumPy and SciPy. The most frequent method is not to invert the matrix, but instead to use the `solve` function, as above." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- Let $f_1 = 1000$ and $f_2 = 2000$. Write the governing equations in matrix form." ] }, { "cell_type": "code", "collapsed": false, "input": [ "T = np.array([[0.5, 1, 0, 0, 0, 0, 0],\n", " [0.866, 0, 0, 0, 0, 0, 0],\n", " [-0.5, 0, 0.5, 1, 0, 0, 0],\n", " [0.866, 0, 0.866, 0, 0, 0, 0],\n", " [0, -1, -0.5, 0, 0.5, 1, 0],\n", " [0, 0, 0.866, 0, 0.866, 0, 0],\n", " [0, 0, 0, -1, -0.5, 0, 0.5]])\n", "print(T)\n", "f1 = 1000\n", "f2 = 2000\n", "f = np.array([f1, -0.433*f1-0.5*f2, -f1, 0, 0, f2, 0])" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- Solve the matrix using inversion (`numpy.linalg.inv`). How long does it take?" ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- Solve the matrix using the `scipy.linalg.solve` function. How long does it take?" ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The larger the matrix, the more pronounced the difference will be due to the accrual of small numerical errors in different methods." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### \u00bb Sparse matrices\n", "\n", "It is apparent that the matrix in the previous example is largely diagonal: most of the values away from the diagonal are zero. Something similar is often observed in the solution of differential equations using linear algebra. Matrices of this form are referred to as _banded_, and it is more efficient to represent them _sparsely_: that is, we only specify the position and value of nonzero entries and presume that all others are zero. For a $7 \\times 7$ matrix, this is marginally efficient, if it helps at all; for a $10,000 \\times 10,000$ matrix, sparse representation is of obvious utility.\n", "\n", "One typical pattern, the _tridiagonal banded matrix_, is shown below.\n", "\n", "![](https://bytebucket.org/davis68/python/raw/fbd606a038f77b680d585e65726b73951d3e81bf/lessons/img/matrix-sparse-numbered.png?token=b5fc1a22027dc7516eb31b82a3f8af5152ebbab3)\n", "\n", "Other patterns are common for sparse matrices; all have in common that they can be more efficiently represented in sparse form than dense form.\n", "\n", "![](https://bytebucket.org/davis68/python/raw/fbd606a038f77b680d585e65726b73951d3e81bf/lessons/img/matrix-sparse.png?token=ab32fc089f9b20835e110097372f631b3e60756c)\n", "\n", "Let's examine the sparse representation of the truss example." ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To recapitulate, for problems requiring large grids, it often becomes economical to use _sparse matrix representations_. This is the case when most entries of a matrix are zero, meaning that it is more compact to store values and indices rather than the entire matrix. These types of matrices often occur in the solution of finite difference, finite element, and nodal discretizations of differential equations.\n", "\n", "SciPy sparse matrix objects support methods such as `min` and `diagonal`; however, conventional slicing does not generally work with these matrices.\n", "\n", "Finally, recall that sparse matrix representation is designed to be memory-efficient primarily. For numerical efficiency, there are several [representations](http://docs.scipy.org/doc/scipy/reference/sparse.html#sparse-matrix-classes) to choose from which have varying features." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Indices\n", "\n", "Indicial notation is convenient when populating banded matrices, such as for finite difference equations." ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "\n", "### Polynomials\n", "\n", "The next examples will illustrate how to use the built-in [`polynomial`](http://docs.scipy.org/doc/numpy/reference/routines.polynomials.polynomial.html) convenience functions (such as Chebyshev, Laguerre, Legendre, Hermite, etc.).\n", "\n", "First, the representation of polynomials: an intuitive power-series representation is used to carry the coefficients, which may be evaluated using `polyval` or integrated and differentiated, etc. (The index gives the power to which the independent variable is raised.)" ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Normally, you'll probably just use a list of coefficients. Occasionally, however, converting to the `Polynomial` class explicitly can be convenient, as you gain access to a number of methods like `linspace` and `roots`." ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### \u00bb Data fitting\n", "\n", "Naturally, the first (and maybe only) thing most of us want to do with polynomials is use them to fit discrete data. [`polyfit`](http://docs.scipy.org/doc/numpy/reference/generated/numpy.polynomial.polynomial.polyfit.html#numpy.polynomial.polynomial.polyfit) fits the bill.\n", "\n", "Let's pull a data set of a few hundred points and try fitting it. The [Energy Information Administration](http://www.eia.gov/totalenergy/data/monthly/#summary) provides monthly summaries of energy consumption going back to January 1973, so let's pull those values up and see if we can fit them reasonably well." ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Load the energy consumption data from online source.\n", "import urllib.request\n", "\n", "url = urllib.request.urlopen(\"https://raw.githubusercontent.com/uiuc-cse/python-sp15/gh-pages/lessons/data/eia-data.csv\").read(20000).decode('utf-8') # read only 20 000 chars\n", "rows = url.split('\\r') # then split it into lines\n", "data_list = [row.split(',') for row in rows]\n", "data_list = [float(pt) for row in data_list for pt in row ]\n", "energy_data = np.array(data_list)\n", "print(energy_data.shape[0])\n", "energy_data = np.reshape(energy_data, (energy_data.shape[0]/2,2))\n", "\n", "ax = plt.subplot(111)\n", "ax.plot(energy_data[:,0], energy_data[:,1])\n", "ax.set_title(r'U.S. Energy Consumption in $10^{15}$ BTU')\n", "plt.show()" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The data are a bit choppy to fit directly with a polynomial (although you could probably throw in a cyclical multiplier on an annual or semiannual period if you were interested in fitting the data that closely). Let's sample periodically, which is easy in NumPy:" ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Okay, now it's actually trivial to carry out a least-squares fit of the data to whatever polynomial degree they support:" ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The usual caveats about curve fitting, extrapolation, etc., of course apply.\n", "\n", "You can also find out information about the residual and other statistics [as documented](http://docs.scipy.org/doc/numpy/reference/generated/numpy.polynomial.polynomial.polyfit.html#numpy.polynomial.polynomial.polyfit):" ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### \u00bb Basis functions\n", "\n", "The [Chebyshev functions](https://en.wikipedia.org/wiki/Chebyshev_polynomials) are a set of orthogonal basis functions used in the solution of differential equations. They can be recursively defined as\n", "\n", "$$T_0(x) = 1 $$\n", "$$T_1(x) = x $$\n", "$$T_{n+1}(x) = 2xT_n(x) - T_{n-1}(x).$$\n", "\n", "The first few functions are:\n", "\n", "$$T_0(x) = 1$$\n", "\n", "$$T_1(x) = x$$\n", "\n", "$$T_2(x) = 2x^2 - 1$$\n", "\n", "$$T_3(x) = 4x^3 - 3x$$\n", "\n", "$$T_4(x) = 8x^4 - 8x^2 + 1$$\n", "\n", "_et cetera_." ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "[Many other basis functions and standard polynomials](http://docs.scipy.org/doc/numpy/reference/routines.polynomials.classes.html) are available as well." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### \u00bb Taylor series expansions\n", "\n", "The functionality for creating Taylor series expansions of existing functions can be found several places in the SciPy suite. For instance, `scipy.interpolate` provides [`approximate_taylor_polynomial`](http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.interpolate.approximate_taylor_polynomial.html) to estimate the Taylor series around a point in an arbitrary function by polynomial fitting." ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Select an arbitrary function to approximate.\n", "from scipy.special import j0 # Bessel function of the first kind, order zero\n", "\n", "# Approximate the function f by the Taylor series expansion p (translated to the origin, so p(0) = f(x0)).\n", "from scipy.interpolate import approximate_taylor_polynomial\n", "x0 = 3.5 # location at which we will approximate\n", "order = 3 # <-- tuning knob for order of approximation\n", "series = approximate_taylor_polynomial(j0, x0, order, 0.01)\n", "print(series)\n", "\n", "# Plot the approximation(s) versus the original function.\n", "x = np.linspace(x0-0.5, x0+0.5, 200)\n", "\n", "mpl.rcParams['figure.figsize']=[15,6]\n", "fig = plt.figure()\n", "ax = fig.add_subplot(111)\n", "ax.plot(x, j0(x), 'r-', lw=2, label=r'$J_0(x)$')\n", "ax.plot(x, series(x-x0))\n", "\n", "# Plot residual in subfigure.\n", "ax = fig.add_subplot(222)\n", "ax.plot(x, j0(x)-series(x-x0), lw=2)\n", "plt.plot()" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You can also utilize the `sympy.series` function:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "from sympy import Symbol, series\n", "from sympy.functions.special.bessel import besselj\n", "\n", "x = Symbol('x')\n", "series(besselj(0,x), x, x0=x0, n=order+1).evalf()" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "\n", "### Random Sampling\n", "\n", "Random numbers are generated in [`numpy.random`](http://docs.scipy.org/doc/numpy/reference/routines.random.html) using a fairly sophisticated pseudorandom number generation algorithm, the [Mersenne Twister](http://en.wikipedia.org/wiki/Mersenne_twister) (incidentally, it's unsuitable for cryptography but acceptable for Monte Carlo simulations).\n", "\n", "[_Dozens_](http://docs.scipy.org/doc/numpy/reference/generated/numpy.random.RandomState.html#numpy.random.RandomState) of distributions are available for your use. When you require a random value, you sample from a specific distribution with its own limits, weights, and moments: uniform [0,1) (`uniform`), gaussian normal (`normal`), $\\chi^{2}$ (`chisquare`), Pareto (`pareto`), and so on." ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "from numpy.random import normal\n", "\n", "n = (1000,1)\n", "x = normal(size=n)\n", "avg = np.mean(x)\n", "std = np.std(x)\n", "\n", "x_avg = np.ones(n)* avg\n", "x_stdl = np.ones(n)*(avg-std)\n", "x_stdh = np.ones(n)*(avg+std)\n", "\n", "plt.plot(x,'bx',x_avg,'r-',x_stdl,'r--',x_stdh,'r--')\n", "plt.title('%d Random Gaussian Numbers' % n[0]); plt.xlabel(r'$n$'); plt.ylabel(r'$U(n)$')\n", "plt.show()\n", "\n", "plt.hist(x,bins=20)\n", "plt.title('Distribution of %d Random Gaussian Numbers' % n[0]); plt.xlabel(r'$U(n)$'); plt.ylabel(r'Frequency of $U$')\n", "plt.show()" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A _seed_ is necessary to start a PRNG. Typically this is the system clock when the internal `RandomState` object is created, although you may provide a seed manually. In this case, the random number sequence will be reproduced each time the program is invoked. _This is particularly helpful when you are trying to create replicable results or debug your code._" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### \u00bb Experimental Randomized Factorial Design\n", "\n", "A [_factorial experiment_](https://en.wikipedia.org/wiki/Factorial_experiment) tries all combinations or ordered pairs of experimental conditions. It is a best practice to run these trials in a random order, thus attempting to minimize bias and environmental factors.\n", "\n", "You can use the `shuffle` function to trivially generate a randomized factorial design for experimentation. Consider a two-factor experiment where each independent variable can assume four states. There are thus sixteen trials to be held:" ] }, { "cell_type": "code", "collapsed": true, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "\n", "## Scientific Python\n", "\n", "[SciPy](http://scipy.org/) is a collection of mathematical routines built on top of NumPy. SciPy also provides convenience functions for scientific computing. (SciPy can refer to either the entire system of modules around NumPy or specifically to the SciPy library; we consistently take the latter sense in this document.)\n", "\n", "You can get help for a component of `scipy` using the `info` function. Similarly, `source` will let you view the source code for a function." ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "\n", "### Special Constants\n", "\n", "[`scipy.constants`](http://docs.scipy.org/doc/scipy/reference/constants.html) provides mathematical constants and physical constants in SI units (unless otherwise stated). It also provides primitive unit conversion support.\n", "\n", "[`scipy.constants.physical_constants`](http://docs.scipy.org/doc/scipy/reference/constants.html#scipy.constants.physical_constants) is a dictionary of named quantities in the format `physical_constants[name] = (value, unit, uncertainty)`." ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- How many seconds does it take a sound wave to traverse 1,000 miles?" ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Make sure you know what you are asking for." ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "\n", "### Integrals & Differential Equations\n", "#### Numerical integration\n", "\n", "A number of strategies exist for integrating raw numerical data and known functions. We will use the following Bessel function as an illustration." ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Some integration functions require numerical data." ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Others can operate directly on functions and yield error estimates as well. `quad` is the recommended integrator for general-purpose ODEs." ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You can also integrate the convenience classes provided in `numpy.polynomial`. These results are analytic." ] }, { "cell_type": "code", "collapsed": false, "input": [ "#from numpy.polynomial import Polynomial as P\n", "from numpy.polynomial.polynomial import polyint\n", "c = [0.667, 4.675, -2.349, 5.602]\n", "print('The base polynomial coefficients are', c)\n", "print('The integrated polynomial coefficients are', polyint(c))\n", "\n", "x = np.linspace(-2, 2, 1001)\n", "plt.plot(x, np.polyval(c, x))\n", "plt.plot(x, np.polyval(polyint(c), x))\n", "plt.ylim((0, 20))\n", "plt.show()" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Ordinary Differential Equations\n", "\n", "`odeint` provides an interface for the solution of first-order differential equations." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### \u00bb First-order linear ODE\n", "\n", "Solve the equation $\\left( 1 + x^2 \\right) \\textrm{d}y = \\left( \\frac{1}{x} - x y \\right) \\textrm{d}x$.\n", "\n", "Rewrite in a canonical form as follows:\n", "\n", "$\\frac{\\textrm{d}y}{\\textrm{d}x} = \\frac{\\left( \\frac{1}{x} - x y \\right)}{1 + x^2} = \\frac{1}{x + x^3} - \\frac{x}{1 + x^2} y$." ] }, { "cell_type": "code", "collapsed": false, "input": [ "from scipy.integrate import odeint\n", "\n", "def dydx(y, x):\n", " if x == 0:\n", " return float('Inf')\n", " return (1/(x+x**3.0)) - x/(1+x*x)*y\n", "\n", "y0 = np.array((1.0))\n", "dx = 1e-5\n", "x = np.arange(1e-2, 10.0, dx)\n", "sol= odeint(dydx, y0, x)\n", "plt.plot(x, sol)\n", "plt.show()" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### \u00bb First-order nonlinear ODE\n", "\n", "Solve the equation $\\frac{\\textrm{d}y}{\\textrm{d}x} = y(x)^2$ subject to the condition $y(0) = A$ for $A \\in \\left\\{0.2, 0.4, ..., 2.0\\right\\}$ over the range $x \\in \\left[0, \\frac{1}{2}\\right]$." ] }, { "cell_type": "code", "collapsed": false, "input": [ "def dydx(y, x):\n", " return y*y\n", "\n", "y0 = np.arange(0.2,2.01,0.2)\n", "x = np.arange(0.0, 0.5, 1e-5)\n", "sol= odeint(dydx, y0, x)\n", "\n", "import matplotlib.cm as cm\n", "for i in np.arange(10):\n", " c = cm.rainbow(i/10.,1)\n", " plt.plot(x[:], sol[:,i], color=c)\n", "plt.ylim((0,10))\n", "plt.show()" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### \u00bb Second-order linear ODE\n", "\n", "Higher-order ODEs may be solved by the expedient of using a system of coupled DEs thus:\n", "\n", "$$\\frac{du}{dt} = v$$\n", "\n", "$$\\frac{dv}{dt} = f(u, v)$$\n", "\n", "Boundary conditions must be handled using the [shooting method](http://en.wikipedia.org/wiki/Shooting_method), as only one boundary condition can be specified at a time with the current `odeint` function." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Solve the 1D steady-state heat equation $g(x) = \\alpha \\left( \\frac{\\partial^2 u}{\\partial x^2} \\right)$ with $\\alpha = 1$ and the boundary conditions $u(x = 0) = 0$, $u(x = 1) = 1$ over the range $x \\in [0, 1]$ with $g(x) = -sin(x)$.\n", "\n", "Rewrite the equation as a system of two first-order equations:\n", "\n", "
$\\frac{\\partial u}{\\partial x} = v(x)$
\n", "\n", "
$\\frac{\\partial v}{\\partial x} = -20 \\sin(x)$
\n", "\n", "In this case (linear equation), given two trial solutions $u_1$ and $u_2$ meeting only the left-hand boundary condition $u(0)$, we can write the actual solution as a linear combination of these two solutions. The coefficients are:\n", "\n", "
$c_1 = \\frac{u(1) - u_1(1)}{u_1(1) - u_2(1)}$
\n", "\n", "
$c_2 = \\frac{u_2(1) - u(1)}{u_1(1) - u_2(1)}$
" ] }, { "cell_type": "code", "collapsed": false, "input": [ "#f = [u, v]\n", "def dfdx(f, x):\n", " return np.array([f[1], -20*np.sin(x)])\n", "\n", "f1 = odeint(dfdx, np.array([0.0, 0.0]), np.linspace(0.0, 1.0, 10000))\n", "f2 = odeint(dfdx, np.array([0.0, 1.0]), np.linspace(0.0, 1.0, 10000))\n", "#Shooting method for converting a BVP to an IVP.\n", "c1 = (1.0 - f2[:,0][-1]) / (f1[:,0][-1] - f2[:,0][-1])\n", "c2 = (f1[:,0][-1] - 1.0) / (f1[:,0][-1] - f2[:,0][-1])\n", "f = c1 * f1[:,0] + c2 * f2[:,0]\n", "\n", "plt.plot(np.linspace(0.0, 1.0, 10000), f1[:,0], 'r--',\n", " np.linspace(0.0, 1.0, 10000), f2[:,0], 'g--',\n", " np.linspace(0.0, 1.0, 10000), f, 'b-')\n", "plt.title(r'Solution of 1D heat equation $u_{xx}(x) = -20 \\sin(x)$')\n", "plt.ylabel(r'$u(x)$')\n", "plt.xlabel(r'$x$')\n", "plt.ylim((-4,2))\n", "plt.show()" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "\n", "### Data Interpolation\n", "\n", "[`scipy.interpolate`](http://docs.scipy.org/doc/scipy/reference/interpolate.html) provides splines, polynomial interpolators, and other resources for univariate and multivariate interpolation.\n", "\n", "[`interp1d`](http://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp1d.html#scipy.interpolate.interp1d) returns a function interpolating data within the range of the data." ] }, { "cell_type": "code", "collapsed": false, "input": [ "from scipy.interpolate import interp1d\n", "\n", "xpts = np.linspace(0, 10, 11)\n", "y = np.sin(-xpts**2/8.0)\n", "\n", "f = interp1d(xpts, y, kind='linear')\n", "f2 = interp1d(xpts, y, kind='cubic')\n", "\n", "x = np.linspace(0, 10, 201)\n", "plt.plot(xpts, y, 'kx',\n", " x, f(x), 'r--',\n", " x, f2(x), 'r-',\n", " x, np.sin(-x**2/8.0), 'k-')\n", "plt.legend(['data', 'linear', 'cubic', 'exact'], loc='upper left', ncol=2)\n", "plt.ylim((-1.5,1.5))\n", "plt.show()\n", "\n", "print('f(%f) = %f'%(2.71828, f(2.71828)))" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "[`griddata`](http://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.griddata.html#scipy.interpolate.griddata) interpolates unstructured $N$-dimensional data according to one of several methods (nearest, linear, and cubic)." ] }, { "cell_type": "code", "collapsed": false, "input": [ "def func(x, y):\n", " return x*(1-x)*np.cos(4*np.pi*x) * np.sin(2*np.pi*np.sqrt(y))\n", "\n", "# Define the basic grid coordinates.\n", "grid_x, grid_y = np.mgrid[0:1:250j, 0:1:250j]\n", "\n", "# Define a random subset of the grid for which we will generate data.\n", "pts = np.random.rand(500,2)\n", "vals = func(pts[:,0], pts[:,1])\n", "\n", "# Load the methods and generate a grid for each approach.\n", "from scipy.interpolate import griddata\n", "grid_z0 = griddata(pts, vals, (grid_x, grid_y), method='nearest')\n", "grid_z1 = griddata(pts, vals, (grid_x, grid_y), method='linear')\n", "grid_z2 = griddata(pts, vals, (grid_x, grid_y), method='cubic')\n", "\n", "plt.subplot(221)\n", "plt.imshow(func(grid_x, grid_y).T, extent=(0,1,0,1), origin='lower', cmap=cm.PiYG)\n", "plt.plot(pts[:,0], pts[:,1], 'k.', ms=1)\n", "plt.title('Original')\n", "\n", "plt.subplot(222)\n", "plt.imshow(grid_z0.T, extent=(0,1,0,1), origin='lower', cmap=cm.PiYG)\n", "plt.title('Nearest')\n", "\n", "plt.subplot(223)\n", "plt.imshow(grid_z1.T, extent=(0,1,0,1), origin='lower', cmap=cm.PiYG)\n", "plt.title('Linear')\n", "\n", "plt.subplot(224)\n", "plt.imshow(grid_z2.T, extent=(0,1,0,1), origin='lower', cmap=cm.PiYG)\n", "plt.title('Cubic')\n", "\n", "plt.gcf().set_size_inches(12, 12)\n", "plt.show()" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Spline interpolations in one and two dimensions are also available." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "\n", "### Optimization\n", "\n", "As an example, we will optimize over the function $\\sqrt[x]x$, which has a global maximum at $x = e = 2.718281828459045...$. This gives us a straightforward benchmark of testing for convergence and accuracy." ] }, { "cell_type": "code", "collapsed": false, "input": [ "from math import e as ee\n", "def f(x):\n", " return x ** (1/x)\n", "x = np.linspace(0,4,401)\n", "y = f(x)\n", "mpl.rcParams['figure.figsize']=[15,3]\n", "plt.plot(x, y, 'b-')\n", "plt.axis((0.0, 4.0, 0.0, 1.5))\n", "plt.show()" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The only wrinkly is that SciPy only provides a [`minimize`](http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html) function, so we have to subtly modify this to work properly." ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "There is a [wide selection of optimization methods](http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html) available in SciPy, from old standbys like the [Nelder\u2013Mead simplex](http://mathworld.wolfram.com/Nelder-MeadMethod.html) and [conjugate gradient](http://mathworld.wolfram.com/ConjugateGradientMethod.html) to [simulated annealing](http://mathworld.wolfram.com/SimulatedAnnealing.html) and the [dog-leg trust-region](http://www.numerical.rl.ac.uk/people/nimg/course/lectures/raphael/lectures/lec7slides.pdf) algorithm." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "\n", "### Special Functions\n", "\n", "Special functions typically refer to particular mathematical functions with well-established notations and roles, although there is no formal definition of the category.\n", "\n", "Some common special functions include the relatively common trigonometric ($\\sin, \\cos, \\tan, \\sinh, \\cosh, \\tanh$) (which are in the `scipy` namespace), Bessel ($J_\\nu, Y_\\nu, I_\\nu, K_\\nu$), Airy ($\\textrm{Ai}$), spherical harmonics ($Y^m_\\ell$), $\\Gamma$, binomial coefficient, and error ($\\textrm{erf}, \\textrm{erfc}$) functions. Many more esoteric functions, such as Kelvin ($\\textrm{ber}, \\textrm{bei}, \\textrm{ker}, \\textrm{kei},$ etc.), Mathieu ($\\textrm{CE}_n, \\textrm{SE}_n$, etc.), and [hypergeometric functions](https://en.wikipedia.org/wiki/Generalized_hypergeometric_function) ($_pF_q, \\textrm{Li}_2$, etc.) are also supported in [`scipy.special`](http://docs.scipy.org/doc/scipy/reference/special.html)." ] }, { "cell_type": "code", "collapsed": false, "input": [ "from scipy.special import j0, jn_zeros\n", "x = np.linspace(0, 16, 201)\n", "plt.plot(x, j0(x), 'r-', jn_zeros(0, 5), [0]*5, 'bo')\n", "plt.ylim((-1.0, 1.0))\n", "plt.show()" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "from scipy.special import airy\n", "(Ai, Aip, Bi, Bip) = airy(1.0)\n", "print('Ai(1.0) = %.4f\\nBi(1.0) = %.4f'%(Ai, Bi))" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "from scipy.special import erf, gamma, eval_legendre\n", "x = np.linspace(0.001, 3, 201)\n", "plt.plot(x, erf(x), 'r-',\n", " x, gamma(x), 'g-',\n", " x, eval_legendre(8, x), 'b-')\n", "plt.ylim((-0.5, 2.5))\n", "plt.show()" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "There are both generic and optimized specific versions of many of these functions, particularly Bessel and Hankel functions." ] }, { "cell_type": "code", "collapsed": false, "input": [ "from scipy.special import jn, jv, j0\n", "\n", "print('j0(1.5) = %s\\n'%j0(1.5))\n", "\n", "# Let's do some time trials to see which is faster. (This is for 1e6 evaluations, not one.)\n", "import timeit\n", "t = timeit.Timer(\"jv(0, 1.5)\",\"from scipy.special import jv\")\n", "print('jv(0, 1.5) @ %.4s s'%t.timeit())\n", "\n", "t = timeit.Timer(\"jn(0, 1.5)\",\"from scipy.special import jn\")\n", "print('jn(0, 1.5) @ %.4s s'%t.timeit())\n", "\n", "t = timeit.Timer(\"j0(1.5)\",\"from scipy.special import j0\")\n", "print('j0(1.5) @ %.4s s'%t.timeit())" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "\n", "### Basic Statistics\n", "\n", "[`scipy.stats`](http://docs.scipy.org/doc/scipy/reference/stats.html) contains many probability distributions and associated convenience functions and statistical functions." ] }, { "cell_type": "code", "collapsed": false, "input": [ "from scipy.stats import alpha, chi2, maxwell, rayleigh, pareto\n", "x = np.linspace(0.01, 10, 1001)\n", "ax = plt.plot(x, alpha.pdf(x, 1.0), lw=2, label=r'$\\alpha$')\n", "ax = plt.plot(x, chi2.pdf(x, 1.0), lw=2, label=r'$\\chi^2$')\n", "ax = plt.plot(x, maxwell.pdf(x, 1.0), lw=2, label=r'Maxwell')\n", "ax = plt.plot(x, rayleigh.pdf(x, 1.0), lw=2, label=r'Rayleigh')\n", "ax = plt.plot(x, pareto.pdf(x, 1.0), lw=2, label=r'Pareto')\n", "\n", "plt.legend(loc='upper right')\n", "plt.ylim(0.0, 1.2)\n", "plt.show()" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "from numpy.random import normal\n", "n = (1000,1)\n", "x = normal(size=n)\n", "plt.plot(x, 'bo')\n", "\n", "from scipy.stats import tmean, skewtest, histogram\n", "x_mean = np.mean(x)\n", "x_tmean= tmean(x, (-1, 1)) #trimmed mean: mean of interval [-1, 1]\n", "\n", "print('x_mean =\\n', x_mean)\n", "print('x_tmean =\\n', x_tmean)\n", "print('is skew significant? the p-value is', skewtest(x)[1])\n", "\n", "h = histogram(x, 5)\n", "print('histogram bins size = %f, low value = %f, occupancy = %s'%(h[2], h[1], h[0]))" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "\n", "## Integrated Examples" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "##### \u00bb Finite-difference heat equation\n", "\n", "Finite-difference models are used throughout engineering to obtain numerical solutions to differential equations. This particular system models the heat equation\n", "\n", "$$ \\frac{1}{\\alpha} \\frac{\\partial u}{\\partial t} = \\frac{\\partial^2 u}{\\partial x^2}$$\n", "\n", "given an initial condition of $u(x,t=0) = \\sin\\left(\\pi x/L\\right)$ and boundary conditions of $u(x=0,t) = 0$ and $u(x=L,t) = 0$.\n", "\n", "To approximate a derivative, the most straightforward way is to take the formal definition\n", "\n", "$$f'(x) = \\frac{f(x+h)-f(x)}{h}$$\n", "\n", "and use a small but nonzero step $h$ in your calculation.\n", "\n", "Application of this principle to the heat equation leads to a statement of the form\n", "\n", "$$ \\frac{1}{\\alpha} \\frac{u^m_i - u^{m-1}_i}{\\Delta t} = \\frac{u^{m-1}_{i-1} - 2 u^{m-1}_{i} + u^{m-1}_{i+1}}{\\Delta x^2} $$\n", "\n", "or $u^m_i = \\frac{\\alpha \\Delta t}{\\Delta x^2}u^{m-1}_{i-1} + \\left[1 - 2\\left(\\frac{\\alpha \\Delta t}{\\Delta x^2}\\right)\\right]u^{m-1}_{i} + \\frac{\\alpha \\Delta t}{\\Delta x^2}u^{m-1}_{i+1}$.\n", "\n", "This clearly yields a way to calculate subsequent time steps point-by-point from the previous time step's data." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Pure Python Version" ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Basic parameters\n", "nt = 60\n", "nx = 10\n", "alpha = 0.1 \n", "length = 1.0 \n", "tmax = 0.5\n", "\n", "# Derived parameters: mesh spacing and time step size\n", "dx = length / (nx-1)\n", "dt = tmax / (nt-1)\n", "\n", "# Create arrays to save data in process.\n", "x = []\n", "t = []\n", "u = []\n", "for i in range(nt):\n", " t.append(i*dt)\n", " ulist = []\n", " for j in range(nx):\n", " ulist.append(0.0)\n", " u.append(ulist)\n", "for i in range(nx):\n", " x.append(i*dx)\n", "\n", "# Set initial and boundary conditions.\n", "from math import sin, pi\n", "for j in range(nx):\n", " u[0][j] = sin(pi*(j*dx)/length)**2\n", "#boundaries are implicitly set by this initial condition\n", "\n", "# Loop through each time step.\n", "r = alpha * dt / (dx*dx)\n", "s = 1 - 2*r\n", "for n in range(1, nt):\n", " for j in range(1, nx-1):\n", " u[n][j] = r*u[n-1][j-1] + s*u[n-1][j] + r*u[n-1][j+1]\n", "\n", "# Plot the results (initial and final conditions).\n", "mpl.rcParams['figure.figsize']=[15,3]\n", "plt.plot(x, u[0], 'k-', lw=2)\n", "plt.plot(x, u[-1], 'b-', lw=2)\n", "\n", "plt.ylim((0,1))\n", "plt.xlim((0,1))\n", "plt.show()" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### NumPy Version #1\n", "\n", "So the code above basically works, but we can make a few improvements by moving to NumPy for more of it.\n", "\n", "First, let's read in the parameters from a separate file. (There's also a version in `coding.ipynb` which uses the command line arguments as input.) We'll also start using NumPy `ndarray`s for the data, which simplifies things a lot." ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Load parameters from disk.\n", "with open('fd-params.cfg', 'r') as config:\n", " for line in config:\n", " variable, value = [word.strip() for word in line.split('=')]\n", " exec(variable + '=' + value)\n", "\n", "# Compute mesh spacing and time step.\n", "dx = length / (nx-1)\n", "dt = tmax / (nt-1)\n", "\n", "# Create arrays to save data in process.\n", "x = np.linspace(0, length+1e-15, nx)\n", "t = np.linspace(0, tmax+1e-15, nt)\n", "u = np.zeros([nx, nt])\n", "\n", "# Set initial and boundary conditions.\n", "u[:, 0] = np.sin(np.pi*x/length)**2\n", "#boundaries are implicitly set by this initial condition\n", "\n", "# Loop through each time step.\n", "r = alpha * dt / (dx*dx)\n", "s = 1 - 2*r\n", "for n in range(1, nt):\n", " for j in range(1, nx-1):\n", " u[j, n] = r*u[j-1, n-1] + s*u[j, n-1] + r*u[j+1, n-1]\n", "\n", "# Plot the results (initial and final conditions).\n", "mpl.rcParams['figure.figsize']=[15,3]\n", "plt.plot(x, u[:,0], 'k-', lw=2)\n", "plt.plot(x, u[:,-1], 'b-', lw=2)\n", "\n", "plt.ylim((0,1))\n", "plt.xlim((0,1))\n", "plt.show()" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### NumPy Version #2\n", "\n", "Now we can try to convert the calculation to use a matrix instead of a step-by-step `for` loop. We also use a few functions to ease things out." ] }, { "cell_type": "code", "collapsed": false, "input": [ "def calc_params(nx, nt, alpha, length, tmax):\n", " # Compute mesh spacing and time step.\n", " dx = length / (nx-1)\n", " dt = tmax / (nt-1)\n", " return (dx, dt)\n", "\n", "def create_arrays(nx, nt):\n", " # Create arrays to save data in process.\n", " x = np.linspace(0, length+1e-15, nx)\n", " t = np.linspace(0, tmax+1e-15, nt)\n", " u = np.zeros([nx, nt])\n", " return (x, t, u)\n", "\n", "def set_ic(x, length):\n", " # Set initial and boundary conditions.\n", " u[:, 0] = np.sin(np.pi*x/length)**2\n", " #boundaries are implicitly set by this initial condition\n", " return u\n", "\n", "def plot_results(u, x):\n", " # Plot the results (initial and final conditions).\n", " mpl.rcParams['figure.figsize']=[15,3]\n", " plt.plot(x, u[:,0], 'k-', lw=2)\n", " plt.plot(x, u[:,-1], 'b-', lw=2)\n", " \n", " #import matplotlib.cm as cm\n", " #for i in range(len(u[:,0:-1:10])-1):\n", " # c = cm.rainbow(i/len(u[:,0:-1:10]), 1)\n", " # plt.plot(x[:], u[:,i], color=c, lw=2)\n", " \n", " plt.ylim((0,1))\n", " plt.xlim((0,1))\n", " plt.show()\n", "\n", "def gen_matrix(nx, alpha, dt, dx):\n", " r = alpha * dt / (dx*dx)\n", " s = 1 - 2*r\n", " A = np.zeros((nx, nx), dtype=np.float128)\n", " i,j = np.indices(A.shape)\n", " A[i==j] = s\n", " A[i==j-1] = r\n", " A[i==j+1] = r\n", " return A\n", "\n", "# Load parameters from disk.\n", "with open('fd-params.cfg', 'r') as config:\n", " for line in config:\n", " variable, value = [word.strip() for word in line.split('=')]\n", " exec(variable + '=' + value)\n", "\n", "(dx, dt) = calc_params(nx, nt, alpha, length, tmax)\n", "(x, t, u ) = create_arrays(nx, nt)\n", "u = set_ic(x, length)\n", "A = gen_matrix(nx, alpha, dt, dx)\n", "\n", "# Loop through each time step.\n", "for n in range(1, nt):\n", " u[:,n] = A.dot(u[:,n-1])\n", "\n", "plot_results(u, x)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now examine `heat_eqn_cl.py` to see how to read parameters off of the command line.\n", "\n", "Another useful trick is to store parameters in a dictionary for easy lookup, rather than cluttering the namespace with dozens of variables.\n", "\n", "You can see how inputting the parameters from a file can let you handle units automatically as well (_e.g._, `dx = 1e-4 m`)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "##### \u00bb Mandelbrot set (fractals)\n", "\n", "The Mandelbrot set is obtained by iteratively assigning a value to each point $c$ in the complex plane according to the formula $z_{n+1} = z_{n}^2 + c$ ($z_{0} = 0$). If the value $z$ goes to $\\infty$ as $n \\rightarrow \\infty$, then that point $c$ is _not_ part of the Mandelbrot set. Conversely, if the value of $z$ remains bounded no matter how large $n$ becomes, the point $c$ is part of the Mandelbrot set.\n", "\n", "From [Wikipedia](https://en.wikipedia.org/wiki/Mandelbrot_set), \"Mandelbrot set images are made by sampling complex numbers and determining for each whether the result tends towards infinity when a particular mathematical operation is iterated on it. Treating the real and imaginary parts of each number as image coordinates, pixels are colored according to how rapidly the sequence diverges, if at all.\"\n", "\n", "The following code implements the _escape-time algorithm_, which tells you how many iterations until a point \"escapes\", or becomes unbounded, if it does so under a certain maximum limit. This version does not utilize complex numbers, although a version could be written which does so." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Version #1\n", "\n", "This version lumps everything into one function." ] }, { "cell_type": "code", "collapsed": false, "input": [ "def mandelbrot(minR, maxR, minI, maxI, samples=51, iters=25):\n", " \"\"\"\n", " Generate the Mandelbrot set within the boundaries of the limits\n", " maxR, minR, minI, maxI with samples and a maximum iteration of iters.\n", " \"\"\"\n", " # Generate a mesh for the set.\n", " setR = np.linspace(minR, maxR, samples)\n", " setI = np.linspace(minI, maxI, samples)\n", "\n", " # Calculate the values at each point of the mesh by the escape-time\n", " # fractal algorithm.\n", " pts = np.zeros([samples, samples])\n", " for ii in range(1, len(setR)):\n", " for jj in range(1, len(setI)):\n", " it = 0\n", " x = 0.0\n", " y = 0.0\n", "\n", " xx = setR[ii]\n", " yy = setI[jj]\n", "\n", " # Look for escape---i.e., does the value settle down in a few\n", " # iterations or does it keep going?\n", " while(x * x + y * y < 4 and it < iters):\n", " xtmp = x * x - y * y + xx\n", " y = 2 * x * y + yy\n", " x = xtmp\n", " it += 1\n", " pts[ii, jj] = it\n", " return setR, setI, pts\n", "\n", "# Plot boundaries\n", "minR = -2.25\n", "maxR = 0.75\n", "minI = -1.5\n", "maxI = 1.5\n", "\n", "samples = 201\n", "iters = 20\n", "\n", "x, y, z = mandelbrot(minR, maxR, minI, maxI, samples, iters)\n", "z = z.transpose()\n", "\n", "mpl.rcParams['figure.figsize']=[8,8]\n", "plt.imshow(z, interpolation='nearest')\n", "plt.show()" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Version #2\n", "\n", "Let's illustrate the use of a parameter dictionary and split up the functions a little better. We'll also show better form and use the `__main__` check to see if this file is being `import`ed or executed.\n", "\n", "The `mand_alg` function is a great candidate for parallelization as well since each point is independent of the others." ] }, { "cell_type": "code", "collapsed": false, "input": [ "def mand_mesh(params):\n", " \"\"\"Generate a mesh for the Mandelbrot set calculation.\"\"\"\n", " setR = np.linspace(params['r.min'], params['r.max'], params['samples'])\n", " setI = np.linspace(params['i.min'], params['i.max'], params['samples'])\n", " pts = np.zeros([params['samples'], params['samples']])\n", " return (setR, setI, pts)\n", "\n", "def mand_escape(params, xx, yy):\n", " # Look for escape---i.e., does the value settle down in a few\n", " # iterations or does it keep going?\n", " x, y = 0.0, 0.0\n", " it = 0\n", " while(x * x + y * y < 4 and it < params['iter.max']):\n", " z_np1 = x * x - y * y + xx\n", " y = 2 * x * y + yy\n", " x = z_np1\n", " it += 1\n", " return it\n", "\n", "def mand_alg(params):\n", " \"\"\"\n", " Generate the Mandelbrot set within the boundaries of the limits\n", " r.max, r.min, i.min, i.max with samples and a maximum iteration of iter.max.\n", " (These are keys in params dictionary.)\n", " \"\"\"\n", " # Calculate the values at each point of the mesh by the escape-time\n", " # fractal algorithm.\n", " setR, setI, pts = mand_mesh(params)\n", " \n", " for ii in range(1, params['samples']):\n", " for jj in range(1, params['samples']):\n", " xx = setR[ii]\n", " yy = setI[jj]\n", " pts[ii, jj] = mand_escape(params, xx, yy)\n", " return setR, setI, pts\n", "\n", "def mand_plot(z):\n", " mpl.rcParams['figure.figsize']=[8,8]\n", " plt.imshow(z, interpolation='nearest')\n", " plt.show()\n", "\n", "if __name__ == \"__main__\":\n", " params = dict()\n", " params['r.min'] = -2.25\n", " params['r.max'] = 0.75\n", " params['i.min'] = -1.50\n", " params['i.max'] = 1.50\n", " params['samples'] = 401\n", " params['iter.max'] = 20\n", " \n", " x, y, z = mand_alg(params)\n", " z = z.transpose()\n", " mand_plot(z)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "\n", "## References\n", "\n", "- Langtangen, Hans Petter. _Python Scripting for Computational Science_, 3ed. Berlin\u2013Heidelberg: Springer\u2013Verlag, 2009.\n", "- Lugo, Michael. [On propagation of errors](http://gottwurfelt.com/2012/03/26/on-propagation-of-errors/). 26 March 2012.\n", "- Warren, Russell. [A Brief Intro to Profiling in Python](https://speakerdeck.com/rwarren/a-brief-intro-to-profiling-in-python). Ottawa Python Authors Group, 28 February 2013." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "\n", "## Credits\n", "\n", "Neal Davis and Lakshmi Rao developed these materials for [Computational Science and Engineering](http://cse.illinois.edu/) at the University of Illinois at Urbana\u2013Champaign.\n", "\n", "\n", "This content is available under a [Creative Commons Attribution 3.0 Unported License](https://creativecommons.org/licenses/by/3.0/).\n", "\n", "[![](https://bytebucket.org/davis68/resources/raw/f7c98d2b95e961fae257707e22a58fa1a2c36bec/logos/baseline_cse_wdmk.png?token=be4cc41d4b2afe594f5b1570a3c5aad96a65f0d6)](http://cse.illinois.edu/)" ] } ], "metadata": {} } ] }