Notebook

Scipy -- Library of scientific algorithms for Python¶

Based on lecture at http://github.com/jrjohansson/scientific-python-lectures.

The SciPy framework builds on top of the low-level NumPy framework for multidimensional arrays, and provides a large number of higher-level scientific algorithms. Some of the topics that SciPy covers are:

Special functions (scipy.special)
Integration (scipy.integrate)
Optimization (scipy.optimize)
Interpolation (scipy.interpolate)
Fourier Transforms (scipy.fftpack)
Signal Processing (scipy.signal)
Linear Algebra (scipy.linalg)
Sparse Eigenvalue Problems (scipy.sparse)
Statistics (scipy.stats)
Multi-dimensional image processing (scipy.ndimage)
File IO (scipy.io)

Each of these submodules provides a number of functions and classes that can be used to solve problems in their respective topics.

In this lecture we will look at how to use some of these subpackages.

To access the SciPy package in a Python program, we start by importing everything from the scipy module.

WARNING: In the new version of python many functionalities are now moved from scipy to numpy, but they are still available in scipy and a deprecated warning is displayed. The work-around is to first import functions from scipy and after that from numpy, to overwrite scipy functions with the same name.

In [1]:

from scipy import *
from numpy import *

If we only need to use part of the SciPy framework we can selectively include only those modules we are interested in. For example, to include the linear algebra package under the name la, we can do:

In [2]:

import scipy.linalg as la

Special functions¶

A large number of mathematical special functions are important for many computional physics problems. SciPy provides implementations of a very extensive set of special functions. For details, see the list of functions in the reference documention at http://docs.scipy.org/doc/scipy/reference/special.html#module-scipy.special.

To demonstrate the typical usage of special functions we will look in more detail at the Bessel functions:

In [3]:

#
# The scipy.special module includes a large number of Bessel-functions
# Here we will use the functions jn and yn, which are the Bessel functions 
# of the first and second kind and real-valued order. We also include the 
# function jn_zeros and yn_zeros that gives the zeroes of the functions jn
# and yn.
#
from scipy.special import jn, yn, jn_zeros, yn_zeros

In [4]:

n = 0    # order
x = 0.0  # x-point

# Bessel function of first kind
print("J_%d(%f) = %f" % (n, x, jn(n, x)))

x = 1.0
# Bessel function of second kind
print("Y_%d(%f) = %f" % (n, x, yn(n, x)))

J_0(0.000000) = 1.000000
Y_0(1.000000) = 0.088257

In [5]:

%matplotlib inline
import matplotlib.pyplot as plt
from IPython.display import Image

In [6]:

x = linspace(0, 10, 100)


for n in range(4):
    plt.plot(x, jn(n, x), label=('$J_%d(x)$' % n))
plt.legend(loc='best');

Integration¶

Numeric integration with general purpose quadrature¶

Numerical evaluation of a function of the type

$\displaystyle \int_a^b f(x) dx$

is called numerical quadrature, or simply quadature. SciPy provides a series of functions for different kind of quadrature, for example the quad, dblquad and tplquad for single, double and triple integrals, respectively.

In [7]:

from scipy.integrate import quad, dblquad, tplquad

The quad function takes a large number of optional arguments, which can be used to fine-tune the behaviour of the function (try help(quad) for details).

The basic usage is as follows:

In [8]:

x_lower, x_upper = 0.0, 1.0
val, abserr = quad( lambda x: x**2, x_lower, x_upper)
print('Value=', val, 'error=', abserr)

Value= 0.33333333333333337 error= 3.700743415417189e-15

In [9]:

L=10*pi
x=linspace(1e-12,L)
plt.plot(x,sin(x)/x)
plt.grid()

In [10]:

y1=quad(lambda x: sin(x)/x, 0,L)
from scipy import special
y2=special.sici(L)
print(y1[0]-y2[0])

4.440892098500626e-16

If we need to pass extra arguments to integrand function we can use the args keyword argument. Let's say we want to evaluate

$f(x) = \displaystyle \int_0^x \frac{j_n(t)+j_m(t)}{t} dt$

In [11]:

f2 = lambda t,n,m: (jn(n,t)+jn(m,t))/t

def f3(t,n,m):
    return (jn(n,t)+jn(m,t))/t

In [12]:

# First specific case for x=1, n=1
# integrating  (j_1(t)+j_2(t))/t
quad(f2, 0, 1, args=(1,2))[0]

Out[12]:

0.5396292385998932

Now simpler: $f_n(x) = n \displaystyle \int_0^x \frac{j_n(t)}{t} dt$

In [13]:

xs = linspace(1e-10,30,300)

for ns in range(1,4):
    fs=[ns * quad(lambda t,n: jn(n, t)/t, 0, x, args=(ns,))[0] for x in xs]
    plt.plot(xs,fs, label='n='+str(ns))
plt.legend(loc='best')

Out[13]:

<matplotlib.legend.Legend at 0x7f7db809b820>

In [14]:

# here we can get away withouth args since function is defined inline

for ns in range(1,4):
    fs=[ns * quad(lambda t: jn(ns, t)/t, 0, x)[0] for x in xs]
    plt.plot(xs,fs, label='n='+str(ns))
plt.legend(loc='best')

Out[14]:

<matplotlib.legend.Legend at 0x7f7dc8308e20>

Higher-dimensional integration works in the same way:

Note that we pass lambda functions for the limits for the y integration, since these in general can be functions of x.

\begin{equation} \int_{x_a}^{x_b} \int_{y_a(x)}^{y_b(x)} f(x,y) dy dx \end{equation}

In [15]:

def integrand(x, y, z):
    return exp(-x**2-y**2-z**2)

x_a,x_b = 0,10
y_a,y_b = 0,10
z_a,z_b = 0,10

val1, abserr = dblquad(integrand, x_a, x_b, lambda x:y_a, lambda x:y_b, args=(0,))

val2, abserr = tplquad(integrand, x_a, x_b, lambda x:y_a, lambda x:y_b, lambda x,y:z_a, lambda x,y:z_b)

print(val1, val2, abserr)

0.7853981633974476 0.6960409996039545 1.4675161390126584e-08

Ordinary differential equations (ODEs)¶

SciPy provides two different ways to solve ODEs: An API based on the function odeint, and object-oriented API based on the class ode. Usually odeint is easier to get started with, but the ode class offers some finer level of control.

Here we will use the odeint functions. For more information about the class ode, try help(ode). It does pretty much the same thing as odeint, but in an object-oriented fashion.

To use odeint, first import it from the scipy.integrate module

In [16]:

from scipy import *
from numpy import *
from scipy.integrate import odeint, ode

A system of ODEs should be formulated in standard form, which is:

$y' = f(y, t)$

where we are searching for function $y(t)$ and where

$y = [y_1(t), y_2(t), ..., y_n(t)]$

and $f$ is some function that gives the derivatives of the function $y_i(t)$. To solve an ODE we need to know the function $f$ and an initial condition, $y(0)$.

Note that higher-order ODEs can always be written in this form by introducing new variables for the intermediate derivatives.

Once we have defined the Python function f and array y_0 (that is $f$ and $y(0)$ in the mathematical formulation), we can use the odeint function as:

y_t = odeint(f, y_0, t)

where t is and array with time-coordinates for which to solve the ODE problem. y_t is an array with one row for each point in time in t, where each column corresponds to a solution y_i(t) at that point in time.

We will see how we can implement f and y_0 in Python code in the examples below.

Example: double pendulum¶

Let's consider a physical example: The double compound pendulum, described in some detail here: http://en.wikipedia.org/wiki/Double_pendulum

In [17]:

from IPython.display import Image
Image(url='http://upload.wikimedia.org/wikipedia/commons/c/c9/Double-compound-pendulum-dimensioned.svg')

Out[17]:

The equations of motion of the pendulum are given on the wiki page:

${\dot \theta_1} = \frac{6}{m\ell^2} \frac{ 2 p_{\theta_1} - 3 \cos(\theta_1-\theta_2) p_{\theta_2}}{16 - 9 \cos^2(\theta_1-\theta_2)}$

${\dot \theta_2} = \frac{6}{m\ell^2} \frac{ 8 p_{\theta_2} - 3 \cos(\theta_1-\theta_2) p_{\theta_1}}{16 - 9 \cos^2(\theta_1-\theta_2)}.$

${\dot p_{\theta_1}} = -\frac{1}{2} m \ell^2 \left [ {\dot \theta_1} {\dot \theta_2} \sin (\theta_1-\theta_2) + 3 \frac{g}{\ell} \sin \theta_1 \right ]$

${\dot p_{\theta_2}} = -\frac{1}{2} m \ell^2 \left [ -{\dot \theta_1} {\dot \theta_2} \sin (\theta_1-\theta_2) + \frac{g}{\ell} \sin \theta_2 \right]$

To make the Python code simpler to follow, let's introduce new variable names and the vector notation: $x = [\theta_1, \theta_2, p_{\theta_1}, p_{\theta_2}]$

${\dot x_1} = \frac{6}{m\ell^2} \frac{ 2 x_3 - 3 \cos(x_1-x_2) x_4}{16 - 9 \cos^2(x_1-x_2)}$

${\dot x_2} = \frac{6}{m\ell^2} \frac{ 8 x_4 - 3 \cos(x_1-x_2) x_3}{16 - 9 \cos^2(x_1-x_2)}$

${\dot x_3} = -\frac{1}{2} m \ell^2 \left [ {\dot x_1} {\dot x_2} \sin (x_1-x_2) + 3 \frac{g}{\ell} \sin x_1 \right ]$

${\dot x_4} = -\frac{1}{2} m \ell^2 \left [ -{\dot x_1} {\dot x_2} \sin (x_1-x_2) + \frac{g}{\ell} \sin x_2 \right]$

In [18]:

help(odeint)

Help on function odeint in module scipy.integrate._odepack_py:

odeint(func, y0, t, args=(), Dfun=None, col_deriv=0, full_output=0, ml=None, mu=None, rtol=None, atol=None, tcrit=None, h0=0.0, hmax=0.0, hmin=0.0, ixpr=0, mxstep=0, mxhnil=0, mxordn=12, mxords=5, printmessg=0, tfirst=False)
    Integrate a system of ordinary differential equations.
    
    .. note:: For new code, use `scipy.integrate.solve_ivp` to solve a
              differential equation.
    
    Solve a system of ordinary differential equations using lsoda from the
    FORTRAN library odepack.
    
    Solves the initial value problem for stiff or non-stiff systems
    of first order ode-s::
    
        dy/dt = func(y, t, ...)  [or func(t, y, ...)]
    
    where y can be a vector.
    
    .. note:: By default, the required order of the first two arguments of
              `func` are in the opposite order of the arguments in the system
              definition function used by the `scipy.integrate.ode` class and
              the function `scipy.integrate.solve_ivp`. To use a function with
              the signature ``func(t, y, ...)``, the argument `tfirst` must be
              set to ``True``.
    
    Parameters
    ----------
    func : callable(y, t, ...) or callable(t, y, ...)
        Computes the derivative of y at t.
        If the signature is ``callable(t, y, ...)``, then the argument
        `tfirst` must be set ``True``.
    y0 : array
        Initial condition on y (can be a vector).
    t : array
        A sequence of time points for which to solve for y. The initial
        value point should be the first element of this sequence.
        This sequence must be monotonically increasing or monotonically
        decreasing; repeated values are allowed.
    args : tuple, optional
        Extra arguments to pass to function.
    Dfun : callable(y, t, ...) or callable(t, y, ...)
        Gradient (Jacobian) of `func`.
        If the signature is ``callable(t, y, ...)``, then the argument
        `tfirst` must be set ``True``.
    col_deriv : bool, optional
        True if `Dfun` defines derivatives down columns (faster),
        otherwise `Dfun` should define derivatives across rows.
    full_output : bool, optional
        True if to return a dictionary of optional outputs as the second output
    printmessg : bool, optional
        Whether to print the convergence message
    tfirst : bool, optional
        If True, the first two arguments of `func` (and `Dfun`, if given)
        must ``t, y`` instead of the default ``y, t``.
    
        .. versionadded:: 1.1.0
    
    Returns
    -------
    y : array, shape (len(t), len(y0))
        Array containing the value of y for each desired time in t,
        with the initial value `y0` in the first row.
    infodict : dict, only returned if full_output == True
        Dictionary containing additional output information
    
        =======  ============================================================
        key      meaning
        =======  ============================================================
        'hu'     vector of step sizes successfully used for each time step
        'tcur'   vector with the value of t reached for each time step
                 (will always be at least as large as the input times)
        'tolsf'  vector of tolerance scale factors, greater than 1.0,
                 computed when a request for too much accuracy was detected
        'tsw'    value of t at the time of the last method switch
                 (given for each time step)
        'nst'    cumulative number of time steps
        'nfe'    cumulative number of function evaluations for each time step
        'nje'    cumulative number of jacobian evaluations for each time step
        'nqu'    a vector of method orders for each successful step
        'imxer'  index of the component of largest magnitude in the
                 weighted local error vector (e / ewt) on an error return, -1
                 otherwise
        'lenrw'  the length of the double work array required
        'leniw'  the length of integer work array required
        'mused'  a vector of method indicators for each successful time step:
                 1: adams (nonstiff), 2: bdf (stiff)
        =======  ============================================================
    
    Other Parameters
    ----------------
    ml, mu : int, optional
        If either of these are not None or non-negative, then the
        Jacobian is assumed to be banded. These give the number of
        lower and upper non-zero diagonals in this banded matrix.
        For the banded case, `Dfun` should return a matrix whose
        rows contain the non-zero bands (starting with the lowest diagonal).
        Thus, the return matrix `jac` from `Dfun` should have shape
        ``(ml + mu + 1, len(y0))`` when ``ml >=0`` or ``mu >=0``.
        The data in `jac` must be stored such that ``jac[i - j + mu, j]``
        holds the derivative of the `i`th equation with respect to the `j`th
        state variable.  If `col_deriv` is True, the transpose of this
        `jac` must be returned.
    rtol, atol : float, optional
        The input parameters `rtol` and `atol` determine the error
        control performed by the solver.  The solver will control the
        vector, e, of estimated local errors in y, according to an
        inequality of the form ``max-norm of (e / ewt) <= 1``,
        where ewt is a vector of positive error weights computed as
        ``ewt = rtol * abs(y) + atol``.
        rtol and atol can be either vectors the same length as y or scalars.
        Defaults to 1.49012e-8.
    tcrit : ndarray, optional
        Vector of critical points (e.g., singularities) where integration
        care should be taken.
    h0 : float, (0: solver-determined), optional
        The step size to be attempted on the first step.
    hmax : float, (0: solver-determined), optional
        The maximum absolute step size allowed.
    hmin : float, (0: solver-determined), optional
        The minimum absolute step size allowed.
    ixpr : bool, optional
        Whether to generate extra printing at method switches.
    mxstep : int, (0: solver-determined), optional
        Maximum number of (internally defined) steps allowed for each
        integration point in t.
    mxhnil : int, (0: solver-determined), optional
        Maximum number of messages printed.
    mxordn : int, (0: solver-determined), optional
        Maximum order to be allowed for the non-stiff (Adams) method.
    mxords : int, (0: solver-determined), optional
        Maximum order to be allowed for the stiff (BDF) method.
    
    See Also
    --------
    solve_ivp : solve an initial value problem for a system of ODEs
    ode : a more object-oriented integrator based on VODE
    quad : for finding the area under a curve
    
    Examples
    --------
    The second order differential equation for the angle `theta` of a
    pendulum acted on by gravity with friction can be written::
    
        theta''(t) + b*theta'(t) + c*sin(theta(t)) = 0
    
    where `b` and `c` are positive constants, and a prime (') denotes a
    derivative. To solve this equation with `odeint`, we must first convert
    it to a system of first order equations. By defining the angular
    velocity ``omega(t) = theta'(t)``, we obtain the system::
    
        theta'(t) = omega(t)
        omega'(t) = -b*omega(t) - c*sin(theta(t))
    
    Let `y` be the vector [`theta`, `omega`]. We implement this system
    in Python as:
    
    >>> def pend(y, t, b, c):
    ...     theta, omega = y
    ...     dydt = [omega, -b*omega - c*np.sin(theta)]
    ...     return dydt
    ...
    
    We assume the constants are `b` = 0.25 and `c` = 5.0:
    
    >>> b = 0.25
    >>> c = 5.0
    
    For initial conditions, we assume the pendulum is nearly vertical
    with `theta(0)` = `pi` - 0.1, and is initially at rest, so
    `omega(0)` = 0.  Then the vector of initial conditions is
    
    >>> y0 = [np.pi - 0.1, 0.0]
    
    We will generate a solution at 101 evenly spaced samples in the interval
    0 <= `t` <= 10.  So our array of times is:
    
    >>> t = np.linspace(0, 10, 101)
    
    Call `odeint` to generate the solution. To pass the parameters
    `b` and `c` to `pend`, we give them to `odeint` using the `args`
    argument.
    
    >>> from scipy.integrate import odeint
    >>> sol = odeint(pend, y0, t, args=(b, c))
    
    The solution is an array with shape (101, 2). The first column
    is `theta(t)`, and the second is `omega(t)`. The following code
    plots both components.
    
    >>> import matplotlib.pyplot as plt
    >>> plt.plot(t, sol[:, 0], 'b', label='theta(t)')
    >>> plt.plot(t, sol[:, 1], 'g', label='omega(t)')
    >>> plt.legend(loc='best')
    >>> plt.xlabel('t')
    >>> plt.grid()
    >>> plt.show()

In [19]:

def dx(x,t):
    """
    The right-hand side of the pendulum ODE
    x=[x1,x2,x3,x4]
    """
    g, L, m = 9.82, 1., 1.
    x1,x2,x3,x4 = x
    c1 = 1/(m*L**2)
    dx1 = 6.*c1 * (2*x3-3*cos(x1-x2)*x4)/(16.-9.*cos(x1-x2)**2)
    dx2 = 6.*c1 * (8*x4-3*cos(x1-x2)*x3)/(16.-9.*cos(x1-x2)**2)
    dx3 = -0.5/c1 * (dx1*dx2 * sin(x1-x2) + 3*g/L * sin(x1))
    dx4 = -0.5/c1 * (-dx1*dx2 * sin(x1-x2)+ g/L * sin(x2))
    return array([dx1,dx2,dx3,dx4])

In [20]:

# choose an initial state
x0 = [pi/4, pi/2, 0, 0]

# time coodinate to solve the ODE for: from 0 to 10 seconds
t = linspace(0, 10, 250)

In [21]:

# solve the ODE problem
x = odeint(dx, x0, t)

In [22]:

# plot the angles as a function of time
fig, axes = plt.subplots(1,2,figsize=(12,4))
axes[0].plot(t, x[:, 0], label="theta1")
axes[0].plot(t, x[:, 1], label="theta2")
axes[0].legend(loc='best')
axes[0].set_xlabel('time')

L = 0.5
x1 =  L * sin(x[:, 0])
y1 = -L * cos(x[:, 0])

x2 = x1 + L * sin(x[:, 1])
y2 = y1 - L * cos(x[:, 1])
    
axes[1].plot(x1, y1, label="pendulum1")
axes[1].plot(x2, y2, label="pendulum2")
axes[1].set_ylim([-1, 0])
axes[1].set_xlim([-1, 1])
axes[1].legend(loc='best');

examples for subplot creation: https://towardsdatascience.com/the-many-ways-to-call-axes-in-matplotlib-2667a7b06e06

matplotlib examples: https://matplotlib.org/3.1.1/gallery/index.html

See animation in pendulum.py.

(To get it work within jupyter seems a bit challenging at the moment.)

Example: Damped harmonic oscillator¶

ODE problems are important in computational physics, so we will look at one more example: the damped harmonic oscillation. This problem is well described on the wiki page: http://en.wikipedia.org/wiki/Damping

The equation of motion for the damped oscillator is:

$\displaystyle \frac{\mathrm{d}^2x}{\mathrm{d}t^2} + 2\zeta\omega_0\frac{\mathrm{d}x}{\mathrm{d}t} + \omega^2_0 x = 0$

where $x$ is the position of the oscillator, $\omega_0$ is the frequency, and $\zeta$ is the damping ratio. To write this second-order ODE on standard form we introduce $p = \frac{\mathrm{d}x}{\mathrm{d}t}$:

$\displaystyle \frac{\mathrm{d}p}{\mathrm{d}t} = - 2\zeta\omega_0 p - \omega^2_0 x$

$\displaystyle \frac{\mathrm{d}x}{\mathrm{d}t} = p$

In the implementation of this example we will add extra arguments to the RHS function for the ODE, rather than using global variables as we did in the previous example. As a consequence of the extra arguments to the RHS, we need to pass an keyword argument args to the odeint function:

In [23]:

def dy(y, t, zeta, w0):
    """
    The right-hand side of the damped oscillator ODE
    """
    x, p = y #y[0], y[1]
    
    dx = p
    dp = -2 * zeta * w0 * p - w0**2 * x

    return array([dx, dp])

In [24]:

# initial state: 
y0 = [1.0, 0.0]

# time coodinate to solve the ODE for
t = linspace(0, 10, 1000)
w0 = 2*pi*1.0

In [25]:

# solve the ODE problem for three different values of the damping ratio

y1 = odeint(dy, y0, t, args=(0.0, w0)) # undamped
y2 = odeint(dy, y0, t, args=(0.2, w0)) # under damped
y3 = odeint(dy, y0, t, args=(1.0, w0)) # critial damping
y4 = odeint(dy, y0, t, args=(5.0, w0)) # over damped

In [26]:

plt.plot(t, y1[:,0], label="undamped", linewidth=0.5)
plt.plot(t, y2[:,0], label="under damped")
plt.plot(t, y3[:,0], label=r"critical damping")
plt.plot(t, y4[:,0], label="over damped")
plt.legend();

Fourier transform¶

Fourier transforms are one of the universal tools in computational physics, which appear over and over again in different contexts. SciPy provides functions for accessing the classic FFTPACK library from NetLib, which is an efficient and well tested FFT library written in FORTRAN. The SciPy API has a few additional convenience functions, but overall the API is closely related to the original FORTRAN library.

To use the fftpack module in a python program, include it using:

In [27]:

from scipy.fftpack import *
from numpy.fft import fftfreq

To demonstrate how to do a fast Fourier transform with SciPy, let's look at the FFT of the solution to the damped oscillator from the previous section:

In [28]:

# calculate the fast fourier transform
# y2 is the solution to the under-damped oscillator from the previous section
F = fft(y2[:,0]) 

# calculate the frequencies for the components in F
w = fftfreq(len(t), t[1]-t[0])

In [29]:

Out[29]:

array([  0.    ,   0.0999,   0.1998,   0.2997,   0.3996,   0.4995,
         0.5994,   0.6993,   0.7992,   0.8991,   0.999 ,   1.0989,
         1.1988,   1.2987,   1.3986,   1.4985,   1.5984,   1.6983,
         1.7982,   1.8981,   1.998 ,   2.0979,   2.1978,   2.2977,
         2.3976,   2.4975,   2.5974,   2.6973,   2.7972,   2.8971,
         2.997 ,   3.0969,   3.1968,   3.2967,   3.3966,   3.4965,
         3.5964,   3.6963,   3.7962,   3.8961,   3.996 ,   4.0959,
         4.1958,   4.2957,   4.3956,   4.4955,   4.5954,   4.6953,
         4.7952,   4.8951,   4.995 ,   5.0949,   5.1948,   5.2947,
         5.3946,   5.4945,   5.5944,   5.6943,   5.7942,   5.8941,
         5.994 ,   6.0939,   6.1938,   6.2937,   6.3936,   6.4935,
         6.5934,   6.6933,   6.7932,   6.8931,   6.993 ,   7.0929,
         7.1928,   7.2927,   7.3926,   7.4925,   7.5924,   7.6923,
         7.7922,   7.8921,   7.992 ,   8.0919,   8.1918,   8.2917,
         8.3916,   8.4915,   8.5914,   8.6913,   8.7912,   8.8911,
         8.991 ,   9.0909,   9.1908,   9.2907,   9.3906,   9.4905,
         9.5904,   9.6903,   9.7902,   9.8901,   9.99  ,  10.0899,
        10.1898,  10.2897,  10.3896,  10.4895,  10.5894,  10.6893,
        10.7892,  10.8891,  10.989 ,  11.0889,  11.1888,  11.2887,
        11.3886,  11.4885,  11.5884,  11.6883,  11.7882,  11.8881,
        11.988 ,  12.0879,  12.1878,  12.2877,  12.3876,  12.4875,
        12.5874,  12.6873,  12.7872,  12.8871,  12.987 ,  13.0869,
        13.1868,  13.2867,  13.3866,  13.4865,  13.5864,  13.6863,
        13.7862,  13.8861,  13.986 ,  14.0859,  14.1858,  14.2857,
        14.3856,  14.4855,  14.5854,  14.6853,  14.7852,  14.8851,
        14.985 ,  15.0849,  15.1848,  15.2847,  15.3846,  15.4845,
        15.5844,  15.6843,  15.7842,  15.8841,  15.984 ,  16.0839,
        16.1838,  16.2837,  16.3836,  16.4835,  16.5834,  16.6833,
        16.7832,  16.8831,  16.983 ,  17.0829,  17.1828,  17.2827,
        17.3826,  17.4825,  17.5824,  17.6823,  17.7822,  17.8821,
        17.982 ,  18.0819,  18.1818,  18.2817,  18.3816,  18.4815,
        18.5814,  18.6813,  18.7812,  18.8811,  18.981 ,  19.0809,
        19.1808,  19.2807,  19.3806,  19.4805,  19.5804,  19.6803,
        19.7802,  19.8801,  19.98  ,  20.0799,  20.1798,  20.2797,
        20.3796,  20.4795,  20.5794,  20.6793,  20.7792,  20.8791,
        20.979 ,  21.0789,  21.1788,  21.2787,  21.3786,  21.4785,
        21.5784,  21.6783,  21.7782,  21.8781,  21.978 ,  22.0779,
        22.1778,  22.2777,  22.3776,  22.4775,  22.5774,  22.6773,
        22.7772,  22.8771,  22.977 ,  23.0769,  23.1768,  23.2767,
        23.3766,  23.4765,  23.5764,  23.6763,  23.7762,  23.8761,
        23.976 ,  24.0759,  24.1758,  24.2757,  24.3756,  24.4755,
        24.5754,  24.6753,  24.7752,  24.8751,  24.975 ,  25.0749,
        25.1748,  25.2747,  25.3746,  25.4745,  25.5744,  25.6743,
        25.7742,  25.8741,  25.974 ,  26.0739,  26.1738,  26.2737,
        26.3736,  26.4735,  26.5734,  26.6733,  26.7732,  26.8731,
        26.973 ,  27.0729,  27.1728,  27.2727,  27.3726,  27.4725,
        27.5724,  27.6723,  27.7722,  27.8721,  27.972 ,  28.0719,
        28.1718,  28.2717,  28.3716,  28.4715,  28.5714,  28.6713,
        28.7712,  28.8711,  28.971 ,  29.0709,  29.1708,  29.2707,
        29.3706,  29.4705,  29.5704,  29.6703,  29.7702,  29.8701,
        29.97  ,  30.0699,  30.1698,  30.2697,  30.3696,  30.4695,
        30.5694,  30.6693,  30.7692,  30.8691,  30.969 ,  31.0689,
        31.1688,  31.2687,  31.3686,  31.4685,  31.5684,  31.6683,
        31.7682,  31.8681,  31.968 ,  32.0679,  32.1678,  32.2677,
        32.3676,  32.4675,  32.5674,  32.6673,  32.7672,  32.8671,
        32.967 ,  33.0669,  33.1668,  33.2667,  33.3666,  33.4665,
        33.5664,  33.6663,  33.7662,  33.8661,  33.966 ,  34.0659,
        34.1658,  34.2657,  34.3656,  34.4655,  34.5654,  34.6653,
        34.7652,  34.8651,  34.965 ,  35.0649,  35.1648,  35.2647,
        35.3646,  35.4645,  35.5644,  35.6643,  35.7642,  35.8641,
        35.964 ,  36.0639,  36.1638,  36.2637,  36.3636,  36.4635,
        36.5634,  36.6633,  36.7632,  36.8631,  36.963 ,  37.0629,
        37.1628,  37.2627,  37.3626,  37.4625,  37.5624,  37.6623,
        37.7622,  37.8621,  37.962 ,  38.0619,  38.1618,  38.2617,
        38.3616,  38.4615,  38.5614,  38.6613,  38.7612,  38.8611,
        38.961 ,  39.0609,  39.1608,  39.2607,  39.3606,  39.4605,
        39.5604,  39.6603,  39.7602,  39.8601,  39.96  ,  40.0599,
        40.1598,  40.2597,  40.3596,  40.4595,  40.5594,  40.6593,
        40.7592,  40.8591,  40.959 ,  41.0589,  41.1588,  41.2587,
        41.3586,  41.4585,  41.5584,  41.6583,  41.7582,  41.8581,
        41.958 ,  42.0579,  42.1578,  42.2577,  42.3576,  42.4575,
        42.5574,  42.6573,  42.7572,  42.8571,  42.957 ,  43.0569,
        43.1568,  43.2567,  43.3566,  43.4565,  43.5564,  43.6563,
        43.7562,  43.8561,  43.956 ,  44.0559,  44.1558,  44.2557,
        44.3556,  44.4555,  44.5554,  44.6553,  44.7552,  44.8551,
        44.955 ,  45.0549,  45.1548,  45.2547,  45.3546,  45.4545,
        45.5544,  45.6543,  45.7542,  45.8541,  45.954 ,  46.0539,
        46.1538,  46.2537,  46.3536,  46.4535,  46.5534,  46.6533,
        46.7532,  46.8531,  46.953 ,  47.0529,  47.1528,  47.2527,
        47.3526,  47.4525,  47.5524,  47.6523,  47.7522,  47.8521,
        47.952 ,  48.0519,  48.1518,  48.2517,  48.3516,  48.4515,
        48.5514,  48.6513,  48.7512,  48.8511,  48.951 ,  49.0509,
        49.1508,  49.2507,  49.3506,  49.4505,  49.5504,  49.6503,
        49.7502,  49.8501, -49.95  , -49.8501, -49.7502, -49.6503,
       -49.5504, -49.4505, -49.3506, -49.2507, -49.1508, -49.0509,
       -48.951 , -48.8511, -48.7512, -48.6513, -48.5514, -48.4515,
       -48.3516, -48.2517, -48.1518, -48.0519, -47.952 , -47.8521,
       -47.7522, -47.6523, -47.5524, -47.4525, -47.3526, -47.2527,
       -47.1528, -47.0529, -46.953 , -46.8531, -46.7532, -46.6533,
       -46.5534, -46.4535, -46.3536, -46.2537, -46.1538, -46.0539,
       -45.954 , -45.8541, -45.7542, -45.6543, -45.5544, -45.4545,
       -45.3546, -45.2547, -45.1548, -45.0549, -44.955 , -44.8551,
       -44.7552, -44.6553, -44.5554, -44.4555, -44.3556, -44.2557,
       -44.1558, -44.0559, -43.956 , -43.8561, -43.7562, -43.6563,
       -43.5564, -43.4565, -43.3566, -43.2567, -43.1568, -43.0569,
       -42.957 , -42.8571, -42.7572, -42.6573, -42.5574, -42.4575,
       -42.3576, -42.2577, -42.1578, -42.0579, -41.958 , -41.8581,
       -41.7582, -41.6583, -41.5584, -41.4585, -41.3586, -41.2587,
       -41.1588, -41.0589, -40.959 , -40.8591, -40.7592, -40.6593,
       -40.5594, -40.4595, -40.3596, -40.2597, -40.1598, -40.0599,
       -39.96  , -39.8601, -39.7602, -39.6603, -39.5604, -39.4605,
       -39.3606, -39.2607, -39.1608, -39.0609, -38.961 , -38.8611,
       -38.7612, -38.6613, -38.5614, -38.4615, -38.3616, -38.2617,
       -38.1618, -38.0619, -37.962 , -37.8621, -37.7622, -37.6623,
       -37.5624, -37.4625, -37.3626, -37.2627, -37.1628, -37.0629,
       -36.963 , -36.8631, -36.7632, -36.6633, -36.5634, -36.4635,
       -36.3636, -36.2637, -36.1638, -36.0639, -35.964 , -35.8641,
       -35.7642, -35.6643, -35.5644, -35.4645, -35.3646, -35.2647,
       -35.1648, -35.0649, -34.965 , -34.8651, -34.7652, -34.6653,
       -34.5654, -34.4655, -34.3656, -34.2657, -34.1658, -34.0659,
       -33.966 , -33.8661, -33.7662, -33.6663, -33.5664, -33.4665,
       -33.3666, -33.2667, -33.1668, -33.0669, -32.967 , -32.8671,
       -32.7672, -32.6673, -32.5674, -32.4675, -32.3676, -32.2677,
       -32.1678, -32.0679, -31.968 , -31.8681, -31.7682, -31.6683,
       -31.5684, -31.4685, -31.3686, -31.2687, -31.1688, -31.0689,
       -30.969 , -30.8691, -30.7692, -30.6693, -30.5694, -30.4695,
       -30.3696, -30.2697, -30.1698, -30.0699, -29.97  , -29.8701,
       -29.7702, -29.6703, -29.5704, -29.4705, -29.3706, -29.2707,
       -29.1708, -29.0709, -28.971 , -28.8711, -28.7712, -28.6713,
       -28.5714, -28.4715, -28.3716, -28.2717, -28.1718, -28.0719,
       -27.972 , -27.8721, -27.7722, -27.6723, -27.5724, -27.4725,
       -27.3726, -27.2727, -27.1728, -27.0729, -26.973 , -26.8731,
       -26.7732, -26.6733, -26.5734, -26.4735, -26.3736, -26.2737,
       -26.1738, -26.0739, -25.974 , -25.8741, -25.7742, -25.6743,
       -25.5744, -25.4745, -25.3746, -25.2747, -25.1748, -25.0749,
       -24.975 , -24.8751, -24.7752, -24.6753, -24.5754, -24.4755,
       -24.3756, -24.2757, -24.1758, -24.0759, -23.976 , -23.8761,
       -23.7762, -23.6763, -23.5764, -23.4765, -23.3766, -23.2767,
       -23.1768, -23.0769, -22.977 , -22.8771, -22.7772, -22.6773,
       -22.5774, -22.4775, -22.3776, -22.2777, -22.1778, -22.0779,
       -21.978 , -21.8781, -21.7782, -21.6783, -21.5784, -21.4785,
       -21.3786, -21.2787, -21.1788, -21.0789, -20.979 , -20.8791,
       -20.7792, -20.6793, -20.5794, -20.4795, -20.3796, -20.2797,
       -20.1798, -20.0799, -19.98  , -19.8801, -19.7802, -19.6803,
       -19.5804, -19.4805, -19.3806, -19.2807, -19.1808, -19.0809,
       -18.981 , -18.8811, -18.7812, -18.6813, -18.5814, -18.4815,
       -18.3816, -18.2817, -18.1818, -18.0819, -17.982 , -17.8821,
       -17.7822, -17.6823, -17.5824, -17.4825, -17.3826, -17.2827,
       -17.1828, -17.0829, -16.983 , -16.8831, -16.7832, -16.6833,
       -16.5834, -16.4835, -16.3836, -16.2837, -16.1838, -16.0839,
       -15.984 , -15.8841, -15.7842, -15.6843, -15.5844, -15.4845,
       -15.3846, -15.2847, -15.1848, -15.0849, -14.985 , -14.8851,
       -14.7852, -14.6853, -14.5854, -14.4855, -14.3856, -14.2857,
       -14.1858, -14.0859, -13.986 , -13.8861, -13.7862, -13.6863,
       -13.5864, -13.4865, -13.3866, -13.2867, -13.1868, -13.0869,
       -12.987 , -12.8871, -12.7872, -12.6873, -12.5874, -12.4875,
       -12.3876, -12.2877, -12.1878, -12.0879, -11.988 , -11.8881,
       -11.7882, -11.6883, -11.5884, -11.4885, -11.3886, -11.2887,
       -11.1888, -11.0889, -10.989 , -10.8891, -10.7892, -10.6893,
       -10.5894, -10.4895, -10.3896, -10.2897, -10.1898, -10.0899,
        -9.99  ,  -9.8901,  -9.7902,  -9.6903,  -9.5904,  -9.4905,
        -9.3906,  -9.2907,  -9.1908,  -9.0909,  -8.991 ,  -8.8911,
        -8.7912,  -8.6913,  -8.5914,  -8.4915,  -8.3916,  -8.2917,
        -8.1918,  -8.0919,  -7.992 ,  -7.8921,  -7.7922,  -7.6923,
        -7.5924,  -7.4925,  -7.3926,  -7.2927,  -7.1928,  -7.0929,
        -6.993 ,  -6.8931,  -6.7932,  -6.6933,  -6.5934,  -6.4935,
        -6.3936,  -6.2937,  -6.1938,  -6.0939,  -5.994 ,  -5.8941,
        -5.7942,  -5.6943,  -5.5944,  -5.4945,  -5.3946,  -5.2947,
        -5.1948,  -5.0949,  -4.995 ,  -4.8951,  -4.7952,  -4.6953,
        -4.5954,  -4.4955,  -4.3956,  -4.2957,  -4.1958,  -4.0959,
        -3.996 ,  -3.8961,  -3.7962,  -3.6963,  -3.5964,  -3.4965,
        -3.3966,  -3.2967,  -3.1968,  -3.0969,  -2.997 ,  -2.8971,
        -2.7972,  -2.6973,  -2.5974,  -2.4975,  -2.3976,  -2.2977,
        -2.1978,  -2.0979,  -1.998 ,  -1.8981,  -1.7982,  -1.6983,
        -1.5984,  -1.4985,  -1.3986,  -1.2987,  -1.1988,  -1.0989,
        -0.999 ,  -0.8991,  -0.7992,  -0.6993,  -0.5994,  -0.4995,
        -0.3996,  -0.2997,  -0.1998,  -0.0999])

In [30]:

plt.plot(w, abs(F));

Properties of Fourier transform of a real signal: \begin{eqnarray} && F(\omega) = \int e^{i\omega t} x(t) dt\\ && F^*(\omega) = F(-\omega)\\ && Re(F(\omega)) = Re(F(-\omega))\\ && Im(F(\omega)) = -Im(F(-\omega)) \end{eqnarray}

In [31]:

plt.plot(w,F.real)
plt.plot(w,F.imag)
plt.xlim(-2,2)

Out[31]:

(-2.0, 2.0)

Since the signal is real, the spectrum is symmetric. We therefore only need to plot the part that corresponds to the postive frequencies. To extract that part of the w and F we can use some of the indexing tricks for NumPy arrays that we saw in Lecture 2:

In [32]:

indices = where(w > 0) # select only indices for elements that corresponds to positive frequencies
w_pos = w[indices]
F_pos = F[indices]

In [33]:

plt.plot(w_pos,F_pos.real);
plt.xlim(0,3);

As expected, we now see a peak in the spectrum that is centered around 1, which is the frequency we used in the damped oscillator example.

Optimization¶

Optimization (finding minima or maxima of a function) is a large field in mathematics, and optimization of complicated functions or in many variables can be rather involved. Here we will only look at a few very simple cases. For a more detailed introduction to optimization with SciPy see: http://scipy-lectures.github.com/advanced/mathematical_optimization/index.html

To use the optimization module in scipy first include the optimize module:

In [34]:

from scipy import optimize

Finding a minima¶

Let's first look at how to find the minima of a simple function of a single variable:

In [35]:

def f(x):
    return 4*x**3 + (x-2)**2 + x**4

In [36]:

x = linspace(-5, 3, 100)
plt.plot(x, f(x));

We can use the fmin_bfgs function to find the minima of a function:

In [37]:

x_min = optimize.fmin_bfgs(f, -2)
x_min 

Optimization terminated successfully.
         Current function value: -3.506641
         Iterations: 5
         Function evaluations: 16
         Gradient evaluations: 8

Out[37]:

array([-2.67298155])

In [38]:

optimize.fmin_bfgs(f, 0.5) 

Optimization terminated successfully.
         Current function value: 2.804988
         Iterations: 3
         Function evaluations: 10
         Gradient evaluations: 5

Out[38]:

array([0.46961745])

We can also use the brent or fminbound functions. They have a bit different syntax and use different algorithms.

In [39]:

optimize.brent(f, brack=(1,2))

Out[39]:

0.4696174340948085

In [40]:

optimize.fminbound(f, -4, 2)

Out[40]:

-2.6729822917513886

Finding a solution to a function, i.e., zeros¶

To find the root for a function of the form $f(x) = 0$ we can use the fsolve function.

It is based on Powell's hybrid method as implemented in MINPACK’s library (hybrd): https://www.extremeoptimization.com/Documentation/Mathematics/Solving-Equations/Solving-Systems-of-Non-Linear-Equations.aspx

Powell's dogleg method, also called Powell's hybrid method, attempts to minimize the sum of the squares of the function values. It does this using a combination of Newton's* method and the steepest descent method. This is a so-called trust region method. This means that every step moves the current point to within a finite region. This makes the method more stable than Newton's method.*

On the other hand, the fact that the method is, in essence, a specialized minimizer means that the algorithm can get stuck in a local minimum that does not correspond to a solution of the system of equations.

In [41]:

Image('https://upload.wikimedia.org/wikipedia/commons/e/e0/NewtonIteration_Ani.gif')

Out[41]:

<IPython.core.display.Image object>

In [132]:

Image('https://upload.wikimedia.org/wikipedia/commons/thumb/f/ff/Gradient_descent.svg/1920px-Gradient_descent.svg.png')

Out[132]:

To use fsolve, we need an initial guess:

In [42]:

omega_c = 3.0
def f(omega):
    # a transcendental equation: resonance frequencies of a low-Q SQUID terminated microwave resonator
    return tan(2*pi*omega) - omega_c/omega

In [43]:

x = linspace(1e-6, 3, 1000)
y = f(x)
plt.plot(x,y)
plt.ylim(-5,5)

Out[43]:

(-5.0, 5.0)

In [44]:

x = linspace(1e-6, 3, 1000)
y = f(x)
mask = where(abs(y) > 50)
x[mask] = y[mask] = NaN # get rid of vertical line when the function flip sign

plt.plot(x, y)
plt.plot([0, 3], [0, 0], 'k')
plt.ylim(-5,5);

In [45]:

optimize.fsolve(f, 0.1)

Out[45]:

array([0.23743014])

In [46]:

optimize.fsolve(f, 0.6)

Out[46]:

array([0.71286972])

In [47]:

optimize.fsolve(f, 1.1)

Out[47]:

array([1.18990285])

Change of sign can occur when there is a zero or a pole. To use bracketing technique, we need to carefully bracket the zeros (and not the poles)

In [48]:

f(0.01),f(0.5)

Out[48]:

(-299.93708533274634, -6.0)

In [49]:

f(0.001),f(0.25),f(0.3),f(0.73)

Out[49]:

(-2999.993716732008,
 1.6331239353195358e+16,
 -13.077683537175254,
 3.806226047209912)

toms748(f, a, b, args=(), k=1, xtol=2e-12, rtol=8.881784197001252e-16, maxiter=100, full_output=False, disp=True)

Finds a zero of the function f on the interval [a , b], where f(a) and f(b) must have opposite signs.

It uses a mixture of inverse cubic interpolation and "Newton-quadratic" steps. [APS1995].

In [50]:

[optimize.toms748(f,0.01,0.25),optimize.toms748(f, 0.3, 0.73)]

Out[50]:

[0.2374301406360339, 0.7128697158579146]

In [51]:

optimize.toms748(f,-0.012,0.01)
optimize.fsolve(f, 1e-16)

Out[51]:

array([0.23743014])

Interpolation¶

Interpolation is simple and convenient in scipy: The interp1d function, when given arrays describing X and Y data, returns and object that behaves like a function that can be called for an arbitrary value of x (in the range covered by X), and it returns the corresponding interpolated y value:

In [52]:

from scipy.interpolate import *
from numpy import random

In [53]:

def f(x):
    return sin(x)

In [54]:

n = linspace(0, 9, 10)  
x = linspace(0, 9, 300)

y_meas = f(n) + 0.05 * random.randn(len(n)) # simulate measurement with noise
y_real = f(x)

linear_interpolation = interp1d(n, y_meas, kind='linear')
y_interp1 = linear_interpolation(x)

cubic_interpolation = interp1d(n, y_meas, kind='cubic')
y_interp2 = cubic_interpolation(x)

cubic_smooth = UnivariateSpline(n, y_meas, s=0.5)
y_interp3 = cubic_smooth(x)

In [55]:

fig, ax = plt.subplots(figsize=(10,4))
ax.plot(n, y_meas, 'bs', label='noisy data')
ax.plot(x, y_real, 'k', lw=2, label='true function')
ax.plot(x, y_interp1, 'r', label='linear interp')
ax.plot(x, y_interp2, 'g', label='cubic interp')
ax.plot(x, y_interp3, 'b', label='cubic smooth')
ax.legend(loc=3);

In [56]:

cubic_interpolation = interp1d(n, f(n), kind='cubic')
diff = cubic_interpolation(x)-f(x)

cubic_smooth = UnivariateSpline(n, f(n), s=0.0, k=4)
diff2 = cubic_smooth(x)-f(x)

plt.plot(x, diff)
plt.plot(x, diff2)

Out[56]:

[<matplotlib.lines.Line2D at 0x7f7da8e657f0>]

interpolate.make_interp_spline allows boundary condition at the two ends.

bc_type2 -- tuple or None

Default is None, which means choosing the boundary conditions automatically. Otherwise, it must be a length-two tuple where the first element sets the boundary conditions at x[0] and the second element sets the boundary conditions at x[-1]. Each of these must be an iterable of pairs (order, value) which gives the values of derivatives of specified orders at the given edge of the interpolation interval. Alternatively, the following string aliases are recognized:`

"clamped": The first derivatives at the ends are zero. This is equivalent to bc_type=([(1, 0.0)], [(1, 0.0)]).
"natural": The second derivatives at ends are zero. This is equivalent to bc_type=([(2, 0.0)], [(2, 0.0)]).
"not-a-knot" (default): The first and second segments are the same polynomial. This is equivalent to having bc_type=None. In this case the spline requires that the third derivative of the spline is continuous at x[0] and x[-1]

Also interp2d exists with similar syntax

UnivariateSpline is object oriented analog for spline interpolation, and offers a bit more functionality: s=Positive smoothing factor

In [82]:

from scipy import interpolate
fx = interpolate.UnivariateSpline(x, sin(x),s=0)
fx1=fx.derivative(1)
fx2=fx.derivative(2)

In [83]:

plt.plot(x,fx(x))
plt.plot(x,sin(x))
plt.plot(x,fx1(x))

Out[83]:

[<matplotlib.lines.Line2D at 0x7f7db824cbb0>]

In [84]:

plt.plot(x,fx1(x))
plt.plot(x,fx2(x))

Out[84]:

[<matplotlib.lines.Line2D at 0x7f7d68169130>]