This chapter is taken from the book A Primer on Scientific Programming with Python by H. P. Langtangen, 5th edition, Springer, 2016.
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.
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 \)).
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.
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
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
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
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.