## Plotting and coding the algorithms we learnt in the last class. ### Plotting # $\newcommand{\R}{\mathbb{R}}$ # ### Numpy and Scipy # *Look at [my webpage](http://www.iiserpune.ac.in/~vmallick/2014s1/mth103/ "installing scipy") for instructions to install these python packages on your laptop.* # # # ### Array/Vectors # Suppose you want to plot $y = x^2$, $-1 \leq x \leq 1$. From our math knowledge, we know this is a curve in $\R^2$. The way one plots such functions in matlab/mathematica, is to find a lot of points lying on the curve and then joining them. We can use a list to store those points. So a basic python code to generate such lists will be as follows # In[2]: def genxy(f, xmin, xmax, N) : """Given a function f, an interval [xmin, xmax] in R and number of points N, This function returns a 2-tuple of two lists: First an equispaced points on x-axis, and the second one is the value of f on those points""" step = float(xmax - xmin)/N listx = [xmin + i * step for i in range(N+1)] listy = [f(x) for x in listx] return (listx, listy) # In[5]: a = genxy(lambda x : x*x, -1, 1, 10) azipped = zip(a[0], a[1]) for x, y in azipped : print "%6.3f\t%6.3f" % (x, y) # Out[5]: # -1.000 1.000 # -0.800 0.640 # -0.600 0.360 # -0.400 0.160 # -0.200 0.040 # 0.000 0.000 # 0.200 0.040 # 0.400 0.160 # 0.600 0.360 # 0.800 0.640 # 1.000 1.000 # # This is not too difficult, but is somewhat tedious. To make tasks simpler, one can use arrays. Arrays are **not** in-built in python. To use them, *you have to install numpy*. # In[9]: # The above code becomes : from numpy import * xs = array(a[0]) ys = array(a[1]) type(xs) # Out[9]: # numpy.ndarray # 1. In an array all elements must be of the same type. If they are not, they are converted to the same type as is shown below (in the examples). # 1. The computations are faster if the size of the array doesn't change during computation. # 1. You can mathematical operation on a whole array and that is faster than doing it by traversing a list. # In[27]: # Array elements are all of the same type v1 = [4, 5, 100, 4, 23.4] w1 = array(v1) print w1 # Out[27]: # [ 4. 5. 100. 4. 23.4] # # In[28]: v2 = ['eggs', 4.3, 3] w2 = array(v2) print w2 # Out[28]: # ['eggs' '4.3' '3'] # # In[29]: v3 = [3, 2, 5] w3 = array(v3) print w3 # Out[29]: # [3 2 5] # # # # In[30]: # doing operations on array def f(x) : return x + x*x v = array(range(5)) print f(v) # Out[30]: # [ 0 2 6 12 20] # # In[33]: # One dimensional arrays are like vectors v = array([1, 4, 5]) w = array([5, 2, 1]) print v + w print 2 * v print v * w # Out[33]: # [6 6 6] # [ 2 8 10] # [5 8 5] # # ### Back to plotting. # We need two things : # 1. We need to generate an array of equidistant points on $x$-axis. For this we have a numpy function called **linspace** # 1. We need to generate another array containing their values. # In[35]: xlist = linspace(-1, 1, 11) print xlist # Out[35]: # [-1. -0.8 -0.6 -0.4 -0.2 0. 0.2 0.4 0.6 0.8 1. ] # # In[36]: def square(x) : return x*x ylist = square(xlist) print ylist # Out[36]: # [ 1. 0.64 0.36 0.16 0.04 0. 0.04 0.16 0.36 0.64 1. ] # # Good, it works. Let us now write it as a function. # In[37]: def array_pts(f, xmin, xmax, N) : xlist = linspace(xmin, xmax, N+1) ylist = f(xlist) return (xlist, ylist) # pretty simple, isn't it # ### Plotting # To plot we have a bunch of functions, I like matplotlib. It has a MatLab like syntax. # In[39]: get_ipython().magic(u'matplotlib inline') # In[47]: from matplotlib.pylab import * a = array_pts(lambda x : x*x, -1, 1, 10) plot(a[0], a[1]) show() b = array_pts(sin, -4*pi, 4*pi, 50) plot(b[0], b[1]) show() # Out[47]: # image file: # image file: # ### Some comments on trigonometric functions # Note that we have not imported math anywhere. The sin we used here is imported from numpy. Actually if we override with the sin from math, we'll get errors : # In[48]: from matplotlib.pylab import * import math as m b = array_pts(m.sin, -4*pi, 4*pi, 50) plot(b[0], b[1]) show() # Out[48]: --------------------------------------------------------------------------- TypeError Traceback (most recent call last) in () 2 import math as m 3 ----> 4 b = array_pts(m.sin, -4*pi, 4*pi, 50) 5 plot(b[0], b[1]) 6 show() in array_pts(f, xmin, xmax, N) 1 def array_pts(f, xmin, xmax, N) : 2 xlist = linspace(xmin, xmax, N+1) ----> 3 ylist = f(xlist) 4 return (xlist, ylist) # pretty simple, isn't it TypeError: only length-1 arrays can be converted to Python scalars ### Accessing array elements # Array elements can be accessed in the same way as list elements # x = array(range(1, 10)) # print x[5] # In[49]: x = array(range(1, 10)) print x[5] # Out[49]: # 6 # # However, you must be careful about assignment. Assignment just reflects the same array, and does not create a new array, just like lists: # In[53]: a = array(range(1, 10)) l = range(1, 10) x = a m = l print a[5], l[5] x[5] = -1 m[5] = -1 print a[5], l[5] print type(a), type(l) # Out[53]: # 6 6 # -1 -1 # # # In[55]: # Out[55]: # a before : [ 1 2 3 4 5 -1 7 8 9] # y before : [ 1 2 3 4 5 -1 7 8 9] # a after : [ 1 2 3 4 5 -1 7 8 9] # y after : [ 1 2 -2 4 5 -1 7 8 9] # # ####Note : # For vectors # # a += b # # is faster than # # a = a + b # # since the later creates a new array to store `a + b` before copying it back to `a`. `a += b` just modifies a directly. Same holds for -=, \*= etc. # ### Putting two graphs in he same picture # One uses the hold funciton. # In[59]: def f(x, s): return 1 / (sqrt(2 * pi) * s) * exp(- x * x / (2 * s * s)) x = linspace(-5, 5, 1000) y = f(x, 1) plot(x, y) hold('on') z = f(x, 2) plot(x, z) show() # Out[59]: # image file: # ### Some decorations # In[68]: x = linspace(-5, 5, 1000) y = f(x, 1) plot(x, y) hold('on') z = f(x, 2) plot(x, z) xlabel('x') ylabel('y') legend(['s=1', 's=2']) title('normal or Gaussian distribution') show() # Out[68]: # image file: