#!/usr/bin/env python # coding: utf-8 #
# #
# # # Exploratory Computing with Python # *Developed by Mark Bakker* # ## Notebook 4: Functions # In this Notebook we learn how to write our own functions, but we start out with a bit about Python packages. # ### A bit about packages # A package is a set of Python functions. When we want to use a function from a package, we need to import it. There are many different ways to import packages. The most basic syntax is # # `import numpy` # # after which any function in `numpy` can be called as `numpy.function()`. If you don't like the name of the package, for example because it is long, you can change the name. The `numpy` package is renamed to `np` by typing # # `import numpy as np` # # after which all functions in `numpy` can be called as `np.function()`. # # Packages can also have subpackages. For example, the `numpy` package has a subpackage called `random`, which has a bunch of functions to deal with random variables. If the `numpy` package is imported with `import numpy as np`, functions in the `random` subpackage can be called as `np.random.function()`. # # If you only need one specific function, you don't have to import the entire package. For example, if you only want the cosine function of the numpy package, you may import it as `from numpy import cos`, after which you can simply call the cosine function as `cos()`. You can even rename functions when you import them. For example, after `from numpy import cos as newname`, you can call the function `newname()` to compute the cosine (I know, pretty silly, but this can become handy). # # In the previous Notebooks we always imported `numpy` and called it `np` and we imported the `matplotlib.pyplot` and called it `plt`. Both are standard names in the Python community. The statement we added before importing `matplotlib` is `%matplotlib inline`. This latter command is an IPython command and not a Python command. It will only work in IPython and is called a magic command. All magic commands are preceded with a `%`. The statement `%matplotlib inline` puts all figures in the Notebook rather than in a separate window. # # Enough about packages for now. Let's start the way we always start. # In[1]: get_ipython().run_line_magic('matplotlib', 'inline') import numpy as np import matplotlib.pyplot as plt # ### Functions # Functions are an essential part of a programming language. # You already used many functions like `plot`, `loadtxt`, and `linspace`. # But you can also define your own functions. # To define a new function, use the `def` command. After `def` follows the name of the function and then between parentheses the arguments of the function and finally a colon. After the colon you indent until you are done with the function. The last line of the function should be `return` followed by what you want to return. For example, consider the following function of $x$: # # $f(x)= \cos(x) \qquad x <0$ # # $f(x) = \exp(-x) \qquad x \ge 0$ # # Let's implement $f(x)$ in a function called `func`. There is one input argument: $x$. # In[2]: def func(x): if x < 0: f = np.cos(x) else: f = np.exp(-x) return f print(func(3)) # Once you define a function in Python, you can call it whenever you want during the session. So we can call it again # In[3]: print(func(-2)) # If you type # # `func(` and then hit [shift-tab] # # and wait a split second, the input arguments of the function pop-up in a little window, just like for other functions we already used. You can also provide additional documentation of your function. Put the documentation at the top of the indented block and put it between triple double quotes (`"""`). Run the code below to define the function `func` with the additional documentation, then in the code cell below type # # `func(` # # and hit [shift][tab] to see the additional documentation. Warning: don't leave a code cell with just `func(` or `func()` as you will get an error on [Kernel][Restart & Run All Cells]. # In[4]: def func(x): """First Python function written by Student X""" if x < 0: f = np.cos(x) else: f = np.exp(-x) return f # In[ ]: # The names of the arguments of a function are the names used inside the function. They have no relationship to the names used outside the function. When using a variable as the argument of a function, only the *value* gets passed to the function. In the example below, the *value* of `y` is passed as the first argument of the function `func`. Inside the function, this value is used for the variable `x`. # In[5]: y = 2 print('func(2):', func(y)) # ### Exercise 1. First function # Write a Python function for the following function: # # $f(x)=e^{-\alpha x}\cos(x)$ # # The function should take `x` and `alpha` as input arguments and return the function value. Give your function a unique name (if you also call it `func` it will overwrite the `func` function that we defined above). Make a plot of $f(x)$ vs. $x$ for $x$ going from 0 to $10\pi$ using two different values of $\alpha$: 0.1 and 0.2. Add a legend and label the axes. # In[ ]: # Answer to Exercise 1 # ### Keyword arguments # Functions may have multiple input arguments followed by keyword arguments. Arguments *must* be entered and must be entered in the order defined. Keyword arguments don't need to be entered. When they are not entered, the default value is used. Keyword arguments may be given in any order as long as they come after the regular arguments. If you specify the keyword arguments in the order they are defined in the argument list, you don't even need to preceed them with the keyword, but it is saver to write the keywords out and it makes your code easier to read. For example, the function $f(x)=A\cos(\pi x+\theta)$ can be written with keyword arguments for $A$ and $\theta$ as follows. # In[6]: def testfunc(x, A=1, theta=0): return A * np.cos(np.pi * x + theta) print(testfunc(1)) # Uses default A=1, theta=0: cos(pi) print(testfunc(1, A=2)) # Now A=2, and theta is still 0: 2*cos(pi) print(testfunc(1, A=2, theta=np.pi / 4)) # Now A=2, theta=pi/4: 2*cos(5pi/4) print(testfunc(1, theta=np.pi / 4, A=2)) # Same as above: 2*cos(5pi/4) print(testfunc(1, theta=np.pi / 4)) # Now theta=pi/4, and A is still 1: cos(5pi/4) # Note that the proper style was applied, as defined in Notebook 1: there are spaces around mathematical symbols, but not around the equal sign of the keyword argument. # ### Local variables # Variables declared inside a function can only be used inside that function. The outside of a function doesn't know about the variables used inside the function, except for the variables that are returned by the function. In the code below, remove the `#` before `print(a)` and you will get an error message, as `a` is a local variable inside the function `localtest` (then put the `#` back, else you get an error when running [Kernel][Restart & Run All Cells]). # In[7]: def localtest(x): a = 3 b = 5 return a * x + b print(localtest(4)) #print(a) # Will cause an error, as 'a' is not known outside function # ### Three types of variables inside a function # There are actually three types of variables inside a function. We already learned about two of them: variables passed to the function through the argument list, like `x` in the function above, and local variables, like `a` and `b` in the function above. The third type are variables defined outside the function but not passed to the function through the argument list. When a variable is used inside a Python function, Python first checks whether the variable has been defined locally. If not, it checks whether the variable is passed to the function through the argument list. And if that is not the case, Python checks whether the variable is defined outside the function, from the place the function was called. If that is not the case either, it will throw an error message. It is considered good coding practice to pass variables to a function when they are needed inside a function, rather than counting on Python to *find* the variable outside the function; it will likely lead to fewer coding errors as well. # # Note that when a variable is defined locally, Python will not check whether that variable is also declared outside the function. It will happily create a new variable with the same name inside the function. It is important to realize the difference between these different types, so let's do a few examples. # In[8]: # This function works properly def test1(x): a = 3 b = 5 return a * x + b print(test1(4)) # This function also works, but it is sloppy coding # since variable a is defined outside the function a = 3 def test2(x): b = 5 return a * x + b print(test2(4)) # In the following function, we define variable `var1` outside the function `test3`. The function `test3` doesn't take any input arguments (but it still needs the parentheses, else Python doesn't know it is a function!), and it creates a local variable `var1`. This local `var1` variable is only known inside the function `test3` and doesn't effect the value of `var1` outside function `test3`. # In[9]: var1 = 8 def test3(): var1 = 4 print('Inside the function test3, var1 equals:', var1) test3() print('Value of var1 outside test3:', var1) # ### Functions are building blocks that need to be tested separately # Functions are the building blocks of a computer code. They represent a well-defined functionality, which means they can *and should* be tested separately. So make it a habit to test whether your function does what you intended it to do. Sometimes it is easy to test a function: you can compare the value to a hand-calculation, for example. Other times it is more difficult, and you need to write some additional code to test the function. It is always worthwhile to do that. If you test your functions well, it will aid you in debugging your code, because you know that the error is not inside the function. # ### Exercise 2, Stream function for flow around a cylinder # Consider two-dimensional inviscid fluid flow (potential flow) around a cylinder. # The origin of the coordinate system is at the center of the cylinder. # The stream function is a function that is constant along stream lines. # The stream function $\psi$ is a function of polar coordinates $r$ and $\theta$. The stream function is constant and equal to zero on the cylinder and doesn't really exist inside the cylinder, so let's make it zero there, like it is on the cylinder. # # $$\begin{split} # \psi &= 0 \qquad r\le R \\ # \psi &= U(r-R^2/r)\sin(\theta) \qquad r\ge R # \end{split}$$ # # where $U$ is the flow in the $x$-direction, $r$ is the radial distance from the center of the cylinder, $\theta$ is the angle, and $R$ is the radius of the cylinder. You may recall it is not always easy to compute the correct angle when given a value of $x$ and $y$, as the regular arctan function returns a value between $-\pi/2$ and $+\pi/2$ (radians), while if $x=-2$ and $y=2$, the angle should be $3\pi/4$. # `numpy` has a very cool function to compute the correct angle between $-\pi$ and $+\pi$ given the $x$ and $y$ coordinates. The function is `arctan2(y, x)`. Note that the function takes as its *first* argument `y` and as its *second* argument `x`. # # Write a function that computes the stream function for flow around a cylinder. The function should take two arguments, `x` and `y`, and two keyword arguments, `U` and `R`, and should return the stream function value. If you write the function correctly, it should give `psi(2, 4, U=2, R=1.5) = 7.1`, and `psi(0.5, 0, U=2, R=1.5) = 0` (inside the cylinder). # In[ ]: # Answer to Exercise 2 # ### Vectorization of a function # Not all functions can be called with an array of values as input argument. For example, the function `func` defined at the beginning of this notebook doesn't work with an array of `x` values. Remove the `#` and try it out. Then put the `#` back # In[10]: def func(x): if x < 0: f = np.cos(x) else: f = np.exp(-x) return f x = np.linspace(-6, 6, 100) #y = func(x) # Run this line after removing the # to see the error that occurs. Then put the # back # The reason this doesn't work is that Python doesn't know what to do with the line # # `if x < 0` # # when `x` contains many values. Hence the error message # # `The truth value of an array with more than one element is ambiguous` # # For some values of `x` the `if` statement may be `True`, for others it may be `False`. A simple way around this problem is to vectorize the function. That means we create a new function, let's call it `funcvec`, that is a vectorized form of `func` and can be called with an array as an argument. This is by far the easiest but not necessarily the computationally fastest way to make sure a function can be called with an array as an argument, and, unfortunately, it won't work for all situations. # In[11]: funcvec = np.vectorize(func) x = np.linspace(-6, 6, 100) y = funcvec(x) plt.plot(x, y); # Back now to the problem of flow around a clinder. Contours of the stream function represent stream lines around the cylinder. To make a contour plot, the function to be contoured needs to be evaluated on a grid of points. The grid of points and an array with the values of the stream function at these points can be passed to a contouring routine to create a contour plot. To create a grid of points, use the function `meshgrid` which takes as input an array of `x` values and an array of `y` values, and returns a grid of `x` values and a grid of `y` values. For example, to have 5 points in the $x$-direction from -1 to +1, and 3 points in y-direction from 0 to 10: # In[12]: x,y = np.meshgrid(np.linspace(-1, 1, 5), np.linspace(0, 10, 3)) print('x values') print(x) print('y values') print(y) # ### Exercise 3, Contour plot for flow around a cylinder # Evaluate the function for the stream function around a cylinder with radius 1.5 on a grid of 100 by 100 points, where `x` varies from -4 to +4, and `y` varies from -3 to 3; use $U=1$. Evaluate the stream function on the entire grid (you need to create a vectorized version of the function you wrote to compute the stream function). Then use the `np.contour` function to create a contour plot (find out how by reading the help of the `contour` function or go to [this demo](http://matplotlib.org/examples/pylab_examples/contour_demo.html)) of the `matplotlib` gallery). You need to use the command `plt.axis('equal')`, so that the scales along the axes are equal and the circle looks like a circle rather than an ellipse. Finally, you may want to add a nice circular patch using the `fill` command and specifying a bunch of $x$ and $y$ values around the circumference of the cylinder. # In[ ]: # Answer to Exercise 3 # ### Return multiple *things* # An assignment can assign values to multiple variables in one statement, for example # In[13]: a, b = 4, 3 print('a:', a) print('b:', b) a, b, c = 27, np.arange(4), 'hello' print('a:', a) print('b:', b) print('c:', c) d, e, f = np.arange(0, 11, 5) print('d:', d) print('e:', e) print('f:', f) # Similarly, a function may return one value or one array. Or a function may return multiple values, multiple arrays, or whatever the programmer decides to return (including nothing, of course). When multiple *things* are returned, they are returned as a tuple. They can be stored as a tuple, or, if the user knows how many *things* are returned, they can be stored in individual variables right away, as in the example below. # In[14]: def newfunc(): dump = 4 * np.ones(5) dump[0] = 100 return 33, dump, 'this works great!' test = newfunc() print(type(test)) print(test[1]) a, b, c = newfunc() print('a:', a) print('b:', b) print('c:', c) # ### Exercise 4, Streamplot of flow around a cylinder # The radial and tangential components of the velocity vector $\vec{v}=(v_r,v_\theta)$ for inviscid fluid flow around a cylinder are given by # # $\begin{split} # v_r&=U(1-R^2/r^2)\cos(\theta) \qquad r\ge R \\ # v_\theta&=-U(1+R^2/r^2)\sin(\theta) \qquad r\ge R # \end{split}$ # # and is zero otherwise. The $x$ and $y$ components of the velocity vector may be obtained from the radial and tangential components as # # $\begin{split} # v_x&=v_r\cos(\theta) - v_\theta\sin(\theta) \\ # v_y &= v_r\sin(\theta) + v_\theta\cos(\theta) # \end{split}$ # # Write a function that returns the $x$ and $y$ components of the velocity vector for fluid flow around a cylinder with $R=1.5$ and $U=2$. # Test your function by making sure that at $(x,y) = (2,3)$ the velocity vector is $(v_x,v_y)=(2.1331, -0.3195)$. # Compute the $x$ and $y$ components of the velocity vector (vectorization won't help here, as your function returns two values, so you need a double loop) on a grid of 50 by 50 points where `x` varies from -4 to +4, and `y` varies from -3 to 3. Create a stream plot using the cool function `plt.streamplot`, which takes four arguments: `x`, `y`, `vx`, `vy`. # In[ ]: # Answer to Exercise 4 # ### Exercise 5, Derivative of a function # The function `func`, which we wrote earlier in this notebook, implements the following function # # $f(x)= \cos(x) \qquad x <0$ # # $f(x) = \exp(-x) \qquad x \ge 0$ # # Derive an analytic expression (by hand) for the first derivative of $f(x)$ and implement it in a Python function. Test your function by comparing its output to a numerical derivative using a central difference scheme # # $\frac{\text{d}f}{\text{d}x}\approx \frac{f(x+d)-f(x-d)}{2d}$ # # where $d$ is a small number. Test your function for both $x<0$ and $x>0$. # In[ ]: # Answer to Exercise 5 # ### Using a function as the argument of another function # So far, we have used single values or arrays as input arguments of functions. But we can also use a function as one of the input arguments of another function. Consider, for example, a function called `takesquare` that takes two input arguments: a function `finput` and a value `x`, and it returns the function `finput` evaluated at `x` and then squared. # In[15]: def takesquare(finput, x): return finput(x) ** 2 # We can now call `takesquare` with any function $f$ that can be called as $f(x)$ and returns one value. For example, we can call it with the cosine function, and we can test right away whether we got the right answer # In[16]: print('takesquare result:', takesquare(np.cos, 2)) print('correct value is: ', np.cos(2) ** 2) # ### Exercise 6. Numerical integration # Numerical integration of a function is a common engineering task. # The `scipy` package has a specific subpackage called `integrate` with a number of numerical integration functions. We will use the `quad` function. Use the `quad` function to integrate the function $f(x)=\text{e}^{-x}$ from 1 till 5 (note that `quad` returns two values; read the doc string to find out what they are). Check that you did it right by doing the integration by hand (which is easy for this function). # # Next, compute the following integral: # # $$\int_1^5 \frac{\text{e}^{-x}}{x}\text{d}x$$ # # This integral is more difficult to do analytically. Perform the integration numerically with the `quad` function and check your answer, for example, at the [wolframalpha website](https://www.wolframalpha.com) where you can simply type: `integrate exp(-x)/x from 1 to 5`. # In[ ]: # Answer to Exercise 6 # ### Interactive functions using `ipywidgets` # The package `ipywidgets` constains widgets for use in Jupyter Notebooks. We will start here by using the simplest form, which is the `interact` function. The `interact` function can be used to, you guessed it, interact with a function. The `interact` function is a bit slow when interacting with a graph, but other than that it works nicely. # # For example, let's write a function that plots a line with length $L$ that makes an angle $\alpha$ with the horizontal. The angle $\alpha$ is an input argument. # In[17]: def plot_line(alpha): L = 20 x = L / 2 * np.cos(np.deg2rad(alpha)) y = L / 2 * np.sin(np.deg2rad(alpha)) plt.plot([-x, x], [-y, y]) plt.axis('scaled') plt.xlim(-20, 20) plt.ylim(-20, 20) plot_line(45) # The `interact` function of `ipywidgets` can now be used to interact with the `plot_line` function we just defined. The `interact` function takes as input arguments the name of the function to interact with, and then the minimum value, maximum value, and step of the input arguments. When you execute the code below, a slider will appear and you can move the slider to change the value of the angle $\alpha$ (and again, it is a bit slow, so don't move it around too quickly). # In[18]: from ipywidgets import interact interact(plot_line, alpha=(-180, 180, 5)); # The `plot_line` function can be modified to also take the color of the line as an input argument. # In[19]: def plot_line(alpha, color='g'): L = 20 x = L / 2 * np.cos(np.deg2rad(alpha)) y = L / 2 * np.sin(np.deg2rad(alpha)) plt.plot([-x, x], [-y, y], color) plt.axis('scaled') plt.xlim(-20, 20) plt.ylim(-20, 20) plot_line(45) # The `interact` function can be used to both change the value of $\alpha$ and the color of the line. Provide all the possible colors as a list, which will appear as a dropdown box when executing the code. # In[20]: from ipywidgets import interact interact(plot_line, alpha=(-180, 180, 5), color=['orange', 'brown', 'fuchsia']); # ### Exercise 7. First wiget # Write a function that plots $a\cos(x)$ and $b\sin(x)$ on the same graph for $x$ going from $0$ to $4\pi$. Set the limits of the vertical axis from -5 to +5. Input arguments of the function are the amplitudes $a$ and $b$, the color of the cosine function, and the color of the sine function (so 4 input arguments in total). Use the interact function to allow $a$ and $b$ to vary from 0 to 5, the color of the cosine function to be orange, pink or red, and the colors of the sine function to be blue, grey or black. # In[ ]: # Answer to Exercise 7 # ### Answers to the exercises # Answer to Exercise 1 # In[21]: def test(x, alpha): return np.exp(-alpha * x) * np.cos(x) x = np.linspace(0, 10 * np.pi, 100) y1 = test(x, 0.1) # This function can be called with an array y2 = test(x, 0.2) plt.plot(x, y1, label=r'$\alpha$=0.1') plt.plot(x, y2, label=r'$\alpha$=0.2') plt.xlabel('x') plt.ylabel('f(x)') plt.legend(); # Back to Exercise 1 # # Answer to Exercise 2 # In[22]: def psi(x, y, U=1, R=1): r = np.sqrt(x ** 2 + y ** 2) if r < R: rv = 0.0 else: theta = np.arctan2(y, x) rv = U * (r - R ** 2 / r) * np.sin(theta) return rv print(psi(2, 4, U=2, R=1.5)) print(psi(0.5, 0, U=2, R=1.5)) # Back to Exercise 2 # # Answer to Exercise 3 # In[23]: x,y = np.meshgrid(np.linspace(-4, 4, 100), np.linspace(-3, 3, 100)) psivec = np.vectorize(psi) R = 1.5 p = psivec(x, y, U=2, R=R) plt.contour(x, y, p, 50) alpha = np.linspace(0, 2 * np.pi, 100) plt.fill(R * np.cos(alpha), R * np.sin(alpha), ec='g', fc='g') plt.axis('scaled') # Back to Exercise 3 # # Answer to Exercise 4 # In[24]: def velocity(x, y, U=1, R=1): r = np.sqrt(x**2 + y**2) theta = np.arctan2(y, x) if r > R: vr = U * (1 - R**2 / r**2) * np.cos(theta) vt = -U * (1 + R**2 / r**2) * np.sin(theta) vx = vr * np.cos(theta) - vt * np.sin(theta) vy = vr * np.sin(theta) + vt * np.cos(theta) else: vx,vy = 0.0,0.0 return vx,vy print('velocity at (2,3): ', velocity(2, 3, U=2, R=1.5)) x,y = np.meshgrid(np.linspace(-4, 4, 50), np.linspace(-3, 3, 50)) vx, vy = np.zeros((50, 50)), np.zeros((50, 50)) R = 1.5 for i in range(50): for j in range(50): vx[i,j], vy[i,j] = velocity(x[i,j], y[i,j], U=2, R=R) alpha = np.linspace(0, 2 * np.pi, 100) plt.fill(R * np.cos(alpha), R * np.sin(alpha), ec='g', fc='g') plt.streamplot(x, y, vx, vy) plt.axis('scaled') plt.xlim(-4, 4) plt.ylim(-3, 3); # Back to Exercise 4 # # Answer to Exercise 5 # In[25]: def dfuncdx(x): if x < 0: rv = -np.sin(x) else: rv = -np.exp(-x) return rv d = 1e-6 x = -1 dfdx = (func(x+d) - func(x-d)) / (2*d) print('True value ', dfuncdx(x)) print('Approx value ', dfdx) x = 1 dfdx = (func(x+d) - func(x-d)) / (2*d) print('True value ', dfuncdx(x)) print('Approx value ', dfdx) # Back to Exercise 5 # # Answer to Exercise 6 # In[26]: def func1(x): return np.exp(-x) def func2(x): return np.exp(-x) / x from scipy.integrate import quad print('func1:') print('numerical integration:', quad(func1, 1, 5)[0]) print('analytic integration:', -np.exp(-5) + np.exp(-1)) print('func2:') print('numerical integration:', quad(func2, 1, 5)[0]) print('wolframalpha result:', 0.218236) # Back to Exercise 6 # # Answer to Exercise 7 # In[27]: def plot_func2(acos=1, asin=1, colorcos='C0', colorsin='C1'): x = np.linspace(0, 4 * np.pi) plt.plot(x, acos * np.cos(x), colorcos) plt.plot(x, asin * np.sin(x), colorsin) plt.ylim(-5, 5) interact(plot_func2, acos=(0, 5, 0.5), asin=(0, 5, 0.5), colorcos=['orange', 'pink', 'red'], colorsin=['blue', 'grey', 'black']); # Back to Exercise 7