Description
Introduction This lab introduces you to the basics of the Python programming language. We will focus on those aspects of the language that make it ideal for numerical calculations. As you learn, you will want to create and share Python programs for your own research and applications. So it is important that, along with the basic facts of the language, you should also pay attention to programming style, that is, how an experience programmer tries to organize code so that it is clear, readable, easy to fix, modify, improve and share. 2 Python files Python can be used interactively, but this should be regarded as a sort of “scratchpad” option, because you really don’t get a good record of what you have done. Programming requires some care, and often the first few attempts at doing a calculation include mistakes. During interactive use, this would mean the entire calculation might have to be retyped again. A better option for serious work is to use Python’s scripting (programming) facility. With sequences of Python commands contained in files, it is easy to see what calculations were done to produce a certain result, and it is easy to show that the correct values were used in producing a graph. It is terribly embarrassing to produce a very nice plot that you show to your advisor only to discover later that you cannot reproduce it or anything like it for similar conditions or parameters. When the commands are in clear text files, with easily read, well-commented code, you have a very good idea of how a particular result was obtained. And you will be able to reproduce it and similar calculations as often as you please. The Python comment character is a “pound” or “hash” sign (#). That is, lines starting with # are not read as Python commands and can contain any text. Similarly, any text on a line following a # can contain textual comments and not Python commands. It is important to include comments in .py files to explain what is being accomplished. # Th is Python s c r i p t w i l l s im pl y say ” h e l l o ” . print ( ” H ell o ! ’ ) # You can a l s o add comments on the same l i n e a s your commands . Listing 1: An example of comments in a Python script. 3 Variables, numpy, and scipy Python uses variable names to represent data. A variable can represent some important value in a program, or it can represent some sort of dummy or temporary value. Important quantities should be given names longer than a few letters, and the names should indicate the meaning of the quantity. For example, if you were using Python to generate a matrix containing a table of squares of numbers, you might name the table tableOfSquares. (The convention I am using here is that the first part of the variable name should be a noun and it should be lower case. Modifying words follow with upper case letters separating the words. This rule comes from the officially recommended naming of Java variables.) Once you have used a variable name, it is bad practice to re-use it to mean something else. To enable the level of scientific computing needed for the class we will make use of several open-source Python libraries. Two of the most important are the packages numpy and scipy. You can think of the numpy 1 package as containing the base data types and operations which we will need. Some examples include the ability to make vectors and matrices, matrix and vector addition, define commonly used functions such as sin(x). scipy builds on top of numpy and is used for actual numerical code such as methods for solving linear systems. In addition, we will be using a library called matplotlib.pyplot to create plots of our calculations. 4 Exercise 1 Start an interactive Python session. 1. Import the numpy library using the command import numpy a s np 2. What value does the variable np.pi return? 3. Print the variable np.pi to 4 decimal places using the command print ( f ”{np . pi : . 4 f }” ) 4. Do a quick Google search for ”python fstrings” and explain what f”{np.pi:.4f}” is doing. 5. In your own words, describe what the float() function is doing. 6. Request and print the machine epsilon value using these commands: ep s = np . f i n f o ( f l o a t ) . ep s print ( f ”machine ep s = { ep s }” ) 7. Set the variable a=1 and the variable b=1+eps. What is the difference in the way that Python displays these values? Can you tell from the form of the printed value that a and b are different? 8. Set the variable c=2 and the variable d=2+eps. Are the values of c and d different? Explain why or why not. 9. Set the variable x to some value that you choose. 10. What is the square of x? Its cube? In Python, the operator ** is used for taking a variable to a power. For example 2 ** 3 = 8. 11. Define the variable theta, which represents an angle. Assign a value to theta. Choose an angle θ and set the variable theta to its 12. What value do np.sin(theta) and np.cos(theta) return? Angles can be measured in degrees or radians. Which of these has Python used? 13. Python variables can also be given “character” or “string” values. A string is a sequence of letters, numbers, spaces, etc., surrounded by single or double quotes quotes (’) or “. In your own words, what is the difference between the following two expressions? a1=’np.sqrt(4)’ a2=np.sqrt(4) 2 5 Vectors and matrices Thoughout this course we will use vectors and matrices. Vectors and matrices are defined as numpy arrays. Below we provide some basic details about these. The easiest way to define a row vector is the np.array() command with the values listed in square brackets rowVec1 = np . a r r a y ( [ 0 , 1 , 3 , −6, np . pi ] ) Listing 2: Define a row vector of 5 values. A column vector can be described similarly, except there is another layer of square brackets for each individual column entry. columnVector1 = np . a r r a y ( [ [ 0 ] , [ 1 ] , [ 3 ] , [ 6 ] , [ np . pi ] ] ) Listing 3: Define a column vector of 5 values. Matrices can be written by simply adding more entries to the inner square brackets appearing in the columVector1 example above. For instance The matrix A = 1 2 3 4 5 6 7 8 9 (1) can be generated with the expression A = np . a r r a y ( [ [ 1 , 2 , 3 ] , [ 4 , 5 , 6 ] , [ 7 , 8 , 9 ] ] ) Often we will need to generate a vector of equally spaced values, which can be useful for plotting and other tasks. This can be done using the numpy function np.linsapce(). This function has the form: np.linspace( firstValue, lastValue, numberOfValues ) We generate 6 evenly spaced numbers in the interval [10,20] using the command: e ven s = np . l i n s p a c e ( 1 0 , 2 0 , 6 ) The vector evens now contains the values [10, 12, 14, 16, 18, 20]. We could generate fifty evenly-spaced points with e ven s = np . l i n s p a c e ( 1 0 , 2 0 , 50 ) 6 Exercise 2 Begin an interactive session with Python. 1. Import the numpy library: import numpy a s np 2. Use the np.linspace() function to create a numpy array called meshPoints containing exactly 1000 values with values evenly spaced between -1 and 1. meshPoints = np . l i n s p a c e ( −1, 1 , 1000 ) 3 3. What expression will yield the value of the 95th element of meshPoints? Keep in mind that Python arrays start at index 0. What is this value? 4. Use the np.size() function to again confirm that the vector meshPoints that you created has length 1000: np . s i z e ( meshPoints ) 5. Use your meshPoints array to plot a sinusoid on the interval [−1, 1]: import m a t pl o tli b . p y pl o t a s p l t p l t . pl o t ( meshPoints , np . s i n ( 2 ∗ np . pi ∗ meshPoints ) ) p l t . show ( ) 6. In your own words describe what each of the above commands does. 7. Create a file named exercise2.py. The first lines of the file should be the following: # Lab 2 , e x e r c i s e 2 # A sample s c r i p t f i l e . # Your name and t h e d a te Follow the header comments with the commands containing exactly the commands you used in the earlier parts of this exercise. Test your script by typing into the terminal python e x e r c i s e 2 . py Remember your terminal needs to be in the directory containing this script file. 7 Vector and matrix operations The Python packages numpy and scipy provide tools for matrix and vector manipulation. We will investigate a few of these. 8 Exercise 3 In an interactive Python session define the following vectors and matrices: import numpy a s np rowVec1 = np . a r r a y ([ −1 , −4 , −9]) c olVec 1 = np . a r r a y ( [ [ 2 ] , [ 8 ] , [ 9 ] ] ) mat1 = np . a r r a y ( [ [ 1 , 3 , 5 ] , [ 7 , 9 , 0 ] , [ 2 , 4 , 6 ] ] ) 1. You can multiply vectors by constants. Compute c olVec 2 = ( np . pi / 4 ) ∗ c olVec 1 2. The np.cos() function can be applied to a vector to yield a vector of cosine values. Compute c olVec 2 = np . c o s ( c olVec 2 ) Print these new values using the command print ( c olVec 2 ) 4 Note that the original values of colVec2 have been overwritten. Are these the values you expect? 3. You can add vectors. Compute c olVec 3 = c olVec 1 + c olVec 2 and print these values. 4. The Euclidean norm of a matrix or a vector is available using np.linalg.norm(). Compute np . l i n a l g . norm ( c olVec 3 ) 5. You can do row-column matrix multiplication using np.dot(). Compute c ol v e c 4 = np . dot (mat1 , c olVec 1 ) print these values. 6. We can transpose a matrix or vector using the command np.transpose() mat1Transpose = np . t r a n s p o s e ( mat1 ) rowVec2 = np . t r a n s p o s e ( c olVec 3 ) print these values. 7. You can, instead, tranpose matrices and vectors by applying the .T method, as follows: mat1Transpose = mat1 .T rowVec2 = c olVec 3 .T 8. Matrix operations such as determinant and trace are available. de te rmin an t = np . l i n a l g . de t ( mat1 ) t r = np . t r a c e ( mat1 ) Print these values. 9. You can pick certain elements out of a vector. Use the following command to find the smallest element in rowVec1. np .min( rowVec1 ) or by using the numpy array .min method: rowVec1 .min( ) 10. For matrices the np.min() and np.max() functions can be set to work along one dimension at a time. Therefore they will produce vectors rather than a single scalar if desired. np .max(mat1 , a x i s =0) np .max(mat1 , a x i s =1) In your own words describe what the axis option does. 11. You can compose a sequence of operations in a single statement. For example, use the following expression to compute the max norm of a vector. np .max( np . abs ( rowVec1 ) ) 5 9 Flow control It is critical to be able to ask questions and to perform repetitive calculations in Python programs. These topics are examples of “flow control” constructs in programming languages. Python provides two basic looping (repetition) constructs: for and while, and the if construct for asking questions. These statements each precede an indented block of Python statements to be repeated as necessary. Python loops are very flexible; generally in this course we will be looping over integers and making use of the range function. In other languages there might be a statement which signifies the end of the loop/if statement, such as the end statement in Matlab. However, in Python, white space and indentation have an actual meaning. In the case of loops and if statements whatever is immediately indented afterwards is considered a part of the loop/if statement. The end of indented statements signifies the end of the loop. The range() command produces a sequence of numbers starting from 0 by default and incrementing by 1 by default. These defaults can be modified. The range function is inclusive of the first value, but stops before the specified end value. Syntax Example for loop for control-variable in sequence: Python statement . . . . . . nFactorial=1; for i in range(1,n+1): nFactorial=nFactorial*i while loop Python statement initializing a control variable while logical condition involving the control variable: Python statement . . . . . . Python statement changing the control variable nFactorial=1 i=1; # initialize i while i <= n: nFactorial=nFactorial*i i=i+1 a simple if if logical condition: Python statement . . . . . . if x ~= 0: # ~ means “not” y = 1/x a compound if if logical condition: Python statement . . . . . . elif: logical condition . . . else: . . . if x ~= 0: y=1/x elif np.sign(x) > 0: y = 1 else: y = -1 10 Exercise 4 The trapezoid rule for the approximate value of the integral of e x on the interval [0, 1] can be written as Z 1 0 e x dx ≈ h 2 e x0 + h N X−1 k=1 e xk + h 2 e xN where h = 1/N and xk = 0, h, 2h, . . . , 1. The following Python code uses the trapezoid rule, with N=40, to carry out this approximation. Use cut-and-paste, or download the file exercise4.py. import numpy a s np 6 # Use t h e t r a p e z o i d r u l e t o approx ima te t h e i n t e g r a l from 0 t o 1 # o f exp ( x ) , u s i n g N i n t e r v a l s # Your name and t h e d a te N = 40 h = 1 / N x = 0. 0 a p p r o x I n t e g r al = 0 fo r k in range ( 0 , N+1) : i f k == 0 or k == N: a p p r o x I n t e g r al = a p p r o x I n t e g r al + ( h / 2 ) ∗ np . exp ( x ) e l s e : a p p r o x I n t e g r al = a p p r o x I n t e g r al + h ∗ np . exp ( x ) x = x + h print ( f ”The approximate i n t e g r a l i s { a p p r o x I n t e g r al }” ) 1. Cut and paste this code into a Python file exer4.py. Run the code and report the outputted values. Is the final value for approxIntegral nearly equal to e 1 − e 0 ? 2. What is the complete sequence of all values taken on by the variable x? 3. How many times is the following statement executed? a p p r o x I n t e g r al=a p p r o x I n t e g r al + ( h / 2 ) ∗np . exp ( x ) 4. How many times is the following statement executed? a p p r o x I n t e g r al=a p p r o x I n t e g r al + h ∗ np . exp ( x ) 11 Scripts, functions and graphics If you have to type everything at the command line, you will not get very far. You need some sort of scripting capability to save the trouble of typing, to make editing easier, and to provide a record of what you have done. You also need the capability of making functions or your scripts will become too long to understand. In this section we will consider first a pure script file and later introduce functions. We will be using graphics in the script file, so you can pick up how graphics can be used in our work. The Fourier series for the function y = x 3 on the interval −1 ≤ x ≤ 1 is y = 2X∞ k=1 (−1)k+1 π 2 k − 6 k 3 sin kx. (2) Computationally, we can only add up a finite number of terms in this series. We can compare the resulting finite Fourier approximation to the original function x 3 . 12 Exercise 5 Use cut-and-paste, or download the file exercise5.py and then answer the questions about the code. 7 import numpy a s np import m a t pl o tli b . p y pl o t a s p l t # compute NTERMS terms o f t h e F ou r ie r S e r i e s f o r y=x ˆ3 # p l o t t h e r e s u l t u s i n g NPOINTS p o i n t s from −1 t o 1 . # Your name and t h e d a te NTERMS = 20 NPOINTS = 1000 x = np . l i n s p a c e ( −1, 1 , NPOINTS) y = np . z e r o s ( np . s i z e ( x ) ) fo r k in range ( 1 , NTERMS + 1 ) : term = 2 ∗ (−1) ∗ ∗( k + 1 ) ∗ ( np . pi ∗∗ 2 / k − 6 / k ∗∗ 3 ) ∗ np . s i n ( k ∗ x ) y = y + term p l t . pl o t ( x , y , ”b” ) # ’ b ’ i s f o r b l u e l i n e p l t . pl o t ( x , x ∗∗ 3 , ”g ” ) # ’ g ’ i s f o r a green l i n e p l t . show ( ) It is always good programming practice to define constants symbolically at the beginning of a program and then to use the symbols within the program. Sometimes these special constants are called “magic numbers.” By convention, symbolic constants are named using all upper case. 1. Add your name and the date to the comments at the beginning of the file. 2. How is the Python variable x related to the dummy variable x in Equation (2)? (Please use no more than one sentence for the answer.) 3. How is the Python statement that begins y=y… inside the loop related to the summation in Equation (2)? (Please use no more than one sentence for the answer.) 4. In your own words, what does the line y = np.zeros(np.size(x)) do? Hint: You can google np.zeros() and np.size() to find the documentation for these functions. 5. Execute the script by typing python exercise5.py at the command line. You should see a plot of two lines, one representing a partial sum of the series and the green line a plot of x 3 , the limit of the partial sums. You do not have to send me this plot. 13 Exercise 6 In this exercise you are going to modify the calculation so that it continues to add terms so long as the largest component of the next term remains larger in absolute value than, say, 0.05. Since the series is of alternating sign, this quantity is a legitimate estimate of the difference between the partial sum and the limit. The while statement is designed for this case. It would be used in the following way, replacing the for statement. TOLERANCE= 0 . 0 5 ; # t h e chosen t o l e r a n c e v al u e <> k=0 term = TOLERANCE + 1 # b i g g e r than TOLERANCE while ( np .max( abs ( term ) ) > TOLERANCE) : k = k + 1 <> print ( f ”Number o f i t e r a t i o n s = {k}” ) <> 8 1. Copy the file exercise5.py to a new file called exercise6.py and change the comments at the beginning of the file to reflect the objective of this exercise. 2. Modify the program by replacing the for loop with a while loop as outlined above. 3. What is the purpose of the statement term= TOLERANCE + 1 4. What is the purpose of the statement k = k + 1 5. Try the script to see how it works. How many iterations are required? Does it generate a plot similar to the one from the previous exercise using a for loop? Hint: If you try this script and it does not quit (it stays “busy”) you can interrupt the calculation by holding the Control key (“CTRL”) and pressing “C”. If you find that your code does not quit, you have a bug. For some reason, the value of the variable term is not getting small. Look at your code carefully to find out why. If you cannot see it, use the debugger to watch execution and see what is happening to the value of term. 14 Exercise 7 The next task is to make a function out of exercise6.py in such a way that it takes as argument the desired tolerance and returns the number of iterations required to meet that tolerance. 1. Copy the file exercise6.py to a new file called exercise7.py. We will turn it into a function by using the line below def e x e r c i s e 7 ( t o l e r a n c e ) : # comments # Your name and t h e d a te All code within the function must be indented. 2. Replace the upper-case TOLERANCE with a lower-case tolerance because it no longer is a constant, and throw away the line giving it a value. 3. Add comments just after the signature and usage lines to indicate what the function does. 4. Delete the lines that create the plot and print so that the function does its work silently. 5. Add the statement return k at the end of the function. This tells Python that when the function is called it should return the value k. 6. Start up an interaction session and type the command from e x e r c i s e 7 import e x e r c i s e 7 this line is saying from the file exercise7.py load the function exercise7() into the Python environment. Next, invoke this function from the command line by choosing a tolerance and calling the function. Using the command 9 num ItsRequired=e x e r c i s e 7 ( 0 . 0 5 ) 7. Print this value. How many iterations are required for a tolerance of 0.05? This value should agree with the value you saw in Exercise 6. 8. To observe convergence, how many iterations are required for tolerances of 0.1, 0.05, 0.025, and 0.0125? Python allows a function to return two or more values instead of a single value. The syntax to accomplish this trick is simple: we just add it to the return statement with a comma. So if our function computes two values, y and ]ttz, we return both by a statement like: return y , z For a function named “funct” that returns two variables depending on a single variable as input, the function call would look like: y , z = f u n c t ( x ) If you have the same function but wish only the first output variable, y, you would write y , = f u n c t ( x ) and if you wish only the second output variable, z, you would write , z = f u n c t ( x ) In either case, the omitted output variable is simply ignored. 15 Exercise 8 1. Copy the file exercise7.py to a new file called exercise8.py. Modify the function so that it returns first the converged value of the sum (a vector) and, second, the number of iterations required. Be sure to change the comments to include a description of all input and output parameters. 2. What form of the command line would return only the (vector) value of the sum to a tolerance of 0.03? 3. What is the norm of this vector to 14 digits of accuracy. Hint: The command print(f’{a:.10f}’) would print out the variable a to 10 digits of accuracy. 4. What form of the command line would return the (vector) value of the sum and the number of iterations required to achieve a tolerance of 0.02? How many iterations were taken, and what is the norm of the (vector) value of the sum, to 14 digits of accuracy? 16 Exercise 9 We will often find it useful in this course to have function names as variables. For example, the sine function in exercise8 could be replaced with an arbitrary function. While the series might not converge for some choices of functions, it would converge for others 1. Copy the file exercise8.py to a new file called exercise9.py. Modify the function so that it accepts a second parameter called func: def e x e r c i s e 9 ( t ol e r a n c e , func ) : 10 and replace the np.sin() function inside the sum with func(). Do not forget to modify the comments in the file to reflect this change. 2. Test that exercise9() and exercise8() give equivalent results when func is really sin by executing the following commands: import numpy a s np from e x e r c i s e 8 import e x e r c i s e 8 from e x e r c i s e 9 import e x e r c i s e 9 y8 , k8=e x e r c i s e 8 ( 0 . 0 2 ) ; y9 , k9=e x e r c i s e 9 ( 0 . 0 2 , np . s i n ) ; #t h e f o l l o w i n g d i f f e r e n c e s h o ul d be sm all np . l i n a l g . norm ( y8−y9 ) 3. Use exercise9 to compute the (vector) sum of the series for exercise9(.02, np.cos) and plot the result. Please include this plot with your summary. You have now completed all the required work for this lab. Send it to me in .tar.gz or .zip file. It should include your summary file, each of the files exercise2.py, exercise5.py, exercise6.py, exercise7.py, exercise8.py and the plot files you created. 11
