Link to booklet is copied to clipboard!

Saved to booklet!

Removed from booklet!

Started capturing to booklet.

Stopped capturing to booklet.

Preferring latest booklet.

No longer prefers latest booklet.

Entered edit mode.

Exited edit mode.

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.
Assume you have a function:

and want to find its roots (values for which ). 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 starting at 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 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, and , are the -values for which gives zero. If we plot
the entire function and indicate those two place with points, and , one can
see that these -values are indeed the ones for which is zero:
Assume that you wish to solve this equation numerically:

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:

So if , you will be looking for the roots of , 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])

Consider this function:

It has several local minima. The -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
If you are annoyed with the *NumPy array* format of the output, you may use a list comprehension when inspecting the result:

`x0=1`

instructs `fmin`

to start the search near .
[float(a) for a in [xmin, f(xmin)]]

[1.426171875000001, 0.22416743026027897]

Numerical integration of a function from to
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

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

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
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:

`4.0`

is the result
and `4.440892098500626e-14`

is an estimate of
the absolute error in the result.
quad(f, 0, 2)[0]

whereby the output now becomes:

`4.0`

.
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 , where is a parameter, then the
length of the curve from to is given by:

for .
Consider as an example, a circle with radius 1:

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

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 .
Consider the Normal distribution with standard deviation and mean :

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

We can evaluate the integral with

`quad`

for different values of :
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 ,
if , but this is clearly an inefficient way of determining . A simple trick, however, turns the problem
into finding a root of a new function, :

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

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 and mean :

In the following integral equation involving , the parameter
is unknown:

and we are effectively asking the question: what must 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 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 . Again, reformulating the problem as:

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 , half of all observations lie between and for . This can be confirmed by
inserting the found value for in the integral:

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

(0.4999999999999999, 5.551115123125781e-15)

Once deleted this booklet is gone forever!

Choose which booklet to go to: