P.D. Nation and J.R. Johansson

For more information about QuTiP see http://qutip.org

Here we import the required functions for later usage.

In [1]:

```
import matplotlib.pyplot as plt
import numpy as np
from scipy.special import factorial
```

**iPython** as a basic calculator. Addition, subtraction, and multiplication, all work in the same way as you would write the equations down on paper

In [2]:

```
10 + 5
```

Out[2]:

15

In [3]:

```
10 - 157
```

Out[3]:

-147

In [4]:

```
4 / 3
```

Out[4]:

1.3333333333333333

In [5]:

```
(50 - 4) * 10 / 5
```

Out[5]:

92.0

However, things like raising a number to a power, $4^{4}$, are written differently.

In [6]:

```
4**4
```

Out[6]:

256

In [7]:

```
4**4.0
```

Out[7]:

256.0

All information stored in a computer must be represented in a binary format consisting of zeros and ones (e.g. $461\rightarrow 111001101$). Each zero or one is called a **bit**, and given $N$ bits, one can store all of the integers in the range $[0,2^{N-1}]$, where the $-1$ is due to the fact that the first bit is reserved for defining if a number is positive or negative

However, given a fixed number of bits, it is impossible to store an arbitrary number exactly. Therefore, if one is given a random number, unless the number is exactly divisible by a factor of two, the conversion between the random number and the binary bit representation ultimately leads to a loss of precision, and what is known as **roundoff error**.

When dealing with numbers inside a computer there are two distinct types of numbers to consider:

**Integers**- (1,2,4,-586,..) Are what are called**fixed-point numbers**, where the term fixed-point means that there is a fixed number of decimal places in the number (zero for integers). These numbers can be stored exactly in a computer.

**Doubles/Floats**- (3.141,0.21,-0.1,..) These are**floating-point numbers**that are the binary equivalent to scientific notation $c=2.99792458\times 10^{8}$. Doubles (also called double-precision numbers) are floating point numbers that are written using 64-bits and, in general, are only accurate to the 15th or 16th decimal place. Floats (or single-precision numbers) use 32-bits, and are good to 6-7 decimal places.**Serious scientific calculations always require a combination of integers and double (64-bit) numbers**.

In [8]:

```
7 + 0.000000000000001
```

Out[8]:

7.000000000000001

In [9]:

```
7 + 0.0000000000000001
```

Out[9]:

7.0

In [10]:

```
0.1 + 0.2
```

Out[10]:

0.30000000000000004

This last example clearly highlights the fact that the computer does not store decimal (floating-point) numbers exactly. The loss of precision in floating-point numbers can be characterized by the **machine precision**, $\epsilon_{\rm m}$, that is defined to be the smallest positive number such that

where the subscript on $1_{\rm c}$ is meant to remind you that this is a computer number. Therefore, for any arbitrary number $N$ is related to its floating-point equivalent $N_{\rm c}$ by

$$N_{\rm c}=N\pm \epsilon, \ \ \forall~|\epsilon|< \epsilon_{\rm m}.$$**Take Home Message** - All double-precision decimal numbers that are not factors of two will have error in the 15th decimal place. This can lead to errors in your numerical solutions if you are not careful.

Python itself has limited support for mathematics outside of simple arithmetic. Therefore, we will use the functions in the NumPy module to do more impressive, and faster, calculations. We have imported NumPy already at the top of this notebook and can use it now by referring to `np`

.

We can now do more impressive calculations:

In [11]:

```
np.exp(2.34)
```

Out[11]:

10.381236562731843

In [12]:

```
np.sqrt(5)
```

Out[12]:

2.23606797749979

In [13]:

```
np.sinc(0.5)
```

Out[13]:

0.6366197723675814

In [14]:

```
radius = 5
area = np.pi * radius**2
area
```

Out[14]:

78.53981633974483

We see that our variables name is defined on the left of the `=`

sign and the value its given is defined on the right. Here we have also used the `pi`

variable that has been predefined by NumPy. Variables can then be used in other expressions.

If a predefined variable is again used on the left side of `=`

then its original value is replaced.

In [15]:

```
x = 10
x = (x**2 + 25) / 10
x
```

Out[15]:

12.5

`=`

sign in a computer program is **not** equivalent to the mathematical equality.

What happens if you try to use a variable without first defining it?

Python would give us an error that the variable is not defined. In addition, there are several words that are reserved by the Python language and cannot be used as variables:

```
and, as, assert, break, class, continue, def, del, elif, else, except,
exec, finally, for, from, global, if, import, in, is, lambda, not, or,
pass, print, raise, return, try, while, with, yield
```

Other than the above reserved words, your variables can be anything that starts with a letter or the underscore character "$\_$" followed by any combination of alphanumeric characters and "$\_$". Note that using upper or lower case letters will give you two different variables.

In [16]:

```
_freq = 8
Oscillator_Energy = 10
_freq * Oscillator_Energy
```

Out[16]:

80

In [17]:

```
speed_of_light = 2.9979 * 10**8
spring_constant = np.sqrt(2 / 5)
```

**strings**. We have already seen one string already in this class:

In [18]:

```
"Hello Class"
```

Out[18]:

'Hello Class'

We can also use single quotes, e.g. `'Hello Class'`

.

If we want to use the quote symbol in the string itself then we need to mix the two types

In [19]:

```
"How was Hwajung's birthday party?"
```

Out[19]:

"How was Hwajung's birthday party?"

In [20]:

```
a = "I like " # There is a blank space at the end of this string.
b = "chicken and HOF"
a + b
```

Out[20]:

'I like chicken and HOF'

`print`

function to accomplish this

In [21]:

```
temp = 23
text = "The temperature right now is"
print(text, temp)
```

The temperature right now is 23

`print`

function automatically takes any number of string, integer, double, or other variables, converts them into strings, and then prints them for the user.

** list** datatype variable.

In [22]:

```
shopping_list = ["eggs", "bread", "milk", "bananas"]
```

**index** that corresponds to the variable inside of square brackets.

In [23]:

```
shopping_list[2]
```

Out[23]:

'milk'

In [24]:

```
shopping_list[0]
```

Out[24]:

'eggs'

In [25]:

```
shopping_list[-1]
```

Out[25]:

'bananas'

In [26]:

```
shopping_list[-2]
```

Out[26]:

'milk'

`len`

function that returns an integer giving the length of the list.

In [27]:

```
len(shopping_list)
```

Out[27]:

4

`append`

and `remove`

, respectively.

In [28]:

```
shopping_list.append("apples")
shopping_list
```

Out[28]:

['eggs', 'bread', 'milk', 'bananas', 'apples']

In [29]:

```
shopping_list.remove("bread")
shopping_list
```

Out[29]:

['eggs', 'milk', 'bananas', 'apples']

In [30]:

```
various_things = [1, "hello", -1.234, [-1, -2, -3]]
various_things
```

Out[30]:

[1, 'hello', -1.234, [-1, -2, -3]]

All of these elements can be accessed in the usual way

In [31]:

```
various_things[0]
```

Out[31]:

1

In [32]:

```
various_things[-1]
```

Out[32]:

[-1, -2, -3]

In [33]:

```
various_things[3][1]
```

Out[33]:

-2

**iteration** and is accomplished in Python using the `for`

command:

In [34]:

```
items = [
"four calling birds",
"three french hens",
"two turtle doves",
"a partridge in a pear tree",
]
for thing in items:
print(thing)
```

four calling birds three french hens two turtle doves a partridge in a pear tree

`print`

function. We are free to call this variable anything we want.

In [35]:

```
for variable in items:
print(variable)
```

four calling birds three french hens two turtle doves a partridge in a pear tree

**block**. If we did not intent the print function then Python would yell at us.

In [36]:

```
for variable in items:
print("My true love gave to me", variable)
```

**slicing** to conveniently access the elements. Slicing can be used on any **sequence** such as lists, strings, and as we will see shortly, arrays. Consider our `shopping_list`

list:

In [37]:

```
shopping_list = ["eggs", "bread", "milk", "bananas", "apples"]
```

To get the first element we used a single index

In [38]:

```
shopping_list[0]
```

Out[38]:

'eggs'

But if we want to get the first three elements in the list we can use:

In [39]:

```
shopping_list[0:3]
```

Out[39]:

['eggs', 'bread', 'milk']

We could also grab the last two elements using:

In [40]:

```
shopping_list[-2:]
```

Out[40]:

['bananas', 'apples']

In [41]:

```
shopping_list[0::2]
```

Out[41]:

['eggs', 'milk', 'apples']

**conditional statements**. The basic operations in boolean logic are "equal" (`==`

), "not equal" (`!=`

), "greater than" (`>`

), "greater than or equal" (`>=`

), "less than" (`<`

), and "less than or equal" (`<=`

). All of these conditionals operate on two variables and return a simple boolean `True`

or `False`

answer. For example

In [42]:

```
a = 5
b = 8
a > b
```

Out[42]:

False

In [43]:

```
c = 0
c <= 0, c >= 0
```

Out[43]:

(True, True)

In [44]:

```
a = 5
b = 6
a == b, a != b
```

Out[44]:

(False, True)

It is important to point out that in Python `1`

and `0`

are the same as `True`

and `False`

, respectively.

In [45]:

```
t = True
f = False
t == 1, f == 0
```

Out[45]:

(True, True)

We can also combine multiple conditional statements

In [46]:

```
a = -1
b = 4
c = 10
d = 11
a < b < c != d
```

Out[46]:

True

These operations can also be used on lists and strings:

In [47]:

```
[4, 5, 6] >= [4, 5, 7]
```

Out[47]:

False

In [48]:

```
[4, 5, 6] <= [4, 5, 7]
```

Out[48]:

True

In [49]:

```
"today" == "Today"
```

Out[49]:

False

`if/else`

and `while`

statements.

In [50]:

```
today = "friday"
if today == "friday":
print("We have class today :(") # this is a code block
else:
print("No class today :)") # this is also a code block
```

We have class today :(

`if`

statement is run only if the conditional `today=='friday'`

returns `True`

. If the conditional is `False`

then the code block inside the `else`

statement is run. We can also check multiple conditions by using the `elif`

statement after `if`

:

In [51]:

```
today = "thursday"
if today == "friday":
print("We have class today :(")
elif today == "thursday":
print("Our assignment is due today :(")
else:
print("No class today :)")
```

Our assignment is due today :(

** while loop** that executes a block of code repeatedly until the conditional statement at the start of the loop is

`False`

.In [52]:

```
n = 0
while n <= 10: # evaluate code block until n>10
print("The current value of n is:", n)
n = n + 1 # increase the value of n by 1
```

`while`

loop you must make sure the conditional is not `True`

forever. Otherwise your program will be in an **infinite loop** that never ends.

Let us determine whether a given number between [1,10] is an even or odd number.

In [53]:

```
for n in [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]:
if np.remainder(n, 2) == 0:
print(n, "is even")
else:
print(n, "is odd")
```

1 is odd 2 is even 3 is odd 4 is even 5 is odd 6 is even 7 is odd 8 is even 9 is odd 10 is even

`range`

that makes creating sequences of integers very easy. For instance, the above example becomes

In [54]:

```
for n in range(1, 11):
if np.remainder(n, 2) == 0:
print(n, "is even")
else:
print(n, "is odd")
```

1 is odd 2 is even 3 is odd 4 is even 5 is odd 6 is even 7 is odd 8 is even 9 is odd 10 is even

*never* part of the generated sequence when using `range`

. If we wanted the `range`

function to start at zero instead of one we could simply write `range(11)`

. We can also make sequences that go in arbitrary steps:

In [55]:

```
for n in range(0, 11, 2):
print(n)
```

0 2 4 6 8 10

`range`

function does not return a list of integers but is something called a **generator**. In general, the `range`

function should only be used in combination with the `for`

command.

In [56]:

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

[0, 1, 1, 2, 3, 5, 8, 13, 21, 34]

We can also write this using a `while`

loop if we wanted to.

In [57]:

```
n = 2
fib = [0, 1]
while n < 10:
fib.append(fib[n - 1] + fib[n - 2])
n = n + 1
print(fib)
```

[0, 1, 1, 2, 3, 5, 8, 13, 21, 34]

**scripts** that contain a collection of constants, variables, data structures, functions, comments, etc., that perform various complicated tasks.

**.py** extension. Python scripts are also called Python **programs**. If we open up any editor, then we are given a blank window that we can enter our Python commands in.

Before we begin to write our scripts, lets first discuss the best format for writing your scripts.

In [58]:

```
# This is an example script for the P461 class
# Here we will calculate the series expansion
# for sin(x) up to an arbitrary order N.
#
# Paul Nation, 02/03/2014
N = 5 # The order of the series expansion
x = np.pi / 4.0 # The point at which we want to evaluate sine
ans = 0.0
for k in range(N + 1):
ans = ans + (-1) ** k * x ** (1 + 2 * k) / factorial(1 + 2 * k)
print("Series approximation:", ans)
print("Error:", np.sin(x) - ans)
```

Series approximation: 0.7071067811796194 Error: 6.928124740568364e-12

**comments** that describe what the script does and when it was created. In python all comments start with the ** #** symbol. Everything after this symbol is ignored by the computer. Second, we have the section of the scripts that load the necessary functions that we need from other packages. Third is a section where we define all of the constants that are going to be used in the script. You should also add comments here that tell us what the constants are. Finally, your main body of code goes after these sections.

**functions**. Functions are blocks of code that accomplish a specific task. Functions usually take "input arguments", perform operations on these inputs, and then "return" one or more results. Functions can be used over and over again, and can also be "called" from the inside of other functions. Let us rewrite our script for $sin(x)$ using a function and then describe each part.

In [59]:

```
N = 5 # The order of the series expansion
x = np.pi / 4.0 # The point at which we want to evaluate sine
def sine_series(x, N):
ans = 0.0
for k in range(N + 1):
ans = ans + (-1) ** k * x ** (1 + 2 * k) / factorial(1 + 2 * k)
return ans
result = sine_series(x, N)
print("Series approximation:", result)
print("Error:", np.sin(x) - result)
```

Series approximation: 0.7071067811796194 Error: 6.928124740568364e-12

`def`

which is short "define", then the name of the function followed by the input arguments in parentheses. After the block of code called by the function, the `return`

keyword specifies what variable(s) and/or data structure(s) are given as the output. So a general functions call is

In [60]:

```
def function_name(arg1, arg2):
"Block of code to run"
"..."
return result
```

Again, everything after the colon (:) that is inside the function must be indented. The beauty of using functions is that we can use the same code over and over, just by changing the constants near the top of our Python script.

Variables that are defined inside of a function are called **local variables** and only defined for the block of code inside of the function. In our previous example, `k`

was a local variable. The input arguments and return arguments are *not* local variables. Once a function is done running, the local variables are erased from memory. Therefore, if you want get something out of a function, your must return the value when your done.

In [61]:

```
N = 100 # Number of points to generate
def random_coordinates(N):
x_coords = []
y_coords = []
for n in range(N):
xnew, ynew = np.random.random(2)
x_coords.append(xnew)
y_coords.append(ynew)
return x_coords, y_coords
xc, yc = random_coordinates(N)
plt.plot(xc, yc, "ro", markersize=8)
plt.show()
```

In [62]:

```
N = 20 # Number of points to generate
def random_coordinates(N):
x_coords = []
y_coords = []
for n in range(N):
xnew, ynew = np.random.random(2)
x_coords.append(xnew)
y_coords.append(ynew)
return x_coords, y_coords
def dist2d(x1, y1, x2, y2):
return np.sqrt((x1 - x2) ** 2 + (y1 - y2) ** 2)
def max_dist(xc, yc):
max_dist = 0.0
num_points = len(xc)
for ii in range(num_points):
for jj in range(num_points):
dist = dist2d(xc[ii], yc[ii], xc[jj], yc[jj])
if dist > max_dist:
max_dist = dist
xvals = [xc[ii], xc[jj]]
yvals = [yc[ii], yc[jj]]
return max_dist, xvals, yvals
xc, yc = random_coordinates(N)
max_dist, pnt1, pnt2 = max_dist(xc, yc)
plt.plot(xc, yc, "ro", markersize=8)
plt.plot(pnt1, pnt2, "b-", lw=2)
plt.show()
```

`max_dist`

function:

In [63]:

```
def max_dist(xc, yc):
"""
Finds the maximum distance between any two points
in a collection of 2D points. The points corresponding
to this distance are also returned.
Parameters
----------
xc : list
List of x-coordinates
yc : list
List of y-coordinates
Returns
-------
max_dist : float
Maximum distance
xvals : list
x-coodinates of two points
yvals : list
y-coordinates of two points
"""
max_dist = 0.0 # initialize max_dist
num_points = len(xc) # number of points in collection
for ii in range(num_points):
for jj in range(num_points):
dist = dist2d(xc[ii], yc[ii], xc[jj], yc[jj])
if dist > max_dist:
max_dist = dist
xvals = [xc[ii], xc[jj]]
yvals = [yc[ii], yc[jj]]
return max_dist, xvals, yvals
```

`"""..."""`

is called a **docstring** and it gives a tells someone who is not familiar with a partiular functions a detailed explaination as to what the function does, what parameters it takes as inputs, and what values it returns. It is also good practice to put some comments next to your local variables so the user knows what each of these is for. Although it seems like a lot of work at first, writing docstrings will make you a much better programmer in the future.