This chapter provides a focused yet comprehensive introduction to Python programming, with a strong emphasis on scientific and engineering applications. It begins with the fundamentals—numbers, variables, logical operations, strings, and user input—then progresses to essential data structures such as lists, tuples, dictionaries, and sets. Core programming concepts like control flow (using loops and conditionals), modular code design with functions, file input/output, and error handling are covered in detail, laying a solid foundation for more advanced tasks. Students are gradually introduced (accompanied by the classwork and lab assignments provided in Part 2) to key scientific libraries that form the backbone of Python-based computation: NumPy for efficient numerical operations and array handling, Matplotlib for creating a wide range of visualizations, and SciPy for specialized scientific routines. The chapter concludes with a brief introduction to object-oriented programming, demonstrating how to structure code in a clean, reusable manner using classes and objects.

Interactive code demonstrations are provided through Trinket.io, allowing students to experiment directly within the browser without needing local setup. By the end of the chapter, students will be prepared to write clear, efficient Python code and apply it confidently to real-world scientific problems.

 Although students are free to use any Python IDE or code editor throughout the course, the classwork and lab assignments (in Part 2) are specifically designed to be completed in a Unix environment using the command line. This approach provides valuable hands-on experience with scripting and execution in a terminal-based workflow, which is common in scientific and high-performance computing contexts. By developing and running Python programs through the command line, students not only deepen their understanding of Python syntax and logic but also build essential skills in navigating Unix systems, managing files, and using shell commands—foundational tools for reproducible and scalable scientific computing. 

2.1 Python Basics

Python is a powerful and flexible programming language developed by Guido van Rossum in 1989. Designed as an interpreted, high-level, object-oriented language with dynamic typing, Python is especially well-suited for rapid development and problem-solving. It is widely adopted in scientific computing and data analysis, making it an excellent tool for implementing and exploring numerical methods. Python is also a core component in many Linux distributions, often used for system administration. For more resources, visit the official Python website: www.python.org.

Why Python?

• Python is an open-source programming language, meaning its source code is freely available and can be used, modified, and distributed at no cost. This accessibility has contributed to its widespread adoption across academia and industry.
• One of Python’s key strengths is its clean and readable syntax, which makes it relatively easy to learn—especially for those new to programming. This simplicity allows students and researchers to focus more on problem-solving and algorithm development rather than getting bogged down in complex language details.
• Python is cross-platform, meaning it runs on virtually all major operating systems including Windows, Linux, and macOS. Programs written in Python are typically portable across these systems with little or no modification, enabling seamless collaboration and deployment in diverse computing environments.
• A major advantage of Python for numerical methods is its rich ecosystem of libraries. The Python Standard Library provides a broad set of built-in modules, and additional tools essential for scientific computing—such as NumPy (numerical arrays), SciPy (scientific algorithms), and Matplotlib (plotting and visualization)—are freely available and widely supported by the community.
• Python is also a flexible, multi-paradigm language, supporting procedural, object-oriented, and functional programming styles. This makes it adaptable to various problem-solving approaches and suitable for building both simple scripts and complex software systems.
• Although Python is generally slower than fully compiled languages like C or Fortran, for many numerical applications the difference in execution  time is outweighed by the faster development cycle, improved readability, and ease of debugging. In practice, Python often serves as a high-level interface to optimized, low-level routines written in C, C++, or Fortran—combining productivity with performance.

To begin, you need to install Python, or verify that you already have it installed. Here I describe how to perform some basic tasks using the command-line interface in a Unix or Unix-like operating system (such as Linux or macOS). Windows users can achieve the same by securely logging into a school Unix account using SSH clients like MobaXterm, PuTTY, or the Windows Subsystem for Linux (WSL). Among these, MobaXterm offers a user-friendly interface with built-in support for terminal emulation, remote file editing, and graphical forwarding, making it a convenient choice for beginners and advanced users alike.

To start a Python shell, open a terminal, type “python” and press Enter.

$ python

You will see system message and the command prompt  >>> as:

Python 3.6.4 |Anaconda custom (64-bit)| (default, Jan 16 2018, 18:10:19) [GCC 7.2.0] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>>

Since we will be working with Python 3^[1], make sure the Python version numbers reported on the first line is python 3.x.y where the minor revisions x,y are not important.

To exit the Python shell, type exit().

>>> exit()

Python operates in two primary modes: normal mode and interactive mode, each serving different purposes in the programming workflow.

Interactive Mode

Python shell provides immediate output for every executed statement. It also stores the output of previously executed statements in active memory. When the new statements are executed by the Python interpreter, the entire sequence of previously executed statements is considered while evaluating the current output. Interactive mode is good for testing small code snippets or debugging.

In a Python shell, try:

$ python

>>> print('Hello, World!’)

Press Enter, you will see the output:

Hello, World! 

In an IPython console, try:

$ ipython

In [1]: print('Hello, World!') 

Press Enter, you will see the output:

Hello, World! 

Normal Mode

Normal mode (also called script mode) is used when running Python programs saved in files, typically with a .py extension. In this mode, the entire script is executed from top to bottom, and the output is displayed after the script runs. This is the standard approach for developing larger programs or automating tasks, where the code is written, saved, and then run as a complete unit.

Create a file with filename ‘hello.py’  and write the statement ‘print(“Hello, World!”)’ to the file. Save the file and run it (from the same directory) with the Python interpreter as:

$ python hello.py

Edit your file hello.py inserting the line ‘#!/user/bin/env python’. Now, you can also run your python program as:

$ chmod u+x hello.py  (What does it do?)

$ ./hello.py

Note: Python is installed in all the computers in this classroom and Math Computer Lab. To install Python in your laptop, I recommend the Anaconda Distribution which include the NumPy, SciPy and Matplotlib libraries. It also comes with Spyder, a popular Integrated Development Environment (IDE) with Matlab like interface. 

Trinket.io 

Trinket.io is a web-based coding platform that supports multiple programing languages, including Python. It uses primarily normal/script mode, but it also provides a web-based interface that can mimic some interactivity during runtime.  Trinket eliminates the need for any local software installation, allowing students to write, run, and modify code seamlessly from any device with internet access, making it ideal for educational settings that prioritize ease of access and quick setup. With Trinket, students can explore algorithm behavior, create convergence plots, and experiment with code in real time—all without setup or compatibility concerns.

Throughout this chapter, we will use Trinket to demonstrate interactive programming directly in the browser. A trinket.io console consists of two main sections: a code editor, where users type their code (e.g., in Python), and an output window, where results, print statements, or error messages appear after execution. Edit the print statement and run it by clicking the arrow button.

Note: There are several excellent platforms that offer true interactive Python experiences, including Jupyter Notebooks, JupyterLab, Google Colab, PythonAnywhere, Deepnote, and Replit.

2.2 Numbers and Variables

Python has three number types: integer (int), floating point numbers (float) , complex numbers (complex).

1. Integers are whole numbers like 1, -11 etc. Integer arithmetic is exact, limitation is restricted only by computer’s memory.

2. Floating point numbers are the representation of real numbers like 1.234 x 10⁵ in computer (not always exact). They are stored in binary with certain precision (15-16 decimal places). A single number containing a dot (’.’) is considered a float. Scientific notation is supported by e or E. 123400., 1.234e5, 0.01234, 1234e-5

3. A complex number is specified by either adding its real part to an imaginary parts (denoted by j suffix) as in 3.0+2.1j, or complex(3.0,2.1).

Try these in an interactive Python shell, Python console,  or Trinket console:

type(10); type(3.5); type(2+3j); int(-3.5); float(10); complex(2,3)

Arithmetic Operators in Python

+        addition. (e.g., 2+3=5)

-         subtraction (e.g., 5-3=2)

*     multiplication (e.g., 2*3=6)

/          floating point division (e.g., 6/3=2)

//       integer division (e.g., 6/2=3, 7//2=3)

%          modulus (remainder)  (e.g., 7%2=1)

**       exponentiation (e.g., 2**3=8, 3**2=9)

Note:  Exponentiation (**) has the highest precedence, followed by multiplication and division (*, /, //, %) and addition and subtraction (+, -). Operators of equal precedence are evaluated left to right with the exception of exponentiation which is evaluated right to left.

• Try these in an interactive mode: 

        Try: 6.5//2; 6.5%2; 8/4/2; 2**2**3; 3*(4/3-1)-1 (what did you get? Was it expected?)

• There are two built-in mathematical functions that do not need to import math library, namely abs and round.

Try: abs(-3.14); abs(2+3j); round(5.54); round(9.5); round(4.5)

• Real/imaginary parts and conjugate of complex numbers are accessed using dot (.) method.

Try: (2+3j).real; (2+3j).imag; (2+3j).conjugate()

• There are many useful mathematical functions in the ‘math module’, which can be accessed using "import math".

Try: import math; math.sqrt(4); math.exp(1); math.log(2); math.log(100,10); math.cos(0); math.atan(1)

• Use import cmath for the complex functions.

Try:  import cmath; cmath.sqrt(-1); cmath.exp(complex(0,1))

• Two useful constants in math module are  math.pi and math.e.

Try:  import math; math.sin(math.pi/2); math.log(math.e)

• You can also import math module using “from math import *" and access all the available functions directly without using “math.function_name“.

Try: from math import * ; sin(pi/2) ; log(e)

For a complete list of the mathematical functions in the math module visit: http://docs.python.org/3/library/math.html

Variables

When an object, such as a number, is created in Python code, memory is allocated for it. The location of this memory is called its ‘address’. It is not convenient to refer to the object by its address, we need a variable name or object identifier.  Here are some rules for a valid variable name.

• Variable names are case sensitive.
• Variable names can contain any letter, underscore, any digits, but cannot start with a digit.
• The following Python reserved keywords or built-in constants cannot be used for variable names: and, assert, break, class, continue, def, del, elif, else, except, False, finally, for, from, global, if, import, in, is, lambda, None, nonlocal, not, or, pass, print, raise, return, True, try, while, yield

Note: Variable names should be meaningful, but not be too long; do not use, I, l, O (similar to 1, 0); use i,j,k for integer counters; use ’max_score’ instead of ’MaxScore’ (reserved for class names).

2.3 Comparisons and Logic

Python Comparison operators

==        equal to

!=         not equal to

>          greater than

<          less than

>=       greater or equal to

<=        less or equal to

• The result of a comparison is a boolean object (type: bool) which is either True or False. 
• Comparison can be modified and combined together with the logic operator keywords not, and, or.

Truth table for not, and and or

2.3 Strings

A Python string object (type: str) is an ordered, immutable sequence of characters.

2.3.1 Defining a string object

To define a variable containing a string of characters, the string should be enclosed by the either single or double quotes:

s1 = "Hello, Sir! How are you doing?"


s2 = "" (an empty string)

Try: print(s1); print(s2)

Strings can be concatenated using either + (plus) operator or by placing them next to each other on the same line.

s3 = "input"+".dat" (or s3 = "input" ".dat")

Strings can be repeated using * operator.

s4 = "abc"*3
s5 = ("aa"+"bb")*2

The built-in function, str converts an object passed as its argument into a string.

s6 = str(1e-3)
Try: print(s6), print(4*s6)

Python has no restriction on the length of a line. For readability, 79 characters is recommended. To break up a string over two or more lines, use the line continuation character, ‘’ or (better) enclose the string in parenthesis.

s7 = "This is an example of" "a long string!"
s8 = ("This is another example of" 
      "a long string!")

The choice of quotes for the strings allows us to include the quote character inside the string: just define it using the other quote:

s9 = "It doesn’t make sense!"

If you need to include both quotes in the string, you should use: \’ for single quotes and \" for double quotes:

s10 = "She said, "She isn’t worried about Python.""

Also try:  \n (linefeed), and \t (horizontal tab).

s11 = "Physics\nChemistry\nMathematics\nBiology" 
s12 = "Physics\tChemistry\tMathematics\tBiology"

All of these snippets are implemented for you to try below.

2.3.2 Indexing and slicing strings

Indexing a string in Python return a single character at a given position, the first character having index 0 and final character in a string of n characters has index n - 1. Python raises an IndexError if you try to index a string out of its length.

Set s = "Embry-Riddle" and print s[0], s[6], s[11], s[-1], s[-6], s[12]

Slicing a string, s[i:j], gives a substring between the characters at two indexes, including the first (i) and excluding the second (j). If the first index is omitted, 0 is assumed; if the second is omitted, the string is sliced to its end.

Try: s[:5], s[1:6], s[6:], s[:]

The optional, third number k in a slice s[i:j:k] specifies the stride. If omitted, default is k = 1. Negative value of k reverses the string.

Try: s[::2], s[1:6:2], s[-1:-6:-1]

  String methods

Python strings are immutable objects. It is not possible to change a string by assignment. Try:

s="orange" 
s[0]="O"
print(s)

s1=s[1:] 
s2="O"+s1
print(s2)

The built-in function len gives the number of characters in a string. Try:

len(s)

String objects come with a large number of methods which are accessed using the dot notation. Try:

s= "EmbrY_RiddLe aeronutical universiTy"

s.upper(), s.lower(), s.title(), s.title().replace("_", "-").center(72)

2.3.3 User input

The function input can be used to obtain information from the user. The simplest use is to assign a string to a variable: s = input()

Try these:

1. test1.py

Note that the input function produces a string value. You cannot use them to do arithmetic operations.

2. test2.py 

3. test3.py 

 4. test4.py 

1. 1. 1. 1. 2. test5.py (Note: Here you should enter only integer values!)

Methods on Strings

• Python strings are immutable objects; it is not possible to change a string by assignment. Try:

s="orange" 
s[0]="O"
print(s)

s1=s[1:] 
s2="O"+s1
print(s2)

• The built-in function len gives the number of characters in a string. Try:

len(s)

• String objects come with a large number of methods which are accessed using the dot notation. Try:

s= "EmbrY_RiddLe aeronutical universiTy"

s.upper(), s.lower(), s.title(), s.title().replace("_", "-").center(72)

1. 2.3.4    String Formatting 

Consider the following code that uses Python’s str.format() method to insert values into a string

print("{} plus {} = {}".format(3, 4, "seven"))
print("{1} - {0} = {2}".format(2, 7, 5))

Python string format provides a powerful way to format the integers, with specifier ‘d’. Try:

Some useful format specifiers for the floating point numbers are: ‘f’ (fixed-point notation), ‘e’ (exponential or scientific notation), and ‘g’ (a general format that uses scientific notation for very small and very big numbers). The desired precision (number of decimal places) is specified as ‘.p’ after the minimum field width. Try:

The print function

• The built-in function print in Python takes a list of objects, and optional arguments end and sep that respectively specify which characters should end the string and which characters should be used to separate the printed objects. Try:

Trinket.io Uses Python2, so most print usage as a function will not work properly!

2.4 Lists, Tuples and Loops

2.4.1 Lists

• Python list is an ordered, mutable array of objects. A list is constructed by specifying the objects, separated by commas, and enclosed in square brackets.

list0 = []                                    (an empty list)
list1 = [0, ’one’, 2.0, -3, 4/5]              (a list of five numbers)

• Python list can contain references to any type of object: strings, various types of numbers, built-in constants, and even other lists.

list2 = [’zero’,1, 3.14, 2+3j, True, list1]

• An item can be retrieved from the list by indexing it (index starts from zero). Try:

list1[0], list1[-1], list2[3], list2[5], list2[5][1], list2[5][2], list2[0][1]

• Python lists are mutable object. i.e., the items can be reassigned and altered. Try:

list2[0] = 0.0; list2

• List can be sliced in the same way as string sequences.

list3 = list2[1:4]; list4 = [1::2]; list5 = [::-1]

List methods

Python lists come with a large number of powerful methods. Since list objects are mutable, they can change shapes without having to copy the contents to a new object, as in the strings. Some useful methods are listed below.

• append(element) : add an item to the end of the list

  list2.append(6)

• extend(list) : add one or more objects by copying them from another list

  list2.extend([11,12,13])

• insert(index,element) : insert an item at the specified index

  list2.insert(1,100)

• pop() : remove and return the last element from the list

  list5.pop()

• remove(element) : remove a specified item from the list

  list2; list2.remove(1)

• reverse() : reverse the list in place

  list1; list1.reverse(); list1

• sort() : sort the list in place

list1.sort(); list1

• sorted(list) : returns a new sorted list, keeping the original unaltered, default is ascending order

sorted(list4); list4; sorted(list4, reverse=True)

• Try the following construction of the data structure known as stack:

• The string method, split generates a list of substrings from a given string, split on a specified separator (default is “white space”).

s=’Sun,Mon,Tue,Wed,Thu,Fri,Sat’
l=s.split(’,’); print(l)

2.4.2 Tuples

• Python tuple objects can be thought of immutable list. Tuples are created by placing items inside parentheses separated by commas.

t0=()                                                                                       an empty tuple !
t1=(’one’,)                                                     a singleton tuple, note the trailing comma !
t2=(’one’, 2, [’a’, ’b’, ’c’], ’four’, 5.0)

• Tuples can be indexed or sliced the same way as lists.

t2[0]; t2[2][2]; t2[1:3]

• As tuples are immutable, they cannot be appended or extended, or have some elements removed. But we can alter the elements within the list inside the tuples. i.e.,

t2[0]=’two’; t2[2]=[’d’, ’e’, ’f’]                                    not possible !
t2[2][0]=’d’; t2[2][1]=’e’; t2[2][2]=’f’                                  possible !

• Tuple packing/unpacking, swapping

t3 = 1, 2, 3                                t3=(1,2,3)

a, b, c = 1, 2, 3                        a=1, b=2, c=3

a, b = b, a                              t=a, a=b, b=t

• Also try the built-in functions:    list(’string’); tuple([1,2,’three’])

2.4.3 for loops

• If you want to loop over the elements in a list, you can use a simple for loop construction, for item in list : to loop over the elements in list, set to item each step.

Math = [’MA305’, ’MA432’, ’MA438’, ’MA448’, ’MA453’, ’MA484’]
for course in Math:        #note the colon (:) at the end 
    print(course)           #note the indentation

• Python has a useful method range of referring to a sequence of numbers to loop over. It can be constructed with up to three arguments: range([a0=0], n, [stride=1]), where the first integer a0 and the stride (which can be negative) are optional with default values 0 and 1 respectively.

n=len(Math)
for i in range (n): 
    print(i, ’:’, Math[i])

• The object created by range can be indexed, cast into lists and tuple and iterated over. Try:

range(5)[2]; range(1,6)[3]; list(range(0,6,2)); list(range(10,0,-2)); tuple(range(5))

• enumerate takes an iterate object and produces, for each item in turn, a tuple (count, item ), consisting of a counting index and the item itself.

Math = [’MA305’, ’MA432’, ’MA438’, ’MA448’, ’MA453’, ’MA484’]
for i, course in enumerate(Math): 
    print(i, ’:’, course)

• Python’s built-in function zip creates an iterate object in which each item is a tuple of items taken in turn from the sequences passed to it. This allows us to iterate over many sequences simultaneously.

list1=[1, 2, 3, 4]
list2=[’a’,’b’,’c’,’d’]
for pair in zip(list1,list2): 
    print(pair)

• The function zip can also be used to unzip sequences of tuples.

list1=[0, 1, 2]
list2=[1, 2, 3] 
list3=[’a’,’b’,’c’] 
z=zip(list1,list2,list3) 
for item in zip(*z):
    print(item)

2.5 Control Flow

2.5.1 if/then/else statements

• Python controls the decisions that a computer can make in a similar way as other programming languages, through if/then/else statements and loop structures with specified termination criteria. Levels of control flow are determined by the indentation level of the code. Four spaces are recommended.

if logical_expression_1: #parenthesis not needed, must end with : 
    statements_1         #must be indented, 4 spaces recommended
elif logical_expression_2: 
    statements_2
else:
    statements_3

• The if statements use standard comparison expressions such as <, >, <=, >=, ==, and ! = to test a truth value. They also incorporate more natural language such as: is, is not, and, or.
• A value is True unless it is 0 (int), 0.0 (float), empty string, ”, empty list, [], empty tuple, (), or None (a special value). Try:

for x in range(10):
    if x % 2: # if x%2=1 (True), x is odd 
        print(x, ’is odd!’)
    else:     # else x%2=0 (False), x is even 
        print(x, ’is even!’)

2.5.2 while loops

• The while loop is used to execute statements as long as some condition holds.

i=0 
while(i<10):
    i +=1
    print(’Hello ’, i) 
print()

• Try the following example that implements Euclid’s algorithm for the greatest common divisor of a and b:

a, b = 112, 24
print(’a=’, a, ’b=’, b)
while b: # the loop is executed till b=0 (’False’) 
    a, b = b, a%b # 1) temp=a%b, 2) a=b, 3) b=temp
print(’greatest common divisor of a and b:’, a)

2.5.3 break/continue/else

• The command break issued inside a loop immediately ends that loop and moves execution to the statements following the loop.

i=0
while True:
    i += 1
    if not (x %15 or x%25): break
print(x, ’is divisible by both 15 and 25’)

• The command continue immediately forces the next iteration of the loop without completing the statement block for the current iteration.

for i in range(1,11): 
    if i % 2:
        continue 
    print(i, ’is even!’)

• A for or while loop may be followed by an else block of statements, which will be executed only if the loop is terminated normally.

list1=[ 0, 1, 3, -4, 5, 6]
for i, item in enumrate(list1): 
    if item < 0:
        break
    else:   # useful if list 1 has no negative items 
        print(item , ’occurs at index’, i)

Exercises

sum_gs=0
for i in range(n): 
sum_gs += 1/2**i 
print(i,sum_gs)

sum_hs=0
for i in range(1,n+1): 
    sum_hs += 1/i 
    print(i,sum_hs)

sum_as=0
for i in range(1,n+1):
    sum_hs += (-1)**(n+1)*(1/i) 
    print(i,sum_hs)

i=0 
sum_gs=0 
while True:
    sum_gs += 1/(2**i) 
    print(i,sum_gs) 
    err=abs(2-sum_gs) 
    if err < 5e-8:
    break
    i += 1

2.6 Dictionaries and Sets

Dictionaries

A dictionary in Python is a type of “associative array” that can contain any objects as its values, but unlike lists and tuples, in which the items are indexed by an integer starting with 0, each item in a dictionary is indexed by a unique key, which may be any immutable object.

• The dictionary can be defined by giving key:value pairs between the curly braces:

score = {'English': 80, 'Mathematics': 95, 'Physics':72, 'Chemistry':92} 
print(score['Mathematics'])
course='Chemistry' 
print(score[course])

• An alternating way of defining a dictionary is to pass a sequence of (key, value) pairs to the dict constructor. If the keys are simple strings, the pairs can also be specified as keyword arguments to this constructor.

score=dict([('English',80), ('Mathematics',95), ('Physics',72), ('Chemistry',92)] 
score=dict(English=80, Mathematics=95,Physics=72,Chemistry=92)

• A for-loop iteration over a dictionary returns the dictionary keys (in no particular order).

for subject in score: 
     print(subject,score[subject])

• The methods keys, values and items return a dictionary’s keys, values and key-value pairs (as tuples) respectively.

print(score.keys()) 
print(score.values()) 
print(score.items()) 

Sets

A set is an unordered collection of unique items. The sets are useful for removing duplicates from a sequence and finding the union, intersection and difference between two collections.

• A set is created y listing its elements between braces ({…}) or by passing an iterable to the set() constructor.

A=set([1,2,3,,a,b,c]) 
B=set([1,1,2,3,2,3,3,4,,a,b,c,b])
print(A) 
print(B)

• As the sets are unordered, the set objects cannot be indexed or sliced, but they can be iterated over, tested for membership, and they support the built it function len.

print( len(A), len(B) ) 
print( 3 in A, x in A ) 
print( A <= B )
print( A | B ) 
print( A & B ) 
print( B - A )

2.7 Functions and Modules

A function in Python is a way of writing reusable code similar to that found in many other programing languages.

Syntax of a function

def name_of_function (arg1, arg2, arg3):
    value = some_expression_with_arg1_arg2_arg3 
    return value

The def statement defines a function name_of_function listing the arguments (if any) arg1, arg2, arg3, and ends with a colon (:). The function’s statements are written in an indented block following this def. The return statement returns the specified value.

Here is an example of defining a function f(x) = 1 + 1/2(e^x + e^-x).

def f(x):
fx = 1+(exp(x)+exp(-x))/2 
return fx

Once defined, the function f can be called any number of times from any place. Try:

a=1 
fa=f(a)
print(’f(’,a,’)’=fa)
print(’f(0)’=f(0))

Python functions can return any type of object using return statement. Functions without the return statement return Python’s special value None. Python functions can also return multiple values by packing them into a tuple.

def roots_quad (a,b,c):
    """Returns the roots of ax^2+bx+c=0.""" 
    d = b**2 -4*a*c
    if d >= 0:
        x1=(-b+sqrt(d))/(2*a)
        x2=(-b-sqrt(d))/(2*a) 
    else:
        x1=complex(-b, sqrt(-d))/(2*a)
        x2=complex(-b, -sqrt(-d))/(2*a) 
    return x1,x2

You can call the above function as:

r1,r2=roots_quad (1,-2,1)
print(’roots of {0: 2d}x^2 + {1: 2d}x + {2: 2d} = 0:n
                {3:0.7g}, {4:0.7g}’.format(a,b,c,r1,r2) )

• Function definitions can appear anywhere in a Python program, but a function cannot be referenced before it is defined.
• Functions can be nested, but a function defined inside another is not (directly) accessible from outside that function.
• The docstring appearing as the first statement of a function within a triple quoted string is the special  __doc__  attribute of the function. It should provide information about how to use the function. Try:

roots_quad.__doc__  

• Python functions can have variable identifiers assigned to them, they can be passed as arguments to other functions, and they can even be returned from other functions.
• It is possible to assign more than one variable name to the same function. Try:

rq=roots_quad 
rq(1,-2,1)

2.7.1 Keyword and Default Arguments

Functions with multiple arguments can be called by passing the arguments in the same order (positional) or by setting them explicitly as keyword arguments. Try:

roots_quad(a=1, c=1, b=-2) 
roots_quad(b=-2, a=1, c=1) 
roots_quad(1, c=1, b=-2)

When mixed, positional argument must come first, the following call is not allowed.

roots_quad(b=-2, 1, 1)

Python functions may also set variable defaults. If the caller does not provide a value for this (optional) argument, a default value is used. The default arguments are set in the function definition.

def newton_method(f, x0, Tol=1.e-5)
    ....
    r = ....
    return r

You can call this function with the first two arguments only or all three arguments with a different value of Tol.

x=newton_method(f, 1)
y=newton_method(f, 1, Tol=1.e-12)

2.7.2 Local and Global Variables

Local variables defined inside a function are not available outside, the variables defined outside the function are available everywhere.

def test_f():
    x=1             # x is a local variable
    print(x,y)      # y is nonlocal, the function can access it if assigned outside

Try:

y=2 
test_f() 
x=3 
test_f() 
print(x,y)

Note that there can be local and global variable with the same name. If you want to modify the global variables from function calls, you need to declare the variables with the keyword global.

def test_g():
    global y 
    y += 1
    print(y)

Try:

y=2 
test_g() 
print(y)

Python’s rule for resolving scope: local → enclosed (nested function) → global → built-ins.

2.7.3 Recursive Functions

A function that can call itself is called a recursive function.

Factorial of a nonnegative integer n:

[latex]n!=\left\{

		\begin{array}{ll}

		n\cdot(n-1)!,  & n\ge 1\\

		1,  & n=0\\

		\end{array}\right.[/latex]

Fibonacci sequence:

[latex]\left\{

		\begin{array}{l}

		\lambda_0=1,\, \, \lambda_1=1\\

		\lambda_n=\lambda_{n-1}+\lambda_{n-2}, \, \,\,\,  n \ge 2

		\end{array}

		\right.

[/latex]

Factorial function

Factorial function can be implemented as:

def factorial ( n ): # assumes n is nonnegative integer 
    if n <= 1:
        factn = 1 
    else:
        factn = n*factorial(n-1) 
    return factn

Fibonacci Sequence

Fibonacci sequence can be generated by appending to a list as:

n=15
fib = [1,1]
for i in range(2,n+1): 
    fib.append(fib[i-1]+fib[i-2])
print(fib)

Here is an example of the use of recursive function:

def fibonacci ( n ): 
    if n <= 1:
       return 1 
    else:
       return fibonacci(n-1)+fibonacci(n-2)
n=15
for i in range(n+1):
    print(’’, fibonacci(i), end=’’)

Classwork

1. Write a function to compute the double factorial of a positive integer.

   [latex]

   $$

   n!! = \left\{

   	\begin{array}{ll}

   	\displaystyle{\prod_{i=1}^{(n+1)/2}} (2i-1) = 1.3.5.7.\cdots.(2n-1) & (n\,\,\mathrm{odd})\\[9pt]

   	\displaystyle{\,\, \, \,\prod_{i=1}^{n/2}} 2i = 2.4.6.8.\cdots.2n & (n\,\,\mathrm{ even})

   	\end{array}\right.    $$

   [/latex]

2. Write a function to determine f(n), the number of trailing zeros in n!, using the special case of de Polignac’s formula:

[latex]\[ f(n) = \sum_{i}^{n}\lfloor \frac{n}{5^i} \rfloor, \][/latex]

where [x] denotes the floor of x, the largest integer less than or equal to x.

2.7.4 lambda Functions

A lambda function in Python is a type of simple in-line function. The executable body of a lambda function must be an expression and not a statement. That is, it may not contain loop blocks, conditionals or print statements.

Examples

1. f(x) = 3x³+ 2x²+ x

f=lambda x: 3*x**3 + 2*x**2 + x 
print(f(2)

    2. f(x, y) = x³+ y³+ xy

f=lambda x,y: x**3 + y**3 + x*y 
print(f(1,2)

2.7.5 Modules

Python has a large library of modules and packages that extends its functionality. Many of these are available as part of the “standard library” provided with the Python interpreter itself. NumPy, SciPy and Matplotlib are very popular packages used in Scientifc Computing.

• The import statement makes reference to the modules that are Python files containing definitions and statements.
• If you write import <module> in a program, the Python interpreter executes the statements in the file <module>.py and enters the module name into the current namespace, so that the functions/procedures it defines are available with the “dotted syntax”: <module>.<function>. Recall how you used math module:

import math                     ( use math.sqrt(2), math.pi )
from math import sqrt, pi       ( use sqrt(2), pi )
from math import *              ( What does this do? )

• You can also define your own module by placing a code within a file module_name.py somewhere in your computer (for small projects, usually in the same directory where you run your code). Note: DO NOT name your module anything that isn’t a valid Python identifer (hyphen, starting with a digit, name of a built-in functions, modules, etc.)

Classwork: Revisit the Bisection code

Let us organize the Bisection code (from Lab 4) using your own modules.

• First, download the code lab4.py from your course Canvas, and put it in a new directory CW7 and make three copies of the file.

$ cp lab4.py bisect.py
$ cp lab4.py funct.py
$ cp lab4.py main.py

• In the main.py delete line 1 to 52.

$ vi main.py
:8,52d

• In the funct.py delete line 17 to end of file.

$ vi funct.py
:17,$d

• In the bisect.py delete line 8 to 17 and line 53 to end of the file.

$ vi bisect.py
:8,17d
:53,$d

• In main.py insert the following lines after the second line.

from funct import F
from bisect import bisection

• Run your code:

$ ./main.py               ( make sure the file has executable permission! )

Also try:

1. Make a copy of main.py to cw7a.py and make the following change.

import funct as mf
from bisect import bisection as bb

and run your new code. What happened? How can you make it run?

    2. In funct.py write the something inside the doc_string. For example,

""" This is my module to keep functions: 
F(x)=x**3+x**2-3.0*x-3.0, g(x)=1-2*x """

and in cw7a.py write the following line and run it.

    print(mf.__name__) 
    print(mf.__doc__) 
    print(bb.__name__) 
    print(bb.__doc__)

What did you see?

• In cw7a.py replace the line if __name__ ==’ __main__’:          by the following line def main(): #if __name__ == ’__main__’: and insert main() in the last line of the code and run it. What did you see?

2.8 File Input/Output

2.8.1 Opening and Closing a File

• You can open a file with a given filename and mode. Filenames can be given as their absolute path or path relative to the directory from which the program is run. Mode is a string with values r (read), w (write), a (append), r+ (read and write) and rb, wb, ab, rb, rb+.

f = open(’out.txt’,’w’) # opens a file ’out.txt’ for text-mode writing

• The write method of a file object writes a string to the file and returns the number of characters written. Run write1.py:

f.write(’This is a test of writing in a output file!’)

• The file objects are closed with the close method:

f.close()

• The built-in function print takes an argument, file, to specify where to redirect its output.
• The following script writes a multiplication table for 1 to 9 in a file “mult_table.txt”. Run write2.py:

f = open(’mult_table.txt’,’w’) 
for i in range(1,10):
    # print , separate by list
    # print(i, 2*i, 3*i, 4*i, 5*i, 6*i, 7*i,8*i, 9*i, sep=’,’, file=f) 
    # print tab separated list
    # print(i, 2*i, 3*i, 4*i, 5*i, 6*i, 7*i,8*i, 9*i, sep=’t’, file=f) 
    # print formatted list
    print(’{0:2d} {1:2d} {2:2d} {3:2d} {4:2d} {5:2d} {6:2d} {7:2d}
           {8:2d}’.format(i,2*i,3*i,4*i,5*i,6*i,7*i,8*i,9*i),file=f)
f.close()

Reading from a file

• f.readline() reads a single line (one at a time) from the file, up to and including the newline character. f.readline() returns an empty string when it reaches the end of file.

• Try the following script read1.py.

  # opens a file ’mult_table.txt’ for text-mode reading 
  f = open(’mult_table.txt’,’r’)
  s1 = f.readline() 
  print(s1)
  s2 = f.readline() 
  print(s2)

• file objects are iterable, and looping over a text file returns its lines one at a time. Try read2.py:

f = open(’mult_table.txt’,’r’) #opens ’mult_table.txt’ for reading 
for line in f:
    print(line, end=’’)        #print(line) produces extra blank line

• The following script read3.py reads data file ‘mult_table.txt’ and prints the multiplication table for 3, 6 and 9 on screen.

f = open(’mult_table.txt’,’r’) 
for i in range(9):
    line =f.readline() 
    Line=line.split() 
    print(Line[2],Line[5],Line[8])
f.close()

Examples

• The following script read4.py reads the data file ‘mult_table.txt’, stores the values of the table for 4 and 8 in the lists four and eight, and prints on screen nicely.

f = open(’mult_table.txt’,’r’) #opens a file ’mult_table.txt’ for reading 
four, eight = [], []
for line in f.readlines(): 
    Line=line.split()
    four.append(int(Line[3])) 
    eight.append(int(Line[7]))
f.close()
for i in range(1,10):
print(’4 x’,i,’=’, four[i-1], ’t 8 x’,i,’=’, eight[i-1])

• The following script read5.py uses f.readline() for the same.

f = open(’mult_table.txt’,’r’) 
four, eight = [], []
for i in range(9):
    line = f.readline() # line is a string one line in the file f 
    Line = line.split() # Line is a list 
    four.append(int(Line[3]))
    eight.append(int(Line[7])) 
f.close()
for i in range(1,10):
print(’4 x’,i,’=’, four[i-1], ’t 8 x’,i,’=’, eight[i-1])

2.8.2 Redirecting Input/Output

Reading stdin (standard input) from a file

• The standard input for any Unix command can be redirected to come from a file (instead of the keyboard) simply using the command < filename.
• In Python, we need to import sys module for this. Let’s try it. Run the script read6.py.

#!/usr/bin/env python3 
import sys

s1 = sys.stdin.readline() 
print(s1)

• Set execution permission to read6.py and run it with the following command.

$ chmod u+x read6.py
$ ./read6.py < mult_table.txt

What did you see?

• How it is different from s1 = input(’Enter a string:’) ?

Sending stdout (standard output) to a file

• Similarly, the standard output (screen) of any Unix command can be redirected to a file with the command > filename
• Let’s try it.

  $ ./read6.py < mult_table.txt > out.txt

  Take a look in the output file out.txt. What happened?

• Also try this:

    $ rm out.txt
  
  $ ./read6.py > out.txt

  What happens? Why? What should you do to make it work?

• And:

$ ./read1.py > out.1txt
$ ./read2.py > out2.txt
$ ./read3.py > out3.txt
$ ./read4.py > out4.txt
$ ./read5.py > out4.txt

Example

Save the following script in a file named test.py. Set execution permission to test.py and run it.

#!/usr/bin/env python3
import sys

#dat = input(’Enter the input data file:’) 
#f = open(dat,’r’)

f = sys.stdin
four, eight = [], [] 
for i in range(9):
line = f.readline() # line is a string one line in the file f 
Line = line.split() # Line is a list 
four.append(int(Line[3]))
eight.append(int(Line[7])) 
f.close()
for i in range(1,10):
print(’4 x’,i,’=’, four[i-1], ’t 8 x’,i,’=’, eight[i-1])

2.9 Errors and Exceptions

There are two types of error in Python: syntax errors (error in the grammar of the language) and exceptions (runtime errors).

• Syntax errors are checked for before the program is executed and caught by the Python compiler and produce a message indicating where the error occurred.
• Exceptions are caused by attempting an invalid operation on an item of data. There are several built-in and user defned mechanism for catching this type of errors and handling the condition gracefully without stopping the programs’s execution.

Some built-in exceptions

• FileNotFoundError: File cannot open, the file does not exist.

• IOError: File cannot open, input/output system error.

• IndexError: Subscript of a list or string is out of range.

• KeyError: Key does not available in the dictionary.

• NameError: Cannot reference a variable, not defined.

• TypeError: Inappropriate type as an argument to a function.

• ValueError: Correct type but incompatible value in function.

• ZeroDivisionError: Attempting to divide by zero.

• SystemExit: Raised by the sys.exit function.

Handling exceptions

• To catch an exception in a block of code, you should write the code within a try: clause and handle any exceptions raised in an except: clause.

try :
    y=1./x
    print(’y=’,y)
except (zeroDivisionError):
    print(’division by zero encountered, x=’,x) 
except (NameError):
    print(’x is not defined’)

Raising exceptions

• The keyword raise allows a program to force a specific exception and customize the message or other data associated with it. For example:

if n%2:
    raise ValueError(’n must be even!’) 
    # continue statement with n even.

• The keyword assert evaluates a conditional expression and raises an AssertionError exception if that expression is not True. For example:

def dotproduct(x,y):
      assert len(x)==len(y),’x, y must be of the same size’
      # continue statement with x, y are of same size
      dp=0
      for i in range(len(x)):
           dp +=x[i]*y[i]
      return dp

2.10 Differences between Python 2 and Python 3

Python 3 was released in 2008 and it is not backward compatible. Python 3 has taken away many redundant features from Python2 2 and simplified some other features and added many new features in it. Due to the extensive documentation, books and a large collection of third party libraries, there are still many users of Python 2. However, most organizations have already started migrating their codes in from Python 2 to Python 3 and almost all new projects in Python are using Python 3. That it why, it is said, “Python 2.x is legacy, Python 3.x is the present and future of the language.”

Here are some noticeable differences between Python 2 and Python 3.

• The print statement in Python 2 is replaced by print() function in Python 3.
• In Python 2, the result of division of two integers is an integer. In Python 3, the result of division of two integers is a float.
• The xrange() function from Python 2 is removed in Python 3. The function range() is used in Python 3.
• In Python 3, the exception arguments should be enclosed in parentheses. In Python 2, it was optional.
• While handling exceptions in Python 3, the keyword as before the the parameter to handle arguments is a must. In Python 2, it was not needed.
• Python 3 does not support old style classes.
• In Python 3, strings are Unicode by default.

2.11 Introduction to NumPy

NumPy is the fundamental package for scientific computing with Python and is the base for many other packages.

• • NumPy has a powerful custom N-dimensional array object for efficient and convenient representation of data.
  • NumPy supports vectorization: a single operation can be carried out on an entire array, no need to loop over the array’s element!
  • The absence of explicit looping and indexing makes the code cleaner, less error-prone and closer to the standard mathematical notation it reflects.
  • NumPy has tools for integration with other programming languages used for scientific programming like C/C++ and Fortran.
  • Moreover, it will introduce you to a whole other set of libraries such as SciPy and Scikit-learn, which you can use to solve almost any problem in scientific computing and data science.

2.11.1 Basic Array Methods

The most powerful construct of NumPy is the ndarray.

• The NumPy array class ndarray provides a generic container for multidimensional array of a homogeneous data which can be sorted, reshaped, written to and read from files and much more.
• That is, unlike Python lists and tuple, the element in NumPy array cannot be of different types. The type in the ndarray is specified by (optional) data type object (dtype).
• Import NumPy package with the statement import numpy as np and then refer to its attributes with the prefix np. Try:

import numpy as np 
x=np.array( [1,2,3] )
y=np.array( [[1, 2, 3],[4,5,6]]] ) 
z=np.array([[[1,2],[3,4]],[[5,6],[7,8]]])

print(x) 
print(y) 
print(z)
print(x[0],x[2], x[-1]) 
print(y[0,1],y[1,2],z[0,0,0],z[1,1,1])

print(type(y)) 
print(y.shape) 
print(y.ndim) 
print(y.dtype) 
print(y.size)

• x.ndim returns the dimension of the array x.
• x.shape returns a tuple giving the number of elements in each dimension of the array x. For example, the output (3,) means an array of 3 elements, (2,3) means 2 rows with 3 columns, (2,2,2) means 2 rows of arrays of 2 rows and 2 columns.
• x.size returns the total number of elements in the array x.
• x.dtype returns the data type of the elements in the array x. The data type is “upcast” to the most general type if they are mixed but compatible types. You can also set the data type in the np.array construct. Try:

import numpy as np
x=np.array( [1,2.,3] ) 
print(x) 
print(x.dtype)
x=np.array( [1, 2, 3+4j] ) 
print(x)
print(x.dtype)

x=np.array( [1,2,3],dtype=complex)
# float,float32,float16,int,int8,int16,int32 
print(x.dtype)
print(x.nbytes)

2.11.2 Creating/Initializing an Array

• np.arange(n) creates a list of 0 to n–1 (like range, but it can generate floating point sequence!)
• np.linspace(a,b,n) creates a sequence (1D array) of n evenly spaced numbers between a and b. np.linspace has a couple of optional boolean arguments.
• Setting retstep=True, returns the step-size. endpoint=False omits the final point in the sequence. Try:

import numpy as np
x=np.arange(10) ; y=np.arange(10,dtype=float)
print(x); print(y)

x=np.linspace(-np.pi,np.pi,20)
# creates 20 evenly spaced numbers between -pi and pi 
x,dx=np.linspace(-np.pi,np.pi,20,retstep=True)
# optional retstep=True returns the step-size 
print(x); print(dx) 
y=np.linspace(0,np.pi,10,endpoint=False) 
print(y)

• Also try: np.geomspace(a,b,n) and np.logspace(a,b,n)

import numpy as np 
u=np.geomspace(1,64,5) 
print(u) 
v=np.logspace(3,4,5) 
print(v)

• You can create an array of a desired shape using np.empty. The initial element values are undefined.
• np.zeros and np.ones create an array of specified shape with elements pre filled with 0 and 1 respectively.
• np.empty, np.zeros and np.ones also take optional dtype argument. The default data type is float.
• You can create new arrays of the same shape of the exiting arrays, use np.empty_like, np.zeros_like or np.ones_like. Try:

import numpy as np
x=np.empty((2,3)) # not initialized
y=np.zeros((3,2)) # prefilled with 0’s 
z=np.ones((3,3)) # prefilled with 1’s 
print(x); print(y); print(z) 
x=np.ones((2,3), dtype=int)
print(x)

a=np.ones_like(y) 
print(a)

b=np.arange(6).reshape(2,3) 
print(b)

2.11.3 Array Operations

NumPy provides many familiar mathematical functions that the math module does, implemented as universal functions that act on each element of the array producing an array in return.

• Thus, no need to loop explicitly over the array elements. Try these operations in one dimensional arrays, print each statement:

import numpy as np 
x=np.linspace(-5,5,4) 
xp1=x+1
xt2=x*2 
xm1=x-1 
x2=x**2
sqrtx=np.sqrt(x) 
sinx=np.sin(x) 
ex=np.exp(x)

y=np.zeros(4)+2 
xmy=x-y
xoy=x/y 
xty=x*y

• Also try the following:

import numpy as np 
M=np.identity(4)
# creates an identity matrix of size 4x4 
N=M+3
P=2*M 
Q=M-N 
R=M-x

2.11.4 Reading and Writing Arrays to a File

In scientific computing, we need to read data frequently from a text file. NumPy provides several functions for reading data from a text flie.

• np.save saves NumPy arrays in (platform-independent) binary file format (.npy) that can be read (on any other operating system) using the function np.load. Try:

M=np.random.rand(100,3) 
np.save(’array1.npy’,M)
# saves the array M in a file array1.npy 
N=np.load(’array1.npy’)
# loads the array from array1.npy as N 
print(M); print(N)

• np.savetxt saves the data in human readable text-files. np.loadtxt reads the data from text-files. Try:

np.savetxt(’array2.txt’,M) 
P=np.loadtxt(’array2.txt’) 
print(P)

• np.loadtxt has several optional arguments to handle the data file that contains descriptive header and columns are aligned in a fixed width format or separated by one or more delimiting characters (spaces, tabs or commas). Syntax:

np.loadtxt(fname, dtype=<class ’float’>, comments=’#’, delimiter=None, 
           skiprows=0, usecols=None, unpack=False)

The description of the arguments of np.loadtxt

• fname: only required argument, name of the file to read data from

• dtype: data type of the array, default is float

• comments: comments in a file usually (# or %) to be ignored

• delimiter: string used to separate columns in a file, default is None, meaning any amount of whitespace (spaces, tabs) delimits the data. To read comma-separated (csv) file, set delimiter=’,’
• skiprows: an integer giving the number of lines at the start of the file to skip over before reading the data, default is 0 (no header)
• usecols: a sequence of column indexes determining which columns of the file to return as data, default is None (all columns)
• unpack: by default data is returned as a single array of rows and columns, unpack=True will transpose this array so that individual columns can be picked o˙ an assigned to different variables.

Revisit Classwork 4 after NumPy

In CW4 you learned how to read data from a file using f.readline(s) and sys.stdin.readline(s). Here we use np.loadtxt for the same.

• First, go to the working directory MA305/CLASSWORK/CW4/. Make a copy of cw4a.py as cw4np.py. Comment out all the lines related to reading from the file.
• Import NumPy package with the statement import numpy as np.
• The file ’cw4a.txt’ has data in four columns, and you want 1st, 3rd and 4th columns as D, X and Y. Also, you do not need to save the first line (header)! Call the NumPy function np.loadtxt as:

  D,X,Y=np.loadtxt(’cw5.txt’,skiprows=1,usecols=(0,2,3),unpack=True) 
  print()

• Print D, X and Y to check if it read the data. What did you observe?
• There are two complains: n is not defned, D is not integer. Do the following:

D,X,Y=np.loadtxt(’cw4a.txt’,skiprows=1,usecols=(0,2,3),unpack=True) 
n=len(D)
#change the data type of D float (default) to integer 
D=D.astype(int)

2.11.5 Some Statistical Methods

First, import NumPy using import numpy as np.

• np.sum returns the sum of all elements in the array. np.min and np.max return minimum and maximum values, np.argmin and np.argmax return location of min and max. Try:

u=np.array([0,1,-2, 4,6,5]) 
np.sum(u); np.mean(u)
np.amin(u); np.argmin(u)
np.amax(u); np.argmax(u)

• np.percentile returns a specified percentile, q, of the data along an axis (or along a fattened version of the array if no axis is given). Try:

x=np.array([[0,1,2,3,4,5],[6,7,8,9,10,11]])
np.percentile(x,50)      #along the whole elements of x 
np.percentile(x,50,axis=0)     #across rows of 2D array 
np.percentile(x,50,axis=1)  #across columns of 2D array 
np.percentile(x,90,axis=1)

• np.mean and np.median return arithmetic mean and median of the array an axis (or along a fattened version of the array if no axis is given). Try:

np.mean(x); np.median(x)        #along the whole elements x 
np.mean(x,axis=1); np.median(x,axis=1) #across columns of x

• np.std returns the sample standard deviation

                [insert latex equation here]

  The value of δ is passed thru the optional argument ddof, default value is δ = 0. Try:

  x=np.array([[0,1,2,3,4,5],[6,7,8,9,10,11]])
  np.std(x) 
  np.std(x,ddof=1) 
  np.std(x,axis=1) 
  np.std(x,axis=1,ddof=1)

• np.var calculates the variance σ² of the values in the array.

np.var(x) 
np.var(x,ddof=1) 
np.var(x,axis=1)

• np.histogram creates a histogram from the values in an array which can be plotted with pylab.hist. Try:

import pylab 
marks=np.array([45,68,90,56,23,60,87,75,59,63,72,82,98,47,76])
np.histogram(marks,bins=5,range=100,density=True) 
pylab.hist(marks,bins=5,range=100,density=True) 
pylab.show()

2.11.6 Polynomials

The polynomial of degree n in NumPy is recognized as

p_n(x) = a₀ + a₁x + a₂x² + · · · a_nxⁿ

In order to work with a polynomial, we need to import numpy.polynomial.Polynomial as

from numpy.polynomial import Polynomial

For example, to define the cubic polynomial p(x) = 3 – x – 3x² + x³, we do the following.

from numpy.polynomial import Polynomial

p=Polynomial([3,-1,-3,1]) 
print(p)

Let us try the following.

• Print the coefficients of the polynomial.

print(p.coef)

• Evaluate the polynomial at x = 2.

print(p(2))

• Find the roots of the polynomial.

print(p.roots())

• Find the derivative p′(x), and evaluate it at x = 1.

print(p.deriv())
print(p.deriv(1))

• Find the integral of p(x).

print(p.integ())

• Evaluate the function f(x) = [insert latex equation here]

f=p.integ(lbnd=1) 
print(f(5))

• Find a polynomial whose roots are x = 1, 2, 3.

q=Polynomial.fromroots([1,2,3]) 
print(q)

2.11.7 Working with Matrices

• To find the transpose of a matrix M, use M.transpose(), (or M.T).

import numpy as np
M=np.linspace(1,12,12).reshape(3,4)
M.tanspose()
M.T

Note: Transposing an one dimensional array returns the array unchanged.

• For matrix multiplication in NumPy, you have to use dot() method instead of *. Try:

M*M 
M.dot(M)

• Two dimensional arrays have an index for each axis. if you want to select every item along a particular axis, replace its index by a single colon (:). Try:

M[3,1]
M[2,:]
M[:,1] 
M[1:-1,1:]

• Also Try:

N=M.flatten() 
P=N.resize(3,3) 
a=np.array([1,1,1])
b=np.array([2,2,2])
c=np.array([3,3,3])
A=np.vstack((a,b,c))
B=np.hstack((a,b,c))
C=np.dstack((a,b,c))

Classwork

1. Create a 4 x 5 matrix, a NumPy array of shape (4,5), defined by A[i, j] = i+j+1. Find B = A^T and print A and B.

   • Method 1

     A=np.zeros(4,5) 
     for i in range(4):
         for j in range(5):
             A[i,j]=i+j+1 
     print(’Matrix A:n’,A)
     
     B=np.zeros(5,4)
     for i in range(5):
         for j in range(4):
             B[i,j]=A[j,i]
     print(’Matrix A-transpose:n’,B)

   • Method 2

     def f(i,j):
         return i+j+1 
     A=np.fromfunction(f,(4,5)) 
     B=A.transpose()

   • Method 3

     A=np.fromfunction(lambda i, j: i+j+1, (4,5)) 
     B=A.T

2. Write a function matrix_multipy() that takes two matrices A and B, check if they are compatible for matrix multiplication, and if compatible returns the product matrix C = AB.

   def matrix_multiply(A,B): 
      l=len(A)
      m=len(A[0])
      n=len(B[0])
      # check compatibility condition for matrix multiplication if m != len(B):
      print(’Multiplication is not possible!’) 
      return
   
      C=np.zeros((l,n)) 
      for i in range(l):
          for j in range(n):
              for k in range(m):
                  C[i,j] += A[i,k]*B[k,j]
   
   return C

3. Test your function with the A and B computed earlier.

   A=np.fromfunction(lambda i, j: i+j+1, (4,5)) 
   B=A.transpose()
   C=matrix_multipy(A,B)
   print(A) 
   print(B) 
   print(C)

4. Also try:

C=A.dot(B); print(C)

2.11.8 Linear Systems

Let A be an n x n coefficient matrix and b be an n x  1 right hand side vector (column matrix). The linear system Ax = b can be solve using the NumPy function numpy.linalg.solve() as follows.

• Import NumPy and set the dimension n of the coefficient matrix A.

  import numpy as np
  
  n=5 
  A=np.zeros((n,n))
  b=np.zeros((n,1))

• Define the n × n coefficient matrix A.

  for i in range(n):
          A[i,i] = 2
          if i>0:
              A[i-1,i] = -1
              A[i,i-1] = -1

• Define the n × 1 right hand side vector b.

  b=A.dot(np.ones((n,1)))

• Solve the system Ax = b. What do you expect the solution x be?

  x=np.linalg.solve(A,b) 
  print(’x=’,x.T)

• Also compute the norm of the residual Ax − b.

res=np.linalg.norm(A*x - b) 
print(’residual=’,res)

2.12 Matplotlib

Matplotlib is the most popular Python package for producing graphical plots that can be embedded in applications, displayed on screen or output as high-quality image files for publication.

• Matplotlib’s is a large package with full fledged object-oriented interface. The matplotlib.pyplot module provides a convenient way of visualizing data.
• To use Matplotlib, import both matplotlib.pyplot and numpy as

  import matplotlib.pyplot as plt 
  import numpy as np

• Let us start with the simplest (x, y) line plot.

  x=[0,1,2,3,4,5,6,7,8,9]
  y=[0,1,4,9,16,25,36,49,64,81]
  plt.plot(x,y)
  # plt.plot(x,y,’-’)
  # plt.plot(x,y,’--’)
  # plt.plot(x,y,’:’)
  # plt.plot(x,y,’o-’) plt.show()

• To plot (x, y) as a scatter plot:

x=np.random.rand(100) 
y=np.random.rand(100) 
plt.scatter(x,y) 
plt.show()

Classwork

Let us plot the following three functions for −2 ≤ x ≤ 2.

f (x) = x³ + x² − 3x − 3, g(x) = 3x² + 2x − 3, h(x) = 6x + 2

• First, define these functions? For simple function you can also use Python’s in-line lambda function as:

f = lambda x: x**3 + x**2 - 3*x -3

or,

def f(x):
    y= x**3 + x**2 - 3*x -3 
    return y

• Then create a list of values x_i, and evaluate the function y_i = f (xi)?

  X=[]
  Y=[]
  dx=4/100
  for i in range(100): 
      x_i=-2+i*dx 
      X.append(xi) 
      Y.append(f(x_i))

• Plot the data (x_i, y_i).

import matplotlib.pyplot as plt 
plot(X,Y)
plt.show()

Using NumPy

import matplotlib.pyplot as plt 
import numpy as np

x=np.linspace(-2,2,100) 
y1=x**3+x**2-3*x-3 
y2=3*x**2+2*x-3 
y3=3*x+2
fig=plt.figure()
plot(x,y1,’-’,x,y2,’--’,x,y3,’:’)
plt.legend([’f’, ’g’, ’h’],loc=’best’) 
plt.xlabel(’x’)
plt.ylabel(’y’)
plt.title(’My plot with Matplotlib’) 
plt.show()
fig.savefig(’plot1.pdf’) # .png, .eps

Classwork

1. Now plot the data from a file data7.txt.

   x,y,z=np.loadtxt(’data7.txt’,skiprows=3,unpack=True) 
   fig=plt.figure()
   plt.plot(x,y,’o-’,x,z,’d-’)
   plt.legend([’men’, ’women’],loc=’best’) 
   plt.xlabel(’Year’)
   plt.ylabel(’Age’)
   plt.title(’Median age at first marriage in the US, 1890-2010’) 
   plt.show()
   fig.savefig(’plot2.pdf’)

2. The normalized Gaussian function with mean µ and standard deviation σ:

   [insert latex equation here]

a. Write a program to calculate and plot the Gaussian functions with µ = 0 and the three values σ = 0.5, 1.0, 1.5. Use a grid of 1,000 points in the interval −10 ≤ x ≤ 10.

b. Verify (by direct summation) that the functions are normalized with area 1.

c. Finally, calculate the first derivative of these functions on the same grid using the following first order forward difference approximation:

g′(x) ≈ g(x+h) – g(x) / h

   for some suitably chosen, small h.

   Note: Using NumPy, it is possible to do this exercise without using a single (Python) loop.

2.13 Object-oriented Programming in Python

Structured programming can be divided into two categories: procedural and object-oriented. The programs we have seen so far in this course are procedural in nature where we define functions that are called, passed data, and which returns values from their calculations. The functions do not hold their own data or remember their state in between they are called, and we do not modify them once they are defined.

The other programming paradigm that is popular through C++ and Java is object-oriented programming where an object represents a concept of some type which holds data about itself (attributes) and defines functions (methods) to act on the data. The action of the function may cause a change in the object’s state (i.e., it may alter the object’s attributes). An object is created (instantiated ) from a detailed plan of something called a class, which dictates its behavior by defining its attributes and methods.

• Python is an object-oriented programming language. In fact, everything in Python is an object. For example, a string in Python is an object that is an instance of the str class. A str object keeps its own data (bunch of characters forming the string) and provides a number of method to manipulate the data.

s="c++ java python" 
s.replace(’ ’, ’, ’) 
s.capitalize() 
s.upper()
s.lower()
s.title() 
s[0]

• The object-oriented programming style is very useful for larger projects where you need to break a problem down into units of data and operations that is appropriate to carry out on that data.

2.13.1 Revisit the Bisection code again

Now, let’s consider the bisection code developed in section [].  Here we rewrite the same code in an object-oriented way.

Suppose, we have to deal with lots of particles and instead of the plot for the particle path we also want to visualize the particles as an animated plot. Then, a natural object-oriented approach of thinking would be to define a Particle class, with attributes such as x, y (coordinates of the position) and ω (angular velocity). We define another class ParticleSimulator where the evolve method will be described that will be responsible for changing the position of the particles obeying the laws of motion (Euler’s method as done in the above code). The particles can be plotted as points using matplotlib.pyplot.plot function as before. We can animate the evolution of particles over time using the matplotlib.animation.FuncAnimation class.

class BisectionSolver:
def __init__(self, func, a, b, tolerance):
self.func = func
self.a = a
self.b = b
self.tolerance = tolerance
self.n = 0
self.error = abs(b – a)
self.root = None

def solve(self):
Fa = self.func(self.a)
Fb = self.func(self.b)

print(‘ n r F(r) error’)

while True:
r = (self.a + self.b) / 2
Fr = self.func(r)

if Fr == 0:
print(f’Got lucky! root = {r}, F(r) = {Fr}’)
print(f’ in {self.n} tries, error <= {self.error / 2}’)
self.root = r
break

if Fa * Fr < 0:
self.b = r
# Fb = Fr (can be reused later if needed)
else:
self.a = r
Fa = Fr

self.error = abs(self.b – self.a)
print(f'{self.n:3d} {r:15.8g} {Fr:15.8g} {self.error:15.8g}’)
self.n += 1

if self.error < self.tolerance:
self.root = r
break

return self.root

##########################################################
if __name__==’__main__’:
import math

def example_func(x):
return x**3 – x – 2

solver = BisectionSolver(example_func, a=1, b=2, tolerance=1e-6)
root = solver.solve()
print(“Approximate root found:”, root)

2.13.2 Revisit the Hailstone Sequence

#!/usr/bin/env python3
import math
import numpy as np
import matplotlib.pyplot as plt

“””
========================================================================
MA305 – cw 7: your name – date
Purpose: To show that every hailstone sequence starting from any positive
integer n converges to 1. This is the famous Collatz conjecture.
========================================================================
“””

class HailstoneSequence:
def __init__(self,N):
self.N=N
self.maxit=10000000

def compute_and_plot(self):
for n in range(1,self.N):
plt.clf()
it=0
n0=n
print()
self.hs=[n0,]
while n != 1:
it += 1
print(n, end=’,’)
if n==1:
break
elif n%2:
n = 3*n+1
else:
n = n//2
self.hs.append(n)
if it > self.maxit:
print(‘The hailstone sequence does not converge for n=’,n0)
break
print(‘\nThe hailstone sequence converged to 1 in’, it, ‘terms.’)
print()
print(len(self.hs))

print(self.hs)
self.x=np.arange(len(self.hs))
plt.plot(self.x,self.hs,’o-‘)
for i in range(len(self.x)):
#plt.text(self.x[i]+0.1,self.hs[i]+0.1,str(self.x[i]+1),fontsize=8)
plt.text(self.x[i]+0.1,self.hs[i]+0.1,str(self.hs[i]),fontsize=8)
plt.title(‘Hailstone sequence of {0:} converges in {1:} steps’.format(n0,len(self.x)))
plt.pause(1)

def main():
sim = HailstoneSequence(100)
sim.compute_and_plot()
plt.show()

if __name__==”__main__”:
main()

2.13.3 Particle Simulation using Euler’s method for ODEs

Let’s consider the problem of designing a computer program to simulate particles positions (x_i(t), y_i(t)), i = 0, 1, · · ·, n which constantly rotates around a central point at various angular speeds ω_i. The law of motion obeys the following equations.

[insert latex equation here]                    (2.1)

Using the Euler’s method, we can approximate the solutions of the differential equations in (2.1), the position of the particles at t = t + Δt as

[insert latex equation here]                    (2.2)

Thus, for a single particle with the known initial position (x(0), y(0)) the particle path can be obtained by

[insert latex equation here]                     (2.3)

where t_n+1 = t_n + Δt for n = 0, 1, 2, · · · Python program to implement the above numerical procedure is shown below.

1  #!/usr/bin/env python3
2  import numpy as np
3  import matplotlib.pyplot as plt
4
5  def evolve (x0, y0, t):
6      n=len(t)
7      x=np.zeros(n); y=np.zeros(n)
8      dt=(t[−1]−t[0])/n
9      x[0]=x0; y[0]=y0
10     for i in range (n−1):
11         norm=np.sqrt(x[i]∗∗2+ y[i]∗∗2)
12         v_x = −ang_speed∗y[i]/norm
13         v_y = ang_speed∗x [i]/norm
14         x[i +1] = x[i]+ dt∗v_x
15         y[i +1] = y[i]+ dt∗v_y
16     return x, y
17
18  if __name__ == ’__main__ ’ :
19      tend=4; nsteps=10000
20      x0 =0.3; y0 =0.5; ang_speed=2
21
22      t=np.linspace(0, tend, nsteps)
23      x,y=evolve(x0, y0, t)
24      fig1 = plt.figure()
25      plt.plot(t, x, ’b−’, t, y, ’r−’)
26      plt.legend([’x’, ’y’], loc=’best’)
27      plt.xlabel(’t’)
28      plt.title(’Explicit Euler Approximation’)
29      plt.grid()
30      plt.show()
31      fig2 = plt.figure()
32      ax = plt.subplot (111, aspect=’equal’)
33      ax.plot(x, y, ’b−’)
34      plt.xlabel(’x’)
35      plt.ylabel(’y’)
36      plt.show()

2.13.2 Object-oriented approach to the Particle Simulator

Now, let’s consider the same problem in an object-oriented way. Suppose, we have to deal with lots of particles and instead of the plot for the particle path we also want to visualize the particles as an animated plot. Then, a natural object-oriented approach of thinking would be to define a Particle class, with attributes such as x, y (coordinates of the position) and ω (angular velocity). We define another class ParticleSimulator where the evolve method will be described that will be responsible for changing the position of the particles obeying the laws of motion (Euler’s method as done in the above code). The particles can be plotted as points using matplotlib.pyplot.plot function as before. We can animate the evolution of particles over time using the matplotlib.animation.FuncAnimation class.

2.13.3 Defining and using classes in Python

A class is defined using the class keyword and indenting the body of statements (attributes and methods) in a block following the declaration. Conventionally the classes names are written in CamelCase. Also, a good habit in programming is to provide docsting describing the purpose of the class after the class statement. The methods in the class are defined using the keyword def as we have used before, but the first argument to each method should be a variable named self. This special name is used to refer to the object itself when it wants to called its own method or refer to attributes.

In the Particle Simulator example, the Particle, the base class can be defined as follows.

1  class Particle:
2     """ Abstract base class representing a particle. """
3      def __init__ (self, x, y, ang_speed):
4         """ Initialize the Particle class with a x, y, ang_speed. """
5          self.x = x
6          self.y = y
7          self.ang_speed = ang_speed
8

Instantiating the object

An instance of a class is created with the syntax object=ClassName (args). An object instantiated from a class should be initialized by setting attributes with appropriate values. To use this simple program, let us save it in a file named simul.py and import it into a new program, or try it in the interactive Python shell with

1  >>> from simul import Particle
2  >>> p=Particle(0.3, 0.5, 2)
3  >>> print(p.x, p.y, p.ang_speed)
4

In the same file simul.py, let us define another class ParticleSimulator and the method evolve which acts on the object particles instantiated from the base class Particle.

1  class ParticleSimulator:
2     """ A class describing particle simulation. """
3
4      def __init__ (self, particles):
5           """ Initialize the ParticleSimulator class with a particle.
      """
6           self.particles = particles
7
8      def evolve(self, dt, nsteps=1):
9            for i in range(nsteps):
10                for p in self.particles:
11 
12                     norm = (p.x∗∗2 + p.y∗∗2)∗∗0.5
13                     v_x = p.ang_speed∗(−p.y)/norm
14                     v_y = p.ang_speed∗p.x/norm
15
16                    d_x = dt ∗ v_x
17                    d_y = dt ∗ v_y
18
19                    p.x += d_x
20                    p.y += d_y
21

Methods and attributes

The class ParticleSimulator has attribute self.particles and the method evolve acts on it. Both attribute and method are accessed using the “dot” notation: object.attr, object.method(). To find the evolution of a particle with x = 0.3, y = 0.5 and ang_speed = 2 after 10 steps with dt = 0.0001, we first instantiate an object p=Particles(0.3,0.5,2) then create a new object simulator=ParticleSimulator(p) and apply the method evolve on it passing the argument dt = 0.1 and nsteps = 10. Try this in an interactive python shell.

1  from simul import Particle, ParticleSimulator
2  particles = [Particle( 0.3, 0.5, +1),
3               Particle( 0.0, −0.5, −1),
4               Particle(−0.1, −0.4, +3)]
5  particles[0].x
6  particles[0].y
7  simulator=ParticleSimulator(particles)
8  simulator.particles[0].x
9  simulator.particles[0].y
10 simulator.evolve(0.0001, 1)
11 simulator.particles[0].x
12 simulator.particles[0].y

We can plot the particle path using matplotlib.pyplot.plot function. Here is a function for plotting the solution.

1  def plot_trajectory(simulator, dt, nsteps):
2      x=[p.x for p in simulator.particles]
3      y=[p.y for p in simulator.particles]
4      x0=[x[0]]; y0=[y[0]]
5      x1=[x[1]]; y1=[y[1]]
6      x2=[x[2]]; y2=[y[2]]
7      for i in range (nsteps):
8          simulator.evolve(dt) # compute new position after one step
9          x=[p.x for p in simulator.particles]
10         y=[p.y for p in simulator.particles]
11         x0.append(x[0]); x1.append(x[1]); x2.append(x[2])
12         y0.append(y[0]); y1.append(y[1]); y2.append(y[2])
13     fig=plt.figure()
14     ax = fig.add_subplot(111, aspect=’equal’)
15     ax.plot(x0, y0, ’r−’, x1, y1, ’g−’, x2, y2, ’b−’)
16     # Axis limits
17     plt.xlim(−1, 1)
18     plt.ylim(−1, 1)
19     plt.show()

Finally, we use the matplotlib.animation.FuncAnimation class to animate the evolution of particles over time.

1  def visualize(simulator):
2      X = [p.x for p in simulator.particles]
3      Y = [p.y for p in simulator.particles]
4
5      fig = plt.figure()
6      ax = plt.subplot(111, aspect=’equal’)
7      line, = ax.plot(X, Y, ’ro’)
8
9      # Axis limits
10     plt.xlim(−1,1)
11     plt.ylim(−1,1)
12
13     # It will be run when the animation starts
14     def init():
15         line.set_data([], [])
16         return line, # The comma is important!
17
18     def animate (i):
19         # We let the particle evolve for 1000 times with dt=0.00001
20         simulator.evolve1(0.00001, 1000)
21         X = [p.x for p in simulator.particles]
22         Y = [p.y for p in simulator.particles]
23
24         line.set_data(X, Y)
25         return line,
26
27     # Call the animate function each 10 ms
28     anim=animation.FuncAnimation (fig,
29                                   animate,
30                                   init_func=init,
31                                   blit=True,
32                                   interval=10)
33     anim.save(’Animation.gif’, writer=’imagemagick’, fps=30)
34     plt.show()

2.14 Introduction to SciPy

SciPy is a library of Python modules for scientific computations. It contains modules for physical/mathematical constants and special functions frequently used in science and engineering, optimization, integration, interpolation and image manipulation. Just like NumPy, underlying algorithms in SciPy are executed as precompiled C program. As such, SciPy functions are fast.