Lesson 4. Functions

2023-01-23

DIY Functions

We have already seen quite a few functions in previous lessons, but so far they have all been made available for us as part of the core Python modules or from a third party module like NumPy. How can we create our own functions?

Objectives

Our aim in this lesson is to learn how to write simple Python functions that evaluate a mathematical expression or expressions. In other words, we are interested in the kind of function we commonly meet in mathematics that maps a subset of \(\mathbf{R}^n\) into \(\mathbf{R}^m\).

By the end of this lesson, you should know

the syntax for defining a Python function;
how arguments are passed to a function;
how values are returned from function;
how to access a function defined in a file;
the meaning of ‘scope’ as it applies to the variables in a function body;

Defining a Function

Suppose that we want to solve the nonlinear equation \(\sin(4x)=\log(x)\).
Algorithms for computing the solution typically require that you first write the equation in the standard form \(f(x)=0\), so in our case we should define \[ f(x) = \sin(4x) - \log(x). \] Use the editor pane in Spyder to type

from numpy import sin, log

def f(x):
    return sin(4*x) - log(x)

and then save to a file myfuncs.py in your working folder (as displayed in the upper right of the Sypder window).

Here, we use the def keyword to indicate that we are defining a function.
After def is a space followed by the name of the function, in this case f, and then a comma-separated list of input arguments, enclosed in parentheses. In this case, there is only one argument so we don’t need any commas. Note the colon at the end of the line, which is an essential part of the syntax.

On the next line we indent four spaces and write the body of the function.
In this example, the body consists only of the single line beginning with the keyword return. The expression following this keyword is evaluated and forms the return value of the function. You probably noticed that the editor automatically indents the code correctly.

At the IPython prompt in the lower right Spyder pane, type

from myfuncs import f
f(0.5)
f(1.0)

You should find that \(f(0.5)\approx 1.6024>0\) and \(f(1.0)\approx-0.7568<0\), so the function must have at least one zero in the interval \((0.5,1.0)\).

The root_scalar function is available from a submodule of SciPy called ‘optimize’. At the IPython prompt, type

from scipy.optimize import root_scalar
?root_scalar

to see the docstring, and then type

sol = root_scalar(f, method='brentq', bracket=(0.5,1.0))

We will discuss the method and bracket arguments below, so don’t worry about them for now. The root_scalar function returns a RootResults object that we have assigned to the variable sol. The numerical value of the root is accessed as sol.root. As a check, we can evaluate the residual, that is, the value of \(f\) at the computed root and observe how close it gets to zero:

f(sol.root)

You should see that the residual is close to the machine epsilon \(2.22\times10^{-16}\).

In this example, to solve our equation we had to write a function f and pass it as an argument to another (much more complicated) function (written by experts). These steps are typical of how Python is used to solve many mathematical problems.

Argument Passing

A Python function can have any number of arguments and any number of return variables. Consider a function \(g:\mathbf{R}^2\to\mathbf{R}^2\) defined by \[ g(x,y)=(2x-y, x^2+\cos(xy)). \] In your myfuncs.py file, type

from numpy import cos

def g(x, y):
    return 2*x-y, x**2 + cos(x*y)

and then save. At the IPython prompt, type

from myfuncs import g
g(1, 3)

and observe that the return value is (-1, 0.010007503399554585). We could also call the function by doing

a, b = g(1, 3)

so that a gets assigned the value -1, and b the value 0.010007503399554585.

We will use the terms dummy arguments and actual arguments to distinguish between the arguments in the definition of a function (here, x and y in our file myfuncs.py) and the arguments supplied when we call the function (here, 1 and 3 in the IPython pane). When the function is called, the kth dummy argument is assigned the value of the kth actual argument, then the code in the body of the function is executed, and any values in the return statement get returned when the function exits.

Note that the names of the actual arguments are irrelevant. We could (if feeling perverse) do

x = 3
y = 1
g(y, x)

and the result would be the same as typing g(1, 3). However, to avoid confusion it is better to do

x = 1
y = 3
g(x, y)

so that the names of the actual arguments match those of the dummy arguments.

Keyword and Optional Arguments

Instead of relying on the order of a function’s arguments, we can pass values by keyword. Continuing with the function g defined above, instead of calling g(1,3) we could type

g(x=1, y=3)

The order of keyword arguments does not matter, so

g(y=3, x=1)

would give the same result.

The function call to root_scalar illustrates another feature of Python functions. The docstring shows that root_scalar really has 12 arguments, but all except the first are optional because they have been assigned a default value. In our call, we changed two of the optional arguments, method and bracket, passing the values by keyword. Here, bracket refers to an interval that brackets the root (with f taking opposite signs at the ends of this interval). Setting method='brentq' tells root_scalar to use an algorithm invented by Prof. Richard Brent from the ANU; the docstring lists seven other possible choices.

Python requires that all non-keyword arguments (also called positional arguments) in a function call must precede any keyword arguments. Thus, in our call to root_scalar the positional argument f comes before the two keyword arguments method and bracket.

Use of an optional argument is convenient if there is a natural default value that will usually be required, but you might sometimes want a different value.
A simple example is a function that computes the normal probability density function with mean mu and standard deviation sigma.

def normal_distribution(x, mu=0.0, sigma=1.0):
    return exp(-((x-mu)/sigma)**2/2) / (sigma*sqrt(2*pi))

Calling normal_distribution(x) returns the value of the standard normal density at x, but we could call normal_distribution(x, 3.0, 0.5) for a normal distribution with mean \(3.0\) and standard deviation \(0.5\). We could also call normal_distribution(x, mu=3.0, sigma=0.5) which might be helpful for someone reading the code.

Local Variables

In the examples above, we wrote functions in which the body consists of just a single line, namely the return statement, but often the body will consist of several steps. For example, recall the problem of computing the solutions of a quadratic equation \(ax^2+bx+c=0\), and add the following in your file myfuncs.py.

from numpy import sqrt

def solve_quadratic(a, b, c):
    sqrt_dscr = sqrt(b**2 - 4*a*c)
    x_plus  = (-b + sqrt_dscr ) / ( 2 * a )
    x_minus = (-b - sqrt_dscr ) / ( 2 * a )
    return x_plus, x_minus

In this example, sqrt_dscr, x_plus and x_minus are said to be local variables in the scoping unit of the function solve_quadratic. The practical significance is that local variables cannot be referenced outside their smallest enclosing scoping unit. Thus, doing

from myfuncs import solve_quadratic

xp, xm = solve_quadratic(1, -5, 2)
sqrt_dscr

will raise a NameError: name 'sqrt_dscr' is not defined. The error is not with the function solve_quadratic, but rather occurs because there is no variable called sqrt_dscr defined in the IPython shell. Similarly, you cannot reference the local variables x_plus and x_minus, although their values can be accessed from xp and xm since the function call effectively results in the assignments xp = x_plus and xm = x_minus.

Limiting the scope of variables in this way is necessary, because otherwise a function call might cause unexpected side effects. For example, you might be using a variable that has the same name as a local variable in the body of the function.

The scoping rules allow a function body to access variables defined in an enclosing scoping unit, as with a in the following example.

a = 5
def h(x):
    return x**2 + a

Indentation

Python is generally relaxed about white space, so that x=y+2 is treated no differently from x = y + 2. However, the start of each line must be indented correctly. For example, doing

def solve_quadratic(a, b, c):
    sqrt_dscr = sqrt(b**2 - 4*a*c)
    x_plus  = (-b + sqrt_dscr ) / ( 2 * a )
   x_minus = (-b - sqrt_dscr ) / ( 2 * a )
    return x_plus, x_minus

raises an IndentationError: unindent does not match any outer indentation level. Fortunately, the Spyder editor takes care of this for you, but care is needed if you copy and paste several lines.

You should also avoid having a line extend beyond 80 columns of text. If you have to evaluate a long, complicated expression that does not fit easily on one line, you can spread it over several lines using a pair of parentheses. For example,

the_answer = ( sin(x+7*y) - exp(-2*z) * ( log(2*x) - 7 ) + tan(3*x)
               - (x**2 + y**2) / ( 3 + sin(pi*z) ) )

Here, the first opening parenthesis ( on the right-hand side has its corresponding closing parenthesis ) at the end of the second line.
Consequently, instead of stopping at the end of the first line, the Python intepretor continues reading to the end of the second line.

Object Methods

A Python object will often have methods associated with it. For example, consider the list

l = [1, 0, -6, 3, 8, 0, 8, -5]

The list type has a sort method, that we can call as follows:

l.sort()

This method updates the list so that its entries are sorted in increasing order; in this case, l will update to [-6, -5, 0, 0, 1, 3, 8, 8]. In IPython, You can use tab completion to find out the names of an object’s methods: type l. followed by the tab key, that is, l followed by a period . and then tab. Notice the method reverse, that reverses the order of the entries. For example, after the command

l = [1, 4, 2, 5, 3, 6]
l.reverse()

the list l will be [6, 3, 5, 2, 4, 1].

Arrays have a large number of methods, some of which are listed below, using as an example

a = array([2, 4, 6, 8])

Method	Action	Example Result
cumprod()	cumulative product	array([2, 8, 48, 384])
cumsum()	cumulative sum	array([2, 6, 12, 20])
max()	maximum entry	8
mean()	arithmetic mean	5.0
min()	minimum entry	2
prod()	product	384
std()	standard deviation	2.23606797749979
sum()	sum	20

Another important array method is copy, which returns a copy of the array: type the commands

a = array([1, 2, 3, 4])
b = a.copy()
c = a
a[2] = -1

You will find that c[2] = -1 but b[2] = 3, because a and c reference the same storage, but b references its own copy of (the original) a. Lists also have a copy method.

Summary

From this lesson, you have learned

the syntax for defining a function;
how positional and keyword arguments are passed to a function;
how values are returned from a function;
about the use of optional arguments;
how the function body is a scoping unit for local variables;
about the use of object methods.