The print function

We're going to be using a bit more of the print function today in our scripts, to help us to see what's happening in our code.

Recall that print() can be used with a value (of any type) inserted directly

print(5)
print("Hello")

It can also print more than one thing at once:

print("Hello", "World")

It can be used to print a variable

x = [1, 5, 7]
print(x)

But what if we want the output to be something like this, where the 5 is a variable?

The value is 5

this is going to be a string with a variable value injected into it. The way this is done is as follows:

a = 5
print("The value is {}".format(a))

{} is a placeholder for the variable that is specified inside format(). If we need to insert more than one variable then we can add them, comma-separated, inside format() as follows

a = 5
b = 6
print("The values are {} and {}".format(a, b))

The placeholders can be labelled so that a variable can be inserted twice in the same string:

a = 5
b = 6
print("The values are {0} and {1}. Here's the first one again...{0}".format(a, b))

Things start to get a little complicated if you want to specify the format of the values you are inserting. Here's an example which prints $\pi$ and $e$ to 2 decimal places:

import numpy as np
print("Some stuff to 2 decimal places... pi = {:.2f} and e = {:.2f}".format(np.pi, np.e))
# or including the labels 
print("Some stuff to 2 decimal places... pi = {0:.2f} and e = {1:.2f}".format(np.pi, np.e))

In {0:.2f}, the 0 is the label (so np.pi is inserted here), the : separates the label from the format specification, .2 is for 2 decimal places and f is for float.

We could spend all day talking about string formatting, and at this stage it will be easier to look up what you want to do when you get to it; here's a guide for reference later, or if you're just really keen.

Logical operators

Recall from Handout 1 that we use logical operators to determine whether a statement is true or false

x == 2
x > 2

Others which can be used to make up such an expression are as follows:

Here are some more examples

a = 3
print(a <= 3)            # True
print(a < 5 or a > 25)   # True
print(a < 5 and a > 25)  # False
print(a % 2 == 0)        # False (is the remainder when divided by 2 zero - in other words is a even)

These will form an important part of the next section on control flow.

Control Flow

This section introduces loops and if statements, part of a family of tools used by programmers and referred to generally as control flow.

for n in range(1,6):
    # do something with n

for n in range(1,6):
# do something with n

for n in range(1, 6):
    print(n**2)

for n in range(1, 6):
    print("The value squared is {}".format(n**2))

for n in [4, 1, 5, 6]:
    print(n)

bears = ["koala", "panda", "sun"]
for x in bears:
    print(x)

# Intitialise a variable
s = 0

# Loop through adding n each time
for n in range(1, 6):
    s += n

# print final value
print(s)

s = 0   # some value we'll add to
n = 0   # this is my counter

while s < 1000:
    n += 1    # increment the counter
    s += n    # add to s  

# Output the results
print("s is equal to {}".format(s))

x = 2
if x >= 2:
    print("that is true")

x = 1
if x >= 2:
    print("that is true")

x = 2
if x < 2:
    print("x is less than 2")
else:
    print("x is not less than 2")

x = 2
if x < 2:
    print("x is less than 2")
elif x == 2:
    print("x is equal to 2")
else:
    print("x is greater than 2")

s = 0
for n in range(100):
    n += 1    # increment the counter
    s += n    # calculate the new sum  
    if s > 1000:
        break

# Output the results
print(s)

x = [2, 5, 8, 5]
print(x[2])

x = [2, 5, 8, 5]
x[2] = 4
print(x)

t = []
for n in range(0, 10):
    t[n] = n**2

t = []
for n in range(0,10):
    t.append(n**2)

import numpy as np
t = np.zeros(10)
for n in range(0,10):
    t[n] = n**2

print(t)

n = np.arange(0, 10)
t = n**2

For loops plot example

import matplotlib.pyplot as plt
import numpy as np

# Array of x values
x = np.linspace(0,8,100)

# Plot f(x)=xe^(ax) for a from 1 to 5
for a in range(1,6):
    f = x*np.exp(-a*x)
    plt.plot(x,f)


plt.show()

Which produces this plot:

Download the full source code for this plot

Functions

Python makes it possible to write your own functions, which take some input and return a value or values, in just the same way as Python's built-in functions. This helps to keep your Python code as modular as possible.

The syntax for creating a function is as follows:

def my_func():
    print("My function prints this")

Note a similar syntax as for control flow: the function begins with the keyword def and then the function name "my_func". This is followed by input arguments inside brackets - for this function there are none, and finally a colon. The contents of the function are then indented.

We can call the function with

my_func()

either from the same file or the Console.

Now let's add an input argument to our function and a more descriptive name:

def zoo_visit(animal):
    print("I went to the zoo and saw a {}".format(animal))

zoo_visit("koala")
zoo_visit("panda")
zoo_visit("sun bear")

We can also set a default value by using argument = ... in the parentheses. This default is used if the input argument is not set.

def zoo_visit(animal="bear"):
    print("I went to the zoo and saw a {}".format(animal))

zoo_visit("koala")
zoo_visit("panda")
zoo_visit("sun bear")
zoo_visit()

Any data type can be sent to a function. Here's a list, with a for loop to print each value - pay careful attention to the indenting:

def print_my_list(list):
    for x in list:
        print(x)

print_my_list([1, 5, 2, 6])

And a simple one with a number

def square_a_number(x):
    print(x**2)

square_a_number(2)

All of the above examples print something. The functions we have been using however return a value which can be assigned to a variable. For example at the moment this does not do what we might like

x = square_a_number(2)

To return a value (or values) from a function we need to use a return statement. This is done as follows:

def square_a_number(x):
    return x**2

x = square_a_number(5)
print(x)

Note that the x in the function argument and the x outside the function are completely unrelated. This is known as the scope of a variable.

Now let's extend this by accepting two input arguments:

def show_me_the_bigger(a, b):
    return max([a, b])

x = show_me_the_bigger(4, 5)
print(x)

def calculate_something(a):
    return ((a + 7)/2.5)**6

a_list = [1, 6, 3.4, 27, 5.12]

x = [calculate_something(a) for a in a_list]
print(x)

Adding help to your function

A comment contained within three quotes """ at the start of our custom function is used to display help. It is known as a docstring (documentation string)

import numpy as np
import matplotlib.pyplot as plt

def sin_plus_cos(x):
    """
    Takes in a value x and 
    returns cos(x)+sin(x)
    """
    return np.cos(x)+np.sin(x)

Test your help with

help(sin_plus_cos)

def  inverse(x):
    return 1/x    

inverse(2)

inverse(0)

Traceback (most recent call last):

  File "<ipython-input-72-1bc64584a5ff>", line 1, in <module>
    runfile('/Users/georgestagg/python/tests/untitled0.py', wdir='/Users/georgestagg/python/tests')

  File "/Users/georgestagg/anaconda3/lib/python3.7/site-packages/spyder_kernels/customize/spydercustomize.py", line 786, in runfile
    execfile(filename, namespace)

  File "/Users/georgestagg/anaconda3/lib/python3.7/site-packages/spyder_kernels/customize/spydercustomize.py", line 110, in execfile
    exec(compile(f.read(), filename, 'exec'), namespace)

  File "/Users/georgestagg/python/tests/untitled0.py", line 6, in <module>
    fraction(1,0)

  File "/Users/georgestagg/python/tests/untitled0.py", line 3, in fraction
    print("as a decimal: {}".format(numerator/denominator))

ZeroDivisionError: division by zero

def inverse(x):
  try:
    return 1/x
  except:
    print("Something went wrong")

inverse(0)
inverse("2")

def inverse(x):
  try:
    return 1/x
  except ZeroDivisionError:
    print("Error: please enter a non-zero value")
  except TypeError:
    print("Error: input argument should be a float or integer")

inverse(0)
inverse("2")
inverse(2)

def inverse(x):
    if x == 0:
        raise ZeroDivisionError("Please enter a non-zero denominator")
    return 1/x

inverse(0)

Writing algorithms

When does a piece of code become an algorithm? For us we're pretty much there already at that point.

An algorithm is a set of step-by-step instructions, in our case carried out by our Python code, to solve a problem.

As an analogy, in the problem of baking a cake, a recipe book lays out the step-by-step instructions to achieve the tasty goal. The detail is such that, if the instructions are followed precisely, a second cake produced using the same recipe will be identical.

Our goal in the next two worked examples will be to identify what the steps are to solving problems presented to us, and then how to code them in Python.

Here are two worked examples:

Worked example: Fibonacci sequence

In the Fibonacci sequence, each number is the sum of the two preceding ones:

where typically the seed values $F_{0}=0,F_{1}=1$ are used.

Our task is to write an algorithm that puts (say the first 15) values of the sequence into a vector.

import numpy as np

# Set up an empty vector F for Fibonacci values
F = np.zeros(15)

# Set the seed values
F[0] = 0
F[1] = 1

# loop through from n = 2 upwards to get other values
for n in range(2, 15):
   F[n] = F[n-1]+F[n-2]

# Take a look at the output
print(F)

Worked example 2: Euler's Method

Spoiler alert: Later in your degree you will be confronted with all sorts of problems involving differential equations, as many physical behaviours can be easily described with them. A differential equation is simply an equation which relates some function(s) to its derivatives, for example

This sort of differential equation (involving the derivative of a function which depends on just one variable) is known as an ordinary differential equation (ODE) and is usually associated with some sort of initial condition, for example the value at $y(0)$, in which case it is known as an initial value problem.

To name but a few differential equations: Newton's second law, Maxwell's equations of electromagnetism, radioactive decay, Newton's law of cooling... some of these you will have already met in some form... others you will meet: the heat equation (thermodynamics), the Schrödinger equation (quantum mechanics), Navier-Stokes (fluid dynamics), and many more...

It is interesting for us to look at solving differential equations for several reasons. One, that it is clearly going to be useful to your later study. But secondly, the oldest algorithm to compute a numerical solution to straight forward problems of this sort is Euler’s method, which involves little more than a for loop.

Suppose that we know the value of $y(t_0)$ = $y_0$ and we want to find $y(t_1)$ = $y_1$ for some equation in the form

If $h = t_1-t_0$ is very small, then we can reasonably expect that $y_1$ will be close to the tangent line to $y(t)$ at $t_0$. The gradient of that tangent line is $m = f(t_0,y_0)$, the right hand side of the ODE at $t_0$. So, from elementary geometry,

And since $t_1 = h + t_0$, we obtain,

Now we can use the same method to find $y_2$ from $y_1$, and so on. So the Euler Method approximates the solution using

Here's a visual of how it works:

The accuracy will clearly depend on choosing $h$ carefully. As we can see, as the stepsize $h$ is reduced, the solution approaches the exact (red curve below)

Unfortunately, small $h$ means more iterations. Even more unfortunately, Euler’s method is often not stable, meaning that the error in the approximation can quickly accumulate to such a size that the numerical solution diverges wildly from the true solution.

The main value in Euler’s method is that it illustrates an important principle. Other, more reliable, methods use the same basic idea to find numerical solutions to differential equations, but use more information about $y(t)$ to move from one point to the next.

Coding the algorithm

The general problem is:

$$\frac{\mathrm{d}y}{\mathrm{d}t}=f(t,y).$$

and we would like to code up

$$y_{n} = y_{n-1} + hf(t_{n-1},y_{n-1}).$$

for some function $f$ which might depend on $t$ and $y$, some stepsize in $t$ given by $h$, and where $y_0$ is known. It's very similar in many ways to the Fibonacci example. In my example above, the differential equation is

$$\frac{\mathrm{d}y}{\mathrm{d}t}=-\frac{y}{2},\quad y(0)=5.$$

So in fact my function will only depend on $y$. Here's my plan:

• Set up an function to do the job of $f(y)=-y/2$
• Choose a step-size for $h$: let's pick $h=0.5$.
• Choose how many values of $y$ I'm going to compute - let's say 20 - and initialise a vector for y using np.zeros.
• Set the initial value y[0] = 0.
• Write a for loop to calculate the remaining values of y[n], based on the value of y[n-1].

Here we go:

import numpy as np
import matplotlib.pyplot as plt

# A function for the RHS of the ODE
def f(y):
    return -y/2

# Time step
h = 0.5

# Initialise y and set value when t = 0
y = np.zeros(20)
y[0] = 5

# For loop to calculate remaining values
for n in range(1, 20):
    y[n] = y[n-1]+h*f(y[n-1])

If we want to plot $y$ versus $t$ we will need a vector for $t$. Given the timestep $h=0.5$, we know that this is going to be

t = np.arange(0, 10, 0.5)

So now we could plot

plt.plot(t,y)

Download the full source code for this plot

Exercise 3.5 (Harder!)

• Use the code above as a template (i.e. copy and paste it!) and change the problem to solve

$$\frac{\mathrm{d}y}{\mathrm{d}t}=0.5y(1-y),\quad y(0)=0.5$$

You should need to change only two lines in the code (the function and the initial value).

• Write a for loop over the whole bit of code above to change the value of $y(0)$ to different values between 0 and 1. It might look something like this:

for y0 in np.linspace(0, 1, 20):
    # Euler's method code as above in here
    # except set the initial condition with y[0] = y0

If you plot all of the solutions then you should get something like this:

The differential equation is a rather famous one, and is used for population modelling. The $y$ here represents the proportion of a population area that is filled. And the $0.5$ in the RHS of the ODE is a parameter which determines the rate of growth/decline (e.g. try changing it to -0.5 to see decline).

Next time

That's it for this week! We now have all the skills to import data, run algorithms on it, and produce beautiful plots of output. We'll expand on the data analysis side of things next week, as we look at some curve fitting techniques.