How to do mathematical analysis

9.1 Introduction

This how-to will demonstrate how to numerically find roots, solve equations, and perform integrals in Python. The scipy package has been chosen for these applications since it has proven very robust.

9.2 Finding roots

Assume you have a function:

$f(x) = (x-2)(x+1)$

and want to find its roots (values $x$ for which $f(x)=0$ ). That can be done using the scipy function, fsolve, that takes two input arguments. The first argument is the function for which you are seeking the roots, and the second argument is a value at which you want fsolve to start its search for a root. Searching for a root for the function $f$ starting at $x=3$ would thus look like this:

from scipy.optimize import fsolve
f = lambda x: (x-2)*(x+1)
fsolve(f,3)

array([2.])

while search for a root starting at $x=-2$ would look like this:

fsolve(f,-2)

array([-1.])

If you like, you could search for several roots like this:

x1 = fsolve(f,3)
x2 = fsolve(f,-2)
x3 = fsolve(f,0)
(x1, x2, x3)

(or in any other more complicated way, e.g. involving a loop construction) that would give this output:

(array([2.]), array([-1.]), array([-1.]))

where one of the roots is found twice.

You may convert from NumPy-arrays to floats with a list comprehension: [float(x) for x in (x1, x2, x3)], thereby getting: [2.0, -1.0, -1.0]. The two values, $2$ and $-1$ , are the $x$ -values for which $f$ gives zero. If we plot the entire function and indicate those two place with points, $(2, 0)$ and $(-1, 0)$ , one can see that these $x$ -values are indeed the ones for which $f(x)$ is zero:

Figure 1

9.3 Solving equations

Assume that you wish to solve this equation numerically:

$\cos(x) = x$

A simple trick lets you solve it as if you were looking for a root of a function. Simply subtract the right-hand side from the left-hand side:

$\cos(x) - x = 0.$

So if $f(x)=\cos(x) - x$ , you will be looking for the roots of $f(x)$ , which looks like this in Python:

from scipy.optimize import fsolve
import numpy as np
f = lambda x : np.cos(x) - x
fsolve(f,0)

array([0.73908513])

Figure 2

9.4 Minimization

Consider this function:

$f(x) = \frac{1}{10}x^2+\cos^2(x)$

It has several local minima. The $x$ -value for one of them can be found with a call of scipy.optimize's fmin function like this:

import scipy.optimize
import numpy as np

f = lambda x: 1/10*x**2 + np.cos(x)**2
xmin = scipy.optimize.fmin(func=f, x0=1)
xmin,f(xmin)

        Optimization terminated successfully.
             Current function value: 0.224167
             Iterations: 14
             Function evaluations: 28
Out[1]: (array([1.42617188]), array([0.22416743]))

where x0=1 instructs fmin to start the search near $x=1$ .

If you are annoyed with the NumPy array format of the output, you may use a list comprehension when inspecting the result:

[float(a) for a in [xmin, f(xmin)]]

[1.426171875000001, 0.22416743026027897]

Figure 3

9.5 Integration

Numerical integration of a function $f(x)$ from $x=a$ to $x=b$ is done in Python with the function quad from the scipy.integrate module. quad takes three arguments:

a function of one variable
the lower limit of integration
the upper limit of integration

As a mnemonic, you may write:

$\mathrm{quad(f,a,b)}=\int_a^b f(x)\,dx.$

For an example, consider that you would like to evaluate this integral numerically:

$\int_0^2x^3dx=\Bigl[\frac{1}{4}x^4\Bigr]_0^2=4$

Then in Python it would look like this:

from scipy.integrate import quad
f = lambda x: pow(x,3)
quad(f, 0, 2)

that gives:

(4.0, 4.440892098500626e-14)

where 4.0 is the result and 4.440892098500626e-14 is an estimate of the absolute error in the result.

If you are only interested in the result of the integration and not the absolute error, you may simply ask for the first entry in the result by writing:

quad(f, 0, 2)[0]

whereby the output now becomes: 4.0.

9.6 Length of curve

Integrals are often interpreted as the area (with sign) of the region between a curve and the horizontal axis. Depending on the integral, it may also have other interpretations. Consider for instance that a curve is given by the parametric expression $\bigl(x(\theta),y(\theta)\bigr)$ , where $\theta$ is a parameter, then the length of the curve from $\left(x(\theta_0),y(\theta_0)\right)$ to $\Bigl(x(\theta_1),y(\theta_1)\Bigr)$ is given by:

$s=\int_{\theta_0}^{\theta_1}\sqrt{\Bigl(\frac{dx}{d\theta}\Bigr)^2+\Bigl(\frac{dy}{d\theta}\Bigr)^2}d\theta$

for $\theta_0 < \theta_1$ .

Consider as an example, a circle with radius 1:

$\begin{array}{lcl} x&=&\cos(\theta)\\ y&=&\sin(\theta) \end{array}$

Before we can use the curve length formula, we need to evaluate analytically the two derivatives:

$\begin{array}{lcl} \displaystyle \frac{dx}{d\theta}&=&-\sin(\theta)\\\\ \displaystyle \frac{dy}{d\theta}&=&\cos(\theta) \end{array}$

We are now ready to evaluate the curve length:

from scipy.integrate import quad
dxdtheta = lambda theta: -np.sin(theta)
dydtheta = lambda theta: np.cos(theta)
quad(lambda theta: np.sqrt(dxdtheta(theta)**2 + dydtheta(theta)**2),0,2*np.pi)[0]

That gives the expected result: 6.283185307179586 which is the numerical representation of $2\pi$ .

9.7 Solving integral equations I

Consider the Normal distribution with standard deviation $1$ and mean $0$ :

$f(x)=\frac{1}{\sqrt{2\pi}}\exp\bigl(-\frac{x^2}{2}\bigr).$

and ask the question: What should $x_0$ be for half of all observations to lie in the interval from $-x_0$ to $x_0$ ? Expressed mathematically, we want to solve:

$\int_{-x_0}^{x_0}f(x)dx=\frac{1}{2}.$

We can evaluate the integral with quad for different values of $x_0$ :

f = lambda x: 1/np.sqrt(2*np.pi) * np.exp(-x**2/2)
quad(f,-0.5,0.5)

(0.3829249225480263, 4.251320657156236e-15)

quad(f,-0.75,0.75)[0]      # the [0] is added so that the output becomes the value of the integral

0.5467452952462637

Thus, by now we know that $\displaystyle 0.3829<\int_{-x_0}^{x_0}f(x)dx<0.5467$ , if $0.5<x_0<0.75$ , but this is clearly an inefficient way of determining $x_0$ . A simple trick, however, turns the problem into finding a root of a new function, $g$ :

$g(x_0)=\frac{1}{2}-\int_{-x_0}^{x_0}f(x)dx=0$

whereby we can use fsolve:

import numpy as np
from scipy.integrate import quad
from scipy.optimize import fsolve

f = lambda x: 1/np.sqrt(2*np.pi) * np.exp(-x**2/2)
x0sol = fsolve(lambda x0: 0.5 - quad(f,-x0,x0)[0],1)[0]
print('x0sol:',x0sol)

which gives:

x0sol: 0.6744897501960817

Figure 4 The area fulfilling $\displaystyle \int_{-0.6745}^{0.6745}f(x)dx=\frac{1}{2}$ .

9.8 Solving integral equations II

We have already solved an integral equation in which the integration bounds were the unknowns. Another example of an integral equation would be one, where a parameter in the integrand is unknown. Consider for instance the Normal distribution function with standard deviation $\sigma$ and mean $0$ :

$f(x,\sigma)=\frac{1}{\sqrt{2\pi}\sigma}\exp\bigl(-\frac{x^2}{2\sigma^2}\bigr).$

In the following integral equation involving $f(x,\sigma)$ , the parameter $\sigma$ is unknown:

$\int_{-1}^{1}f(x,\sigma)dx=\frac{1}{2}.$

and we are effectively asking the question: what must $\sigma$ be for the integral on the left hand side to evaluate to the value on the right hand side?

Coding functions of two variables is discussed elsewhere. For the normal distribution, the code may look like this:

f = lambda x,sigma: 1/np.sqrt(2*np.pi)/sigma * np.exp(-np.power(x,2)/2/np.power(sigma,2))

We may then of course do some trial and error in an attempt to find the $\sigma$ that fulfill the equation:

print('For sigma=1   ',quad(lambda x: f(x,1),-1,1)[0])
print('For sigma=1.25',quad(lambda x: f(x,1.25),-1,1)[0])
print('For sigma=1.5 ',quad(lambda x: f(x,1.5),-1,1)[0])

(recall, the [0] is selecting the first value in the tuple provided by quad).

For sigma=1    0.682689492137086
For sigma=1.25 0.5762892028332066
For sigma=1.5  0.4950149249061543

and find that $\sigma\approx 1.5$ . Again, reformulating the problem as:

$g(\sigma)=\frac{1}{2}-\int_{-1}^{1}f(x,\sigma)dx=0$

the problem can be solved automatically using fsolve:

f = lambda x,sigma: 1/np.sqrt(2*np.pi)/sigma * np.exp(-np.power(x,2)/2/np.power(sigma,2))
sigma_sol = fsolve(lambda sigma: 0.5 - quad(lambda x: f(x,sigma),-1,1)[0],1)[0]
sigma_sol

1.4826022185056023

Thus for $\sigma=1.4826$ , half of all observations lie between $-x_0$ and $x_0$ for $x_0=1$ . This can be confirmed by inserting the found value for $\sigma$ in the integral:

quad(lambda x: f(x,sigma_sol),-1,1)

(0.4999999999999999, 5.551115123125781e-15)