# Ch.3: Functions and branching

Hans Petter Langtangen [1, 2]

[1] Simula Research Laboratory
[2] University of Oslo, Dept. of Informatics

## We have used many Python functions

Mathematical functions:

from math import *
y = sin(x)*log(x)


Other functions:

n = len(somelist)
integers = range(5, n, 2)


Functions used with the dot syntax (called methods):

C = [5, 10, 40, 45]
i = C.index(10)        # result: i=1
C.append(50)
C.insert(2, 20)


What is a function? So far we have seen that we put some objects in and sometimes get an object (result) out of functions. Now it is time to write our own functions!

## Functions are one of the most import tools in programming

• Function = a collection of statements we can execute wherever and whenever we want
• Function can take input objects (arguments) and produce output objects (returned results)
• Functions help to organize programs, make them more understandable, shorter, reusable, and easier to extend

## Python function for implementing a mathematical function

The mathematical function $$F(C)={9\over5}C+32$$

can be implemented in Python as follows:

def F(C):
return (9.0/5)*C + 32


Note:

• Functions start with def, then the name of the function, then a list of arguments (here C) - the function header
• Inside the function: statements - the function body
• Wherever we want, inside the function, we can "stop the function" and return as many values/variables we want

## Functions must be called

A function does not do anything before it is called.

Note:

The call F(C) produces (returns) a float object, which means that F(C) is replaced by this float object. We can therefore make the call F(C) everywhere a float can be used.

## Functions can have as many arguments as you like

Make a Python function of the mathematical function $$y(t) = v_0t- \frac{1}{2}gt^2$$

## Function arguments become local variables

Local vs global variables.

When calling yfunc(t, 3), all these statements are in fact executed:

t = 0.6  # arguments get values as in standard assignments
v0 = 3
g = 9.81
return v0*t - 0.5*g*t**2


Inside yfunc, t, v0, and g are local variables, not visible outside yfunc and desroyed after return.

Outside yfunc (in the main program), t, v0, and y are global variables, visible everywhere.

## Functions may access global variables

The yfunc(t,v0) function took two arguments. Could implement $$y(t)$$ as a function of $$t$$ only:

>>> def yfunc(t):
...     g = 9.81
...     return v0*t - 0.5*g*t**2
...
>>> t = 0.6
>>> yfunc(t)
...
NameError: global name 'v0' is not defined


Problem: v0 must be defined in the calling program program before we call yfunc!

>>> v0 = 5
>>> yfunc(0.6)
1.2342


Note: v0 and t (in the main program) are global variables, while the t in yfunc is a local variable.

## Local variables hide global variables of the same name

Test this:

Question.

What gets printed?

## Global variables can be changed if declared global

What gets printed?

1. v0: 2
4.0608
2. v0: 9


What happens if we comment out global v0?

1. v0: 2
4.0608
2. v0: 2


v0 in yfunc becomes a local variable (i.e., we have two v0)

## Functions can return multiple values

Say we want to compute $$y(t)$$ and $$y'(t)=v_0-gt$$:

def yfunc(t, v0):
g = 9.81
y = v0*t - 0.5*g*t**2
dydt = v0 - g*t
return y, dydt

# call:
position, velocity = yfunc(0.6, 3)


Separate the objects to be returned by comma, assign to variables separated by comma. Actually, a tuple is returned:

>>> def f(x):
...     return x, x**2, x**4
...
>>> s = f(2)
>>> s
(2, 4, 16)
>>> type(s)
<type 'tuple'>
>>> x, x2, x4 = f(2)   # same syntax as x, y = (obj1, obj2)


## Example: Compute a function defined as a sum

The function $$L(x;n) = \sum_{i=1}^n {1\over i}\left( {x\over 1+x}\right)^{i}$$ is an approximation to $$\ln (1+x)$$ for a finite $$n$$ and $$x\geq 1$$.

Corresponding Python function for $$L(x;n)$$:

def L(x, n):
x = float(x)  # ensure float division below
s = 0
for i in range(1, n+1):
s += (1.0/i)*(x/(1+x))**i
return s

x = 5
from math import log as ln
print L(x, 10), L(x, 100), ln(1+x)


## Returning errors as well from the L(x, n) function

We can return more: 1) the first neglected term in the sum and 2) the error ($$\ln (1+x) - L(x;n)$$):

def L2(x, n):
x = float(x)
s = 0
for i in range(1, n+1):
s += (1.0/i)*(x/(1+x))**i
value_of_sum = s
first_neglected_term = (1.0/(n+1))*(x/(1+x))**(n+1)
from math import log
exact_error = log(1+x) - value_of_sum
return value_of_sum, first_neglected_term, exact_error

# typical call:
x = 1.2; n = 100
value, approximate_error, exact_error = L2(x, n)


## Functions do not need to return objects

def somefunc(obj):
print obj

return_value = somefunc(3.4)


Here, return_value becomes None because if we do not explicitly return something, Python will insert return None.

## Example on a function without return value

Make a table of $$L(x;n)$$ vs. $$\ln (1+x)$$:

def table(x):
print '\nx=%g, ln(1+x)=%g' % (x, log(1+x))
for n in [1, 2, 10, 100, 500]:
value, next, error = L2(x, n)
print 'n=%-4d %-10g  (next term: %8.2e  '\
'error: %8.2e)' % (n, value, next, error)


No need to return anything here - the purpose is to print.

x=10, ln(1+x)=2.3979
n=1    0.909091    (next term: 4.13e-01  error: 1.49e+00)
n=2    1.32231     (next term: 2.50e-01  error: 1.08e+00)
n=10   2.17907     (next term: 3.19e-02  error: 2.19e-01)
n=100  2.39789     (next term: 6.53e-07  error: 6.59e-06)
n=500  2.3979      (next term: 3.65e-24  error: 6.22e-15)


## Keyword arguments are useful to simplify function calls and help document the arguments

Functions can have arguments of the form name=value, called keyword arguments:

def somefunc(arg1, arg2, kwarg1=True, kwarg2=0):
print arg1, arg2, kwarg1, kwarg2


## Examples on calling functions with keyword arguments

>>> def somefunc(arg1, arg2, kwarg1=True, kwarg2=0):
>>>     print arg1, arg2, kwarg1, kwarg2

>>> somefunc('Hello', [1,2])   # drop kwarg1 and kwarg2
Hello [1, 2] True 0            # default values are used

>>> somefunc('Hello', [1,2], kwarg1='Hi')
Hello [1, 2] Hi 0              # kwarg2 has default value

>>> somefunc('Hello', [1,2], kwarg2='Hi')
Hello [1, 2] True Hi           # kwarg1 has default value

>>> somefunc('Hello', [1,2], kwarg2='Hi', kwarg1=6)
Hello [1, 2] 6 Hi              # specify all args


If we use name=value for all arguments in the call, their sequence can in fact be arbitrary:

>>> somefunc(kwarg2='Hello', arg1='Hi', kwarg1=6, arg2=[2])
Hi [2] 6 Hello


## How to implement a mathematical function of one variable, but with additional parameteres?

Consider a function of $$t$$, with parameters $$A$$, $$a$$, and $$\omega$$: $$f(t; A,a, \omega) = Ae^{-at}\sin (\omega t)$$

Possible implementation.

Python function with $$t$$ as positional argument, and $$A$$, $$a$$, and $$\omega$$ as keyword arguments:

from math import pi, exp, sin

def f(t, A=1, a=1, omega=2*pi):
return A*exp(-a*t)*sin(omega*t)

v1 = f(0.2)
v2 = f(0.2, omega=1)
v2 = f(0.2, 1, 3)  # same as f(0.2, A=1, a=3)
v3 = f(0.2, omega=1, A=2.5)
v4 = f(A=5, a=0.1, omega=1, t=1.3)
v5 = f(t=0.2, A=9)
v6 = f(t=0.2, 9)   # illegal: keyword arg before positional


## Doc strings are used to document the usage of a function

Important Python convention:

Document the purpose of a function, its arguments, and its return values in a doc string - a (triple-quoted) string written right after the function header.

def C2F(C):
"""Convert Celsius degrees (C) to Fahrenheit."""
return (9.0/5)*C + 32

def line(x0, y0, x1, y1):
"""
Compute the coefficients a and b in the mathematical
expression for a straight line y = a*x + b that goes
through two points (x0, y0) and (x1, y1).

x0, y0: a point on the line (floats).
x1, y1: another point on the line (floats).
return: a, b (floats) for the line (y=a*x+b).
"""
a = (y1 - y0)/(x1 - x0)
b = y0 - a*x0
return a, b


## Python convention: input is function arguments, output is returned

• A function can have three types of input and output data:
• input data specified through positional/keyword arguments
• input/output data given as positional/keyword arguments that will be modified and returned
• output data created inside the function
• All output data are returned, all input data are arguments

def somefunc(i1, i2, i3, io4, io5, i6=value1, io7=value2):
# modify io4, io5, io7; compute o1, o2, o3
return o1, o2, o3, io4, io5, io7


The function arguments are

• pure input: i1, i2, i3, i6
• input and output: io4, io5, io7

## The main program is the set of statements outside functions

from math import *          # in main

def f(x):                   # in main
e = exp(-0.1*x)
s = sin(6*pi*x)
return e*s

x = 2                       # in main
y = f(x)                    # in main
print 'f(%g)=%g' % (x, y)   # in main


The execution starts with the first statement in the main program and proceeds line by line, top to bottom.

def statements define a function, but the statements inside the function are not executed before the function is called.

## Python functions as arguments to Python functions

• Programs doing calculus frequently need to have functions as arguments in other functions, e.g.,
• numerical integration: $$\int_a^b f(x)dx$$
• numerical differentiation: $$f'(x)$$
• numerical root finding: $$f(x)=0$$
• All three cases need $$f$$ as a Python function f(x)

Example: numerical computation of $$f''(x)$$.

$$f''(x) \approx {f(x-h) - 2f(x) + f(x+h)\over h^2}$$

def diff2(f, x, h=1E-6):
r = (f(x-h) - 2*f(x) + f(x+h))/float(h*h)
return r


No difficulty with f being a function (more complicated in Matlab, C, C++, Fortran, Java, ...).

## Application of the diff2 function (read the output!)

Code:

def g(t):
return t**(-6)

# make table of g''(t) for 13 h values:
for k in range(1,14):
h = 10**(-k)
print 'h=%.0e: %.5f' % (h, diff2(g, 1, h))


Output ($$g''(1)=42$$):

h=1e-01: 44.61504
h=1e-02: 42.02521
h=1e-03: 42.00025
h=1e-04: 42.00000
h=1e-05: 41.99999
h=1e-06: 42.00074
h=1e-07: 41.94423
h=1e-08: 47.73959
h=1e-09: -666.13381
h=1e-10: 0.00000
h=1e-11: 0.00000
h=1e-12: -666133814.77509
h=1e-13: 66613381477.50939


## Round-off errors caused nonsense values in the table

• For $$h < 10^{-8}$$ the results are totally wrong!
• We would expect better approximations as $$h$$ gets smaller
• Problem 1: for small $$h$$ we subtract numbers of approx equal size and this gives rise to round-off errors
• Problem 2: for small $$h$$ the round-off errors are multiplied by a big number
• Remedy: use float variables with more digits
• Python has a (slow) float variable (decimal.Decimal) with arbitrary number of digits
• Using 25 digits gives accurate results for $$h \leq 10^{-13}$$
• Is this really a problem? Quite seldom - other uncertainies in input data to a mathematical computation makes it usual to have (e.g.) $$10^{-2}\leq h \leq 10^{-6}$$

## Lambda functions for compact inline function definitions

def f(x):
return x**2 - 1


The lambda construction can define this function in one line:

f = lambda x: x**2 - 1


In general,

somefunc = lambda a1, a2, ...: some_expression


is equivalent to

def somefunc(a1, a2, ...):
return some_expression


Lambda functions can be used directly as arguments in function calls:

value = someotherfunc(lambda x, y, z: x+y+3*z, 4)


## Example on using a lambda function to save typing

Verbose standard code:

def g(t):
return t**(-6)

dgdt = diff2(g, 2)
print dgdt


More compact code with lambda:

dgdt = diff2(lambda t: t**(-6), 2)
print dgdt


## If tests for branching the flow of statements

Sometimes we want to peform different actions depending on a condition. Example: $$f(x) = \left\lbrace\begin{array}{ll} \sin x, & 0\leq x\leq \pi\\ 0, & \hbox{otherwise} \end{array}\right.$$

A Python implementation of $$f$$ needs to test on the value of $$x$$ and branch into two computations:

## The general form of if tests

if-else (the else block can be skipped):

if condition:
<block of statements, executed if condition is True>
else:
<block of statements, executed if condition is False>


Multiple if-else.

if condition1:
<block of statements>
elif condition2:
<block of statements>
elif condition3:
<block of statements>
else:
<block of statements>
<next statement>


## Example on multiple branching

A piecewisely defined function.

$$N(x) = \left\lbrace\begin{array}{ll} 0, & x < 0\\ x, & 0\leq x < 1\\ 2-x, & 1\leq x < 2\\ 0, & x \geq 2 \end{array}\right.$$

Python implementation with multiple if-else-branching.

def N(x):
if x < 0:
return 0
elif 0 <= x < 1:
return x
elif 1 <= x < 2:
return 2 - x
elif x >= 2:
return 0


## Inline if tests for shorter code

Common construction:

if condition:
variable = value1
else:
variable = value2


More compact syntax with one-line if-else:

variable = (value1 if condition else value2)


Example:

def f(x):
return (sin(x) if 0 <= x <= 2*pi else 0)


## We shall write special test functions to verify functions

def double(x):            # some function
return 2*x

def test_double():        # associated test function
"""Call double(x) to check that it works."""
x = 4                 # some chosen x value
expected = 8          # expected result from double(x)
computed = double(x)
success = computed == expected  # boolean value: test passed?
msg = 'computed %s, expected %s' % (computed, expected)
assert success, msg


Rules for test functions:

• name begins with test_
• no arguments
• must have an assert success statement, where success is True if the test passed and False otherwise (assert success, msg prints msg on failure)
The optional msg parameter writes a message if the test fails.

## Test functions with many tests

def double(x):            # some function
return 2*x

def test_double():        # associated test function
tol = 1E-14           # tolerance for float comparison
x_values =        [3, 7,  -2, 0, 4.5, 'hello']
expected_values = [6, 14, -4, 0,   9, 'hellohello']
for x, expected in zip(x_values, expected_values):
computed = double(x)
msg = '%s != %s' % (computed, expected)
assert abs(expected - computed) < tol, msg


A test function will run silently if all tests pass. If one test above fails, assert will raise an AssertionError.

## Why write test functions according to these rules?

• Easy to recognize where functions are verified
• Test frameworks, like nose and pytest, can automatically run all your test functions (in a folder tree) and report if any bugs have sneaked in
• This is a very well established standard

Terminal> py.test -s .
Terminal> nosetests -s .


We recommend py.test - it has superior output.

Unit tests.

A test function as test_double() is often referred to as a unit test since it tests a small unit (function) of a program. When all unit tests work, the whole program is supposed to work.

• Many find test functions to be a difficult topic
• The idea is simple: make problem where you know the answer, call the function, compare with the known answer
• Just write some test functions and it will be easy
• The fact that a successful test function runs silently is annoying - can (during development) be convenient to insert some print statements so you realize that the statements are run

## Summary of if tests and functions

If tests:

if x < 0:
value = -1
elif x >= 0 and x <= 1:
value = x
else:
value = 1


User-defined functions:

def quadratic_polynomial(x, a, b, c):
value = a*x*x + b*x + c
derivative = 2*a*x + b
return value, derivative

# function call:
x = 1
p, dp = quadratic_polynomial(x, 2, 0.5, 1)
p, dp = quadratic_polynomial(x=x, a=-4, b=0.5, c=0)


Positional arguments must appear before keyword arguments:

def f(x, A=1, a=1, w=pi):
return A*exp(-a*x)*sin(w*x)


## A summarizing example for Chapter 3; problem

An integral $$\int_a^b f(x)dx$$ can be approximated by Simpson's rule: \begin{align*} \int_a^b f(x)dx \approx {b-a\over 3n}\biggl( & f(a) + f(b) + 4\sum_{i=1}^{n/2} f(a + (2i-1)h)\\ & + 2\sum_{i=1}^{n/2-1} f(a+2ih)\biggr) \end{align*}

Problem: make a function Simpson(f, a, b, n=500) for computing an integral of f(x) by Simpson's rule. Call Simpson(...) for $${3\over2}\int_0^\pi\sin^3x dx$$ (exact value: 2) for $$n=2,6,12,100,500$$.

## The program: function for computing the formula

def Simpson(f, a, b, n=500):
"""
Return the approximation of the integral of f
from a to b using Simpson's rule with n intervals.
"""

h = (b - a)/float(n)

sum1 = 0
for i in range(1, n/2 + 1):
sum1 += f(a + (2*i-1)*h)

sum2 = 0
for i in range(1, n/2):
sum2 += f(a + 2*i*h)

integral = (b-a)/(3*n)*(f(a) + f(b) + 4*sum1 + 2*sum2)
return integral


## The program: function, now with test for possible errors

def Simpson(f, a, b, n=500):

if a > b:
print 'Error: a=%g > b=%g' % (a, b)
return None

# Check that n is even
if n % 2 != 0:
print 'Error: n=%d is not an even integer!' % n
n = n+1  # make n even

# as before...
...
return integral


## The program: application (and main program)

def h(x):
return (3./2)*sin(x)**3

from math import sin, pi

def application():
print 'Integral of 1.5*sin^3 from 0 to pi:'
for n in 2, 6, 12, 100, 500:
approx = Simpson(h, 0, pi, n)
print 'n=%3d, approx=%18.15f, error=%9.2E' % \
(n, approx, 2-approx)

application()


## The program: verification (with test function)

Property of Simpson's rule: 2nd degree polynomials are integrated exactly!

def test_Simpson():      # rule: no arguments
"""Check that quadratic functions are integrated exactly."""
a = 1.5
b = 2.0
n = 8
g = lambda x: 3*x**2 - 7*x + 2.5       # test integrand
G = lambda x: x**3 - 3.5*x**2 + 2.5*x  # integral of g
exact = G(b) - G(a)
approx = Simpson(g, a, b, n)
success = abs(exact - approx) < 1E-14  # tolerance for floats
msg = 'exact=%g, approx=%g' % (exact, approx)
assert success, msg


Can either call test_Simpson() or run nose or pytest:

Terminal> nosetests -s Simpson.py
Terminal> py.test -s Simpson.py
...
Ran 1 test in 0.005s

OK