$$ \newcommand{\tp}{\thinspace .} $$

This chapter is taken from the book A Primer on Scientific Programming with Python by H. P. Langtangen, 5th edition, Springer, 2016.

Evaluating standard mathematical functions

Mathematical formulas frequently involve functions such as sin, cos, tan, sinh, cosh, exp, log, etc. On a pocket calculator you have special buttons for such functions. Similarly, in a program you also have ready-made functionality for evaluating these types of mathematical functions. One could in principle write one's own program for evaluating, e.g., the \( \sin(x) \) function, but how to do this in an efficient way is a non-trivial topic. Experts have worked on this problem for decades and implemented their best recipes in pieces of software that we should reuse. This section tells you how to reach sin, cos, and similar functions in a Python context.

Example: Using the square root function

Problem

Consider the formula for the height \( y \) of a ball in vertical motion, with initial upward velocity \( v_0 \): $$ \begin{equation*} y_c = v_0t - \frac{1}{2}gt^2,\end{equation*} $$ where \( g \) is the acceleration of gravity and \( t \) is time. We now ask the question: How long time does it take for the ball to reach the height \( y_c \)? The answer is straightforward to derive. When \( y=y_c \) we have $$ \begin{equation*} y_c = v_0t - \frac{1}{2}gt^2 \tp \end{equation*} $$ We recognize that this equation is a quadratic equation, which we must solve with respect to \( t \). Rearranging, $$ \begin{equation*} \frac{1}{2}gt^2 - v_0t + y_c = 0,\end{equation*} $$ and using the well-known formula for the two solutions of a quadratic equation, we find $$ \begin{equation} t_1 = \left({v_0 - \sqrt{v_0^2 - 2gy_c}}\right)/g,\quad t_2 = \left({v_0 + \sqrt{v_0^2 - 2gy_c}}\right)/g\tp \tag{4} \end{equation} $$ There are two solutions because the ball reaches the height \( y_c \) on its way up \( (t=t_1 \)) and on its way down (\( t=t_2 > t_1 \)).

The program

To evaluate the expressions for \( t_1 \) and \( t_2 \) from (4) in a computer program, we need access to the square root function. In Python, the square root function and lots of other mathematical functions, such as sin, cos, sinh, exp, and log, are available in a module called math. We must first import the module before we can use it, that is, we must write import math. Thereafter, to take the square root of a variable a, we can write math.sqrt(a). This is demonstrated in a program for computing \( t_1 \) and \( t_2 \): 

v0 = 5
g = 9.81
yc = 0.2
import math
t1 = (v0 - math.sqrt(v0**2 - 2*g*yc))/g
t2 = (v0 + math.sqrt(v0**2 - 2*g*yc))/g
print 'At t=%g s and %g s, the height is %g m.' % (t1, t2, yc)
The output from this program becomes 

At t=0.0417064 s and 0.977662 s, the height is 0.2 m.
You can find the program as the file ball_yc.py in the src/formulas folder.

Two ways of importing a module

The standard way to import a module, say math, is to write 

import math
and then access individual functions in the module with the module name as prefix as in 

x = math.sqrt(y)
People working with mathematical functions often find math.sqrt(y) less pleasing than just sqrt(y). Fortunately, there is an alternative import syntax that allows us to skip the module name prefix. This alternative syntax has the form from module import function. A specific example is 

from math import sqrt
Now we can work with sqrt directly, without the math. prefix. More than one function can be imported: 

from math import sqrt, exp, log, sin
Sometimes one just writes 

from math import *
to import all functions in the math module. This includes sin, cos, tan, asin, acos, atan, sinh, cosh, tanh, exp, log (base \( e \)), log10 (base 10), sqrt, as well as the famous numbers e and pi. Importing all functions from a module, using the asterisk (*) syntax, is convenient, but this may result in a lot of extra names in the program that are not used. It is in general recommended not to import more functions than those that are really used in the program. Nevertheless, the convenience of the compact from math import * syntax occasionally wins over the general recommendation among practitioners - and in this document.

With a from math import sqrt statement we can write the formulas for the roots in a more pleasing way: 

t1 = (v0 - sqrt(v0**2 - 2*g*yc))/g
t2 = (v0 + sqrt(v0**2 - 2*g*yc))/g

Import with new names

Imported modules and functions can be given new names in the import statement, e.g., 

import math as m
# m is now the name of the math module
v = m.sin(m.pi)

from math import log as ln
v = ln(5)

from math import sin as s, cos as c, log as ln
v = s(x)*c(x) + ln(x)
In Python, everything is an object, and variables refer to objects, so new variables may refer to modules and functions as well as numbers and strings. The examples above on new names can also be coded by introducing new variables explicitly: 

m = math
ln = m.log
s = m.sin
c = m.cos

Example: Computing with \( \sinh x \)

Our next examples involve calling some more mathematical functions from the math module. We look at the definition of the \( \sinh (x) \) function: $$ \begin{equation} \sinh (x) = \frac{1}{2}\left(e^{x} - e^{-x}\right)\tp \tag{5} \end{equation} $$ We can evaluate \( \sinh(x) \) in three ways: i) by calling math.sinh, ii) by computing the right-hand side of (5), using math.exp, or iii) by computing the right-hand side of (5) with the aid of the power expressions math.e**x and math.e**(-x). A program doing these three alternative calculations is found in the file 3sinh.py. The core of the program looks like this: 

from math import sinh, exp, e, pi
x = 2*pi
r1 = sinh(x)
r2 = 0.5*(exp(x) - exp(-x))
r3 = 0.5*(e**x - e**(-x))
print r1, r2, r3
The output from the program shows that all three computations give identical results: 

267.744894041 267.744894041 267.744894041

A first glimpse of rounding errors

The previous example computes a function in three different yet mathematically equivalent ways, and the output from the print statement shows that the three resulting numbers are equal. Nevertheless, this is not the whole story. Let us try to print out r1, r2, r3 with 16 decimals: 

print '%.16f %.16f %.16f' % (r1,r2,r3)
This statement leads to the output 

267.7448940410164369 267.7448940410164369 267.7448940410163232
Now r1 and r2 are equal, but r3 is different! Why is this so?

Our program computes with real numbers, and real numbers need in general an infinite number of decimals to be represented exactly. The computer truncates the sequence of decimals because the storage is finite. In fact, it is quite standard to keep only 17 digits in a real number on a computer. Exactly how this truncation is done is not explained in this document, but you read more on Wikipedia. For now the purpose is to notify the reader that real numbers on a computer often have a small error. Only a few real numbers can be represented exactly, the rest of the real numbers are only approximations.

For this reason, most arithmetic operations involve inaccurate real numbers, resulting in inaccurate calculations. Think of the following two calculations: \( 1/49\cdot 49 \) and \( 1/51\cdot 51 \). Both expressions are identical to 1, but when we perform the calculations in Python, 

print '%.16f %.16f' % (1/49.0*49, 1/51.0*51)
the result becomes 

0.9999999999999999 1.0000000000000000
The reason why we do not get exactly 1.0 as answer in the first case is because 1/49 is not correctly represented in the computer. Also 1/51 has an inexact representation, but the error does not propagate to the final answer.

To summarize, errors in floating-point numbers may propagate through mathematical calculations and result in answers that are only approximations to the exact underlying mathematical values. The errors in the answers are commonly known as rounding errors. As soon as you use Python interactively as explained in the next section, you will encounter rounding errors quite often.

Python has a special module decimal and the SymPy package has an alternative module mpmath, which allow real numbers to be represented with adjustable accuracy so that rounding errors can be made as small as desired (an example appears the document Functions and branching [5]). However, we will hardly use such modules because approximations implied by many mathematical methods applied throughout this document normally lead to (much) larger errors than those caused by rounding.