In [1]:
%matplotlib inline

In [2]:
import numpy as np
import matplotlib.pyplot as plt

SciPy

SciPy is a collection of numerical algorithms with python interfaces. In many cases, these interfaces are wrappers around standard numerical libraries that have been developed in the community and are used with other languages. Usually detailed references are available to explain the implementation.

There are many subpackages. Generally, you load the subpackages separately, e.g.

from scipy import linalg, optimize

then you have access to the methods in those namespaces

Numerical Methods

One thing to keep in mind -- all numerical methods have strengths and weaknesses, and make assumptions. You should always do some research into the method to understand what it is doing.

It is also always a good idea to run a new method on some test where you know the answer, to make sure it is behaving as expected.

Integration

we'll do some integrals of the form

$$I = \int_a^b f(x) dx$$

We can imagine two situations:

  • our function $f(x)$ is given by an analytic expression. This gives us the freedom to pick our integration points, and in general can allow us to optimize our result and get high accuracy
  • our function $f(x)$ is defined on at a set of (possibly regular spaced) points.

In numerical analysis, the term quadrature is used to describe any integration method that represents the integral as the weighted sum of a discrete number of points.


In [3]:
from scipy import integrate
help(integrate)


Help on package scipy.integrate in scipy:

NAME
    scipy.integrate

DESCRIPTION
    =============================================
    Integration and ODEs (:mod:`scipy.integrate`)
    =============================================
    
    .. currentmodule:: scipy.integrate
    
    Integrating functions, given function object
    ============================================
    
    .. autosummary::
       :toctree: generated/
    
       quad          -- General purpose integration
       dblquad       -- General purpose double integration
       tplquad       -- General purpose triple integration
       nquad         -- General purpose n-dimensional integration
       fixed_quad    -- Integrate func(x) using Gaussian quadrature of order n
       quadrature    -- Integrate with given tolerance using Gaussian quadrature
       romberg       -- Integrate func using Romberg integration
       quad_explain  -- Print information for use of quad
       newton_cotes  -- Weights and error coefficient for Newton-Cotes integration
       IntegrationWarning -- Warning on issues during integration
    
    Integrating functions, given fixed samples
    ==========================================
    
    .. autosummary::
       :toctree: generated/
    
       trapz         -- Use trapezoidal rule to compute integral.
       cumtrapz      -- Use trapezoidal rule to cumulatively compute integral.
       simps         -- Use Simpson's rule to compute integral from samples.
       romb          -- Use Romberg Integration to compute integral from
                     -- (2**k + 1) evenly-spaced samples.
    
    .. seealso::
    
       :mod:`scipy.special` for orthogonal polynomials (special) for Gaussian
       quadrature roots and weights for other weighting factors and regions.
    
    Integrators of ODE systems
    ==========================
    
    .. autosummary::
       :toctree: generated/
    
       odeint        -- General integration of ordinary differential equations.
       ode           -- Integrate ODE using VODE and ZVODE routines.
       complex_ode   -- Convert a complex-valued ODE to real-valued and integrate.
       solve_bvp     -- Solve a boundary value problem for a system of ODEs.

PACKAGE CONTENTS
    _bvp
    _dop
    _ode
    _odepack
    _quadpack
    _test_multivariate
    _test_odeint_banded
    lsoda
    odepack
    quadpack
    quadrature
    setup
    vode

CLASSES
    builtins.UserWarning(builtins.Warning)
        scipy.integrate.quadpack.IntegrationWarning
    builtins.object
        scipy.integrate._ode.ode
            scipy.integrate._ode.complex_ode
    
    class IntegrationWarning(builtins.UserWarning)
     |  Warning on issues during integration.
     |  
     |  Method resolution order:
     |      IntegrationWarning
     |      builtins.UserWarning
     |      builtins.Warning
     |      builtins.Exception
     |      builtins.BaseException
     |      builtins.object
     |  
     |  Data descriptors defined here:
     |  
     |  __weakref__
     |      list of weak references to the object (if defined)
     |  
     |  ----------------------------------------------------------------------
     |  Methods inherited from builtins.UserWarning:
     |  
     |  __init__(self, /, *args, **kwargs)
     |      Initialize self.  See help(type(self)) for accurate signature.
     |  
     |  __new__(*args, **kwargs) from builtins.type
     |      Create and return a new object.  See help(type) for accurate signature.
     |  
     |  ----------------------------------------------------------------------
     |  Methods inherited from builtins.BaseException:
     |  
     |  __delattr__(self, name, /)
     |      Implement delattr(self, name).
     |  
     |  __getattribute__(self, name, /)
     |      Return getattr(self, name).
     |  
     |  __reduce__(...)
     |      helper for pickle
     |  
     |  __repr__(self, /)
     |      Return repr(self).
     |  
     |  __setattr__(self, name, value, /)
     |      Implement setattr(self, name, value).
     |  
     |  __setstate__(...)
     |  
     |  __str__(self, /)
     |      Return str(self).
     |  
     |  with_traceback(...)
     |      Exception.with_traceback(tb) --
     |      set self.__traceback__ to tb and return self.
     |  
     |  ----------------------------------------------------------------------
     |  Data descriptors inherited from builtins.BaseException:
     |  
     |  __cause__
     |      exception cause
     |  
     |  __context__
     |      exception context
     |  
     |  __dict__
     |  
     |  __suppress_context__
     |  
     |  __traceback__
     |  
     |  args
    
    class complex_ode(ode)
     |  A wrapper of ode for complex systems.
     |  
     |  This functions similarly as `ode`, but re-maps a complex-valued
     |  equation system to a real-valued one before using the integrators.
     |  
     |  Parameters
     |  ----------
     |  f : callable ``f(t, y, *f_args)``
     |      Rhs of the equation. t is a scalar, ``y.shape == (n,)``.
     |      ``f_args`` is set by calling ``set_f_params(*args)``.
     |  jac : callable ``jac(t, y, *jac_args)``
     |      Jacobian of the rhs, ``jac[i,j] = d f[i] / d y[j]``.
     |      ``jac_args`` is set by calling ``set_f_params(*args)``.
     |  
     |  Attributes
     |  ----------
     |  t : float
     |      Current time.
     |  y : ndarray
     |      Current variable values.
     |  
     |  Examples
     |  --------
     |  For usage examples, see `ode`.
     |  
     |  Method resolution order:
     |      complex_ode
     |      ode
     |      builtins.object
     |  
     |  Methods defined here:
     |  
     |  __init__(self, f, jac=None)
     |      Initialize self.  See help(type(self)) for accurate signature.
     |  
     |  integrate(self, t, step=0, relax=0)
     |      Find y=y(t), set y as an initial condition, and return y.
     |  
     |  set_initial_value(self, y, t=0.0)
     |      Set initial conditions y(t) = y.
     |  
     |  set_integrator(self, name, **integrator_params)
     |      Set integrator by name.
     |      
     |      Parameters
     |      ----------
     |      name : str
     |          Name of the integrator
     |      integrator_params
     |          Additional parameters for the integrator.
     |  
     |  set_solout(self, solout)
     |      Set callable to be called at every successful integration step.
     |      
     |      Parameters
     |      ----------
     |      solout : callable
     |          ``solout(t, y)`` is called at each internal integrator step,
     |          t is a scalar providing the current independent position
     |          y is the current soloution ``y.shape == (n,)``
     |          solout should return -1 to stop integration
     |          otherwise it should return None or 0
     |  
     |  ----------------------------------------------------------------------
     |  Data descriptors defined here:
     |  
     |  y
     |  
     |  ----------------------------------------------------------------------
     |  Methods inherited from ode:
     |  
     |  set_f_params(self, *args)
     |      Set extra parameters for user-supplied function f.
     |  
     |  set_jac_params(self, *args)
     |      Set extra parameters for user-supplied function jac.
     |  
     |  successful(self)
     |      Check if integration was successful.
     |  
     |  ----------------------------------------------------------------------
     |  Data descriptors inherited from ode:
     |  
     |  __dict__
     |      dictionary for instance variables (if defined)
     |  
     |  __weakref__
     |      list of weak references to the object (if defined)
    
    class ode(builtins.object)
     |  A generic interface class to numeric integrators.
     |  
     |  Solve an equation system :math:`y'(t) = f(t,y)` with (optional) ``jac = df/dy``.
     |  
     |  *Note*: The first two arguments of ``f(t, y, ...)`` are in the
     |  opposite order of the arguments in the system definition function used
     |  by `scipy.integrate.odeint`.
     |  
     |  Parameters
     |  ----------
     |  f : callable ``f(t, y, *f_args)``
     |      Right-hand side of the differential equation. t is a scalar,
     |      ``y.shape == (n,)``.
     |      ``f_args`` is set by calling ``set_f_params(*args)``.
     |      `f` should return a scalar, array or list (not a tuple).
     |  jac : callable ``jac(t, y, *jac_args)``, optional
     |      Jacobian of the right-hand side, ``jac[i,j] = d f[i] / d y[j]``.
     |      ``jac_args`` is set by calling ``set_jac_params(*args)``.
     |  
     |  Attributes
     |  ----------
     |  t : float
     |      Current time.
     |  y : ndarray
     |      Current variable values.
     |  
     |  See also
     |  --------
     |  odeint : an integrator with a simpler interface based on lsoda from ODEPACK
     |  quad : for finding the area under a curve
     |  
     |  Notes
     |  -----
     |  Available integrators are listed below. They can be selected using
     |  the `set_integrator` method.
     |  
     |  "vode"
     |  
     |      Real-valued Variable-coefficient Ordinary Differential Equation
     |      solver, with fixed-leading-coefficient implementation. It provides
     |      implicit Adams method (for non-stiff problems) and a method based on
     |      backward differentiation formulas (BDF) (for stiff problems).
     |  
     |      Source: http://www.netlib.org/ode/vode.f
     |  
     |      .. warning::
     |  
     |         This integrator is not re-entrant. You cannot have two `ode`
     |         instances using the "vode" integrator at the same time.
     |  
     |      This integrator accepts the following parameters in `set_integrator`
     |      method of the `ode` class:
     |  
     |      - atol : float or sequence
     |        absolute tolerance for solution
     |      - rtol : float or sequence
     |        relative tolerance for solution
     |      - lband : None or int
     |      - uband : None or int
     |        Jacobian band width, jac[i,j] != 0 for i-lband <= j <= i+uband.
     |        Setting these requires your jac routine to return the jacobian
     |        in packed format, jac_packed[i-j+uband, j] = jac[i,j]. The
     |        dimension of the matrix must be (lband+uband+1, len(y)).
     |      - method: 'adams' or 'bdf'
     |        Which solver to use, Adams (non-stiff) or BDF (stiff)
     |      - with_jacobian : bool
     |        This option is only considered when the user has not supplied a
     |        Jacobian function and has not indicated (by setting either band)
     |        that the Jacobian is banded.  In this case, `with_jacobian` specifies
     |        whether the iteration method of the ODE solver's correction step is
     |        chord iteration with an internally generated full Jacobian or
     |        functional iteration with no Jacobian.
     |      - nsteps : int
     |        Maximum number of (internally defined) steps allowed during one
     |        call to the solver.
     |      - first_step : float
     |      - min_step : float
     |      - max_step : float
     |        Limits for the step sizes used by the integrator.
     |      - order : int
     |        Maximum order used by the integrator,
     |        order <= 12 for Adams, <= 5 for BDF.
     |  
     |  "zvode"
     |  
     |      Complex-valued Variable-coefficient Ordinary Differential Equation
     |      solver, with fixed-leading-coefficient implementation.  It provides
     |      implicit Adams method (for non-stiff problems) and a method based on
     |      backward differentiation formulas (BDF) (for stiff problems).
     |  
     |      Source: http://www.netlib.org/ode/zvode.f
     |  
     |      .. warning::
     |  
     |         This integrator is not re-entrant. You cannot have two `ode`
     |         instances using the "zvode" integrator at the same time.
     |  
     |      This integrator accepts the same parameters in `set_integrator`
     |      as the "vode" solver.
     |  
     |      .. note::
     |  
     |          When using ZVODE for a stiff system, it should only be used for
     |          the case in which the function f is analytic, that is, when each f(i)
     |          is an analytic function of each y(j).  Analyticity means that the
     |          partial derivative df(i)/dy(j) is a unique complex number, and this
     |          fact is critical in the way ZVODE solves the dense or banded linear
     |          systems that arise in the stiff case.  For a complex stiff ODE system
     |          in which f is not analytic, ZVODE is likely to have convergence
     |          failures, and for this problem one should instead use DVODE on the
     |          equivalent real system (in the real and imaginary parts of y).
     |  
     |  "lsoda"
     |  
     |      Real-valued Variable-coefficient Ordinary Differential Equation
     |      solver, with fixed-leading-coefficient implementation. It provides
     |      automatic method switching between implicit Adams method (for non-stiff
     |      problems) and a method based on backward differentiation formulas (BDF)
     |      (for stiff problems).
     |  
     |      Source: http://www.netlib.org/odepack
     |  
     |      .. warning::
     |  
     |         This integrator is not re-entrant. You cannot have two `ode`
     |         instances using the "lsoda" integrator at the same time.
     |  
     |      This integrator accepts the following parameters in `set_integrator`
     |      method of the `ode` class:
     |  
     |      - atol : float or sequence
     |        absolute tolerance for solution
     |      - rtol : float or sequence
     |        relative tolerance for solution
     |      - lband : None or int
     |      - uband : None or int
     |        Jacobian band width, jac[i,j] != 0 for i-lband <= j <= i+uband.
     |        Setting these requires your jac routine to return the jacobian
     |        in packed format, jac_packed[i-j+uband, j] = jac[i,j].
     |      - with_jacobian : bool
     |        *Not used.*
     |      - nsteps : int
     |        Maximum number of (internally defined) steps allowed during one
     |        call to the solver.
     |      - first_step : float
     |      - min_step : float
     |      - max_step : float
     |        Limits for the step sizes used by the integrator.
     |      - max_order_ns : int
     |        Maximum order used in the nonstiff case (default 12).
     |      - max_order_s : int
     |        Maximum order used in the stiff case (default 5).
     |      - max_hnil : int
     |        Maximum number of messages reporting too small step size (t + h = t)
     |        (default 0)
     |      - ixpr : int
     |        Whether to generate extra printing at method switches (default False).
     |  
     |  "dopri5"
     |  
     |      This is an explicit runge-kutta method of order (4)5 due to Dormand &
     |      Prince (with stepsize control and dense output).
     |  
     |      Authors:
     |  
     |          E. Hairer and G. Wanner
     |          Universite de Geneve, Dept. de Mathematiques
     |          CH-1211 Geneve 24, Switzerland
     |          e-mail:  ernst.hairer@math.unige.ch, gerhard.wanner@math.unige.ch
     |  
     |      This code is described in [HNW93]_.
     |  
     |      This integrator accepts the following parameters in set_integrator()
     |      method of the ode class:
     |  
     |      - atol : float or sequence
     |        absolute tolerance for solution
     |      - rtol : float or sequence
     |        relative tolerance for solution
     |      - nsteps : int
     |        Maximum number of (internally defined) steps allowed during one
     |        call to the solver.
     |      - first_step : float
     |      - max_step : float
     |      - safety : float
     |        Safety factor on new step selection (default 0.9)
     |      - ifactor : float
     |      - dfactor : float
     |        Maximum factor to increase/decrease step size by in one step
     |      - beta : float
     |        Beta parameter for stabilised step size control.
     |      - verbosity : int
     |        Switch for printing messages (< 0 for no messages).
     |  
     |  "dop853"
     |  
     |      This is an explicit runge-kutta method of order 8(5,3) due to Dormand
     |      & Prince (with stepsize control and dense output).
     |  
     |      Options and references the same as "dopri5".
     |  
     |  Examples
     |  --------
     |  
     |  A problem to integrate and the corresponding jacobian:
     |  
     |  >>> from scipy.integrate import ode
     |  >>>
     |  >>> y0, t0 = [1.0j, 2.0], 0
     |  >>>
     |  >>> def f(t, y, arg1):
     |  ...     return [1j*arg1*y[0] + y[1], -arg1*y[1]**2]
     |  >>> def jac(t, y, arg1):
     |  ...     return [[1j*arg1, 1], [0, -arg1*2*y[1]]]
     |  
     |  The integration:
     |  
     |  >>> r = ode(f, jac).set_integrator('zvode', method='bdf')
     |  >>> r.set_initial_value(y0, t0).set_f_params(2.0).set_jac_params(2.0)
     |  >>> t1 = 10
     |  >>> dt = 1
     |  >>> while r.successful() and r.t < t1:
     |  ...     print(r.t+dt, r.integrate(r.t+dt))
     |  (1, array([-0.71038232+0.23749653j,  0.40000271+0.j        ]))
     |  (2.0, array([ 0.19098503-0.52359246j,  0.22222356+0.j        ]))
     |  (3.0, array([ 0.47153208+0.52701229j,  0.15384681+0.j        ]))
     |  (4.0, array([-0.61905937+0.30726255j,  0.11764744+0.j        ]))
     |  (5.0, array([ 0.02340997-0.61418799j,  0.09523835+0.j        ]))
     |  (6.0, array([ 0.58643071+0.339819j,  0.08000018+0.j      ]))
     |  (7.0, array([-0.52070105+0.44525141j,  0.06896565+0.j        ]))
     |  (8.0, array([-0.15986733-0.61234476j,  0.06060616+0.j        ]))
     |  (9.0, array([ 0.64850462+0.15048982j,  0.05405414+0.j        ]))
     |  (10.0, array([-0.38404699+0.56382299j,  0.04878055+0.j        ]))
     |  
     |  References
     |  ----------
     |  .. [HNW93] E. Hairer, S.P. Norsett and G. Wanner, Solving Ordinary
     |      Differential Equations i. Nonstiff Problems. 2nd edition.
     |      Springer Series in Computational Mathematics,
     |      Springer-Verlag (1993)
     |  
     |  Methods defined here:
     |  
     |  __init__(self, f, jac=None)
     |      Initialize self.  See help(type(self)) for accurate signature.
     |  
     |  integrate(self, t, step=0, relax=0)
     |      Find y=y(t), set y as an initial condition, and return y.
     |  
     |  set_f_params(self, *args)
     |      Set extra parameters for user-supplied function f.
     |  
     |  set_initial_value(self, y, t=0.0)
     |      Set initial conditions y(t) = y.
     |  
     |  set_integrator(self, name, **integrator_params)
     |      Set integrator by name.
     |      
     |      Parameters
     |      ----------
     |      name : str
     |          Name of the integrator.
     |      integrator_params
     |          Additional parameters for the integrator.
     |  
     |  set_jac_params(self, *args)
     |      Set extra parameters for user-supplied function jac.
     |  
     |  set_solout(self, solout)
     |      Set callable to be called at every successful integration step.
     |      
     |      Parameters
     |      ----------
     |      solout : callable
     |          ``solout(t, y)`` is called at each internal integrator step,
     |          t is a scalar providing the current independent position
     |          y is the current soloution ``y.shape == (n,)``
     |          solout should return -1 to stop integration
     |          otherwise it should return None or 0
     |  
     |  successful(self)
     |      Check if integration was successful.
     |  
     |  ----------------------------------------------------------------------
     |  Data descriptors defined here:
     |  
     |  __dict__
     |      dictionary for instance variables (if defined)
     |  
     |  __weakref__
     |      list of weak references to the object (if defined)
     |  
     |  y

FUNCTIONS
    cumtrapz(y, x=None, dx=1.0, axis=-1, initial=None)
        Cumulatively integrate y(x) using the composite trapezoidal rule.
        
        Parameters
        ----------
        y : array_like
            Values to integrate.
        x : array_like, optional
            The coordinate to integrate along.  If None (default), use spacing `dx`
            between consecutive elements in `y`.
        dx : int, optional
            Spacing between elements of `y`.  Only used if `x` is None.
        axis : int, optional
            Specifies the axis to cumulate.  Default is -1 (last axis).
        initial : scalar, optional
            If given, uses this value as the first value in the returned result.
            Typically this value should be 0.  Default is None, which means no
            value at ``x[0]`` is returned and `res` has one element less than `y`
            along the axis of integration.
        
        Returns
        -------
        res : ndarray
            The result of cumulative integration of `y` along `axis`.
            If `initial` is None, the shape is such that the axis of integration
            has one less value than `y`.  If `initial` is given, the shape is equal
            to that of `y`.
        
        See Also
        --------
        numpy.cumsum, numpy.cumprod
        quad: adaptive quadrature using QUADPACK
        romberg: adaptive Romberg quadrature
        quadrature: adaptive Gaussian quadrature
        fixed_quad: fixed-order Gaussian quadrature
        dblquad: double integrals
        tplquad: triple integrals
        romb: integrators for sampled data
        ode: ODE integrators
        odeint: ODE integrators
        
        Examples
        --------
        >>> from scipy import integrate
        >>> import matplotlib.pyplot as plt
        
        >>> x = np.linspace(-2, 2, num=20)
        >>> y = x
        >>> y_int = integrate.cumtrapz(y, x, initial=0)
        >>> plt.plot(x, y_int, 'ro', x, y[0] + 0.5 * x**2, 'b-')
        >>> plt.show()
    
    dblquad(func, a, b, gfun, hfun, args=(), epsabs=1.49e-08, epsrel=1.49e-08)
        Compute a double integral.
        
        Return the double (definite) integral of ``func(y, x)`` from ``x = a..b``
        and ``y = gfun(x)..hfun(x)``.
        
        Parameters
        ----------
        func : callable
            A Python function or method of at least two variables: y must be the
            first argument and x the second argument.
        a, b : float
            The limits of integration in x: `a` < `b`
        gfun : callable
            The lower boundary curve in y which is a function taking a single
            floating point argument (x) and returning a floating point result: a
            lambda function can be useful here.
        hfun : callable
            The upper boundary curve in y (same requirements as `gfun`).
        args : sequence, optional
            Extra arguments to pass to `func`.
        epsabs : float, optional
            Absolute tolerance passed directly to the inner 1-D quadrature
            integration. Default is 1.49e-8.
        epsrel : float, optional
            Relative tolerance of the inner 1-D integrals. Default is 1.49e-8.
        
        Returns
        -------
        y : float
            The resultant integral.
        abserr : float
            An estimate of the error.
        
        See also
        --------
        quad : single integral
        tplquad : triple integral
        nquad : N-dimensional integrals
        fixed_quad : fixed-order Gaussian quadrature
        quadrature : adaptive Gaussian quadrature
        odeint : ODE integrator
        ode : ODE integrator
        simps : integrator for sampled data
        romb : integrator for sampled data
        scipy.special : for coefficients and roots of orthogonal polynomials
    
    fixed_quad(func, a, b, args=(), n=5)
        Compute a definite integral using fixed-order Gaussian quadrature.
        
        Integrate `func` from `a` to `b` using Gaussian quadrature of
        order `n`.
        
        Parameters
        ----------
        func : callable
            A Python function or method to integrate (must accept vector inputs).
        a : float
            Lower limit of integration.
        b : float
            Upper limit of integration.
        args : tuple, optional
            Extra arguments to pass to function, if any.
        n : int, optional
            Order of quadrature integration. Default is 5.
        
        Returns
        -------
        val : float
            Gaussian quadrature approximation to the integral
        none : None
            Statically returned value of None
        
        
        See Also
        --------
        quad : adaptive quadrature using QUADPACK
        dblquad : double integrals
        tplquad : triple integrals
        romberg : adaptive Romberg quadrature
        quadrature : adaptive Gaussian quadrature
        romb : integrators for sampled data
        simps : integrators for sampled data
        cumtrapz : cumulative integration for sampled data
        ode : ODE integrator
        odeint : ODE integrator
    
    newton_cotes(rn, equal=0)
        Return weights and error coefficient for Newton-Cotes integration.
        
        Suppose we have (N+1) samples of f at the positions
        x_0, x_1, ..., x_N.  Then an N-point Newton-Cotes formula for the
        integral between x_0 and x_N is:
        
        :math:`\int_{x_0}^{x_N} f(x)dx = \Delta x \sum_{i=0}^{N} a_i f(x_i)
        + B_N (\Delta x)^{N+2} f^{N+1} (\xi)`
        
        where :math:`\xi \in [x_0,x_N]`
        and :math:`\Delta x = \frac{x_N-x_0}{N}` is the average samples spacing.
        
        If the samples are equally-spaced and N is even, then the error
        term is :math:`B_N (\Delta x)^{N+3} f^{N+2}(\xi)`.
        
        Parameters
        ----------
        rn : int
            The integer order for equally-spaced data or the relative positions of
            the samples with the first sample at 0 and the last at N, where N+1 is
            the length of `rn`.  N is the order of the Newton-Cotes integration.
        equal : int, optional
            Set to 1 to enforce equally spaced data.
        
        Returns
        -------
        an : ndarray
            1-D array of weights to apply to the function at the provided sample
            positions.
        B : float
            Error coefficient.
        
        Notes
        -----
        Normally, the Newton-Cotes rules are used on smaller integration
        regions and a composite rule is used to return the total integral.
    
    nquad(func, ranges, args=None, opts=None, full_output=False)
        Integration over multiple variables.
        
        Wraps `quad` to enable integration over multiple variables.
        Various options allow improved integration of discontinuous functions, as
        well as the use of weighted integration, and generally finer control of the
        integration process.
        
        Parameters
        ----------
        func : callable
            The function to be integrated. Has arguments of ``x0, ... xn``,
            ``t0, tm``, where integration is carried out over ``x0, ... xn``, which
            must be floats.  Function signature should be
            ``func(x0, x1, ..., xn, t0, t1, ..., tm)``.  Integration is carried out
            in order.  That is, integration over ``x0`` is the innermost integral,
            and ``xn`` is the outermost.
            If performance is a concern, this function may be a ctypes function of
            the form::
        
                f(int n, double args[n])
        
            where ``n`` is the number of extra parameters and args is an array
            of doubles of the additional parameters.  This function may then
            be compiled to a dynamic/shared library then imported through
            ``ctypes``, setting the function's argtypes to ``(c_int, c_double)``,
            and the function's restype to ``(c_double)``.  Its pointer may then be
            passed into `nquad` normally.
            This allows the underlying Fortran library to evaluate the function in
            the innermost integration calls without callbacks to Python, and also
            speeds up the evaluation of the function itself.
        ranges : iterable object
            Each element of ranges may be either a sequence  of 2 numbers, or else
            a callable that returns such a sequence.  ``ranges[0]`` corresponds to
            integration over x0, and so on.  If an element of ranges is a callable,
            then it will be called with all of the integration arguments available,
            as well as any parametric arguments. e.g. if 
            ``func = f(x0, x1, x2, t0, t1)``, then ``ranges[0]`` may be defined as
            either ``(a, b)`` or else as ``(a, b) = range0(x1, x2, t0, t1)``.
        args : iterable object, optional
            Additional arguments ``t0, ..., tn``, required by `func`, `ranges`, and
            ``opts``.
        opts : iterable object or dict, optional
            Options to be passed to `quad`.  May be empty, a dict, or
            a sequence of dicts or functions that return a dict.  If empty, the
            default options from scipy.integrate.quad are used.  If a dict, the same
            options are used for all levels of integraion.  If a sequence, then each
            element of the sequence corresponds to a particular integration. e.g.
            opts[0] corresponds to integration over x0, and so on. If a callable, 
            the signature must be the same as for ``ranges``. The available
            options together with their default values are:
        
              - epsabs = 1.49e-08
              - epsrel = 1.49e-08
              - limit  = 50
              - points = None
              - weight = None
              - wvar   = None
              - wopts  = None
        
            For more information on these options, see `quad` and `quad_explain`.
        
        full_output : bool, optional
            Partial implementation of ``full_output`` from scipy.integrate.quad. 
            The number of integrand function evaluations ``neval`` can be obtained 
            by setting ``full_output=True`` when calling nquad.
        
        Returns
        -------
        result : float
            The result of the integration.
        abserr : float
            The maximum of the estimates of the absolute error in the various
            integration results.
        out_dict : dict, optional
            A dict containing additional information on the integration. 
        
        See Also
        --------
        quad : 1-dimensional numerical integration
        dblquad, tplquad : double and triple integrals
        fixed_quad : fixed-order Gaussian quadrature
        quadrature : adaptive Gaussian quadrature
        
        Examples
        --------
        >>> from scipy import integrate
        >>> func = lambda x0,x1,x2,x3 : x0**2 + x1*x2 - x3**3 + np.sin(x0) + (
        ...                                 1 if (x0-.2*x3-.5-.25*x1>0) else 0)
        >>> points = [[lambda x1,x2,x3 : 0.2*x3 + 0.5 + 0.25*x1], [], [], []]
        >>> def opts0(*args, **kwargs):
        ...     return {'points':[0.2*args[2] + 0.5 + 0.25*args[0]]}
        >>> integrate.nquad(func, [[0,1], [-1,1], [.13,.8], [-.15,1]],
        ...                 opts=[opts0,{},{},{}], full_output=True)
        (1.5267454070738633, 2.9437360001402324e-14, {'neval': 388962})
        
        >>> scale = .1
        >>> def func2(x0, x1, x2, x3, t0, t1):
        ...     return x0*x1*x3**2 + np.sin(x2) + 1 + (1 if x0+t1*x1-t0>0 else 0)
        >>> def lim0(x1, x2, x3, t0, t1):
        ...     return [scale * (x1**2 + x2 + np.cos(x3)*t0*t1 + 1) - 1,
        ...             scale * (x1**2 + x2 + np.cos(x3)*t0*t1 + 1) + 1]
        >>> def lim1(x2, x3, t0, t1):
        ...     return [scale * (t0*x2 + t1*x3) - 1,
        ...             scale * (t0*x2 + t1*x3) + 1]
        >>> def lim2(x3, t0, t1):
        ...     return [scale * (x3 + t0**2*t1**3) - 1,
        ...             scale * (x3 + t0**2*t1**3) + 1]
        >>> def lim3(t0, t1):
        ...     return [scale * (t0+t1) - 1, scale * (t0+t1) + 1]
        >>> def opts0(x1, x2, x3, t0, t1):
        ...     return {'points' : [t0 - t1*x1]}
        >>> def opts1(x2, x3, t0, t1):
        ...     return {}
        >>> def opts2(x3, t0, t1):
        ...     return {}
        >>> def opts3(t0, t1):
        ...     return {}
        >>> integrate.nquad(func2, [lim0, lim1, lim2, lim3], args=(0,0),
        ...                 opts=[opts0, opts1, opts2, opts3])
        (25.066666666666666, 2.7829590483937256e-13)
    
    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)
        Integrate a system of ordinary differential equations.
        
        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, t0, ...)
        
        where y can be a vector.
        
        *Note*: The first two arguments of ``func(y, t0, ...)`` are in the
        opposite order of the arguments in the system definition function used
        by the `scipy.integrate.ode` class.
        
        Parameters
        ----------
        func : callable(y, t0, ...)
            Computes the derivative of y at t0.
        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.
        args : tuple, optional
            Extra arguments to pass to function.
        Dfun : callable(y, t0, ...)
            Gradient (Jacobian) of `func`.
        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
        
        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
        --------
        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 it initially at rest, so
        `omega(0)` = 0.  Then the vector of initial conditions is
        
        >>> y0 = [np.pi - 0.1, 0.0]
        
        We generate a solution 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()
    
    quad(func, a, b, args=(), full_output=0, epsabs=1.49e-08, epsrel=1.49e-08, limit=50, points=None, weight=None, wvar=None, wopts=None, maxp1=50, limlst=50)
        Compute a definite integral.
        
        Integrate func from `a` to `b` (possibly infinite interval) using a
        technique from the Fortran library QUADPACK.
        
        Parameters
        ----------
        func : function
            A Python function or method to integrate.  If `func` takes many
            arguments, it is integrated along the axis corresponding to the
            first argument.
            If the user desires improved integration performance, then f may
            instead be a ``ctypes`` function of the form:
        
                f(int n, double args[n]),
        
            where ``args`` is an array of function arguments and ``n`` is the
            length of ``args``. ``f.argtypes`` should be set to
            ``(c_int, c_double)``, and ``f.restype`` should be ``(c_double,)``.
        a : float
            Lower limit of integration (use -numpy.inf for -infinity).
        b : float
            Upper limit of integration (use numpy.inf for +infinity).
        args : tuple, optional
            Extra arguments to pass to `func`.
        full_output : int, optional
            Non-zero to return a dictionary of integration information.
            If non-zero, warning messages are also suppressed and the
            message is appended to the output tuple.
        
        Returns
        -------
        y : float
            The integral of func from `a` to `b`.
        abserr : float
            An estimate of the absolute error in the result.
        infodict : dict
            A dictionary containing additional information.
            Run scipy.integrate.quad_explain() for more information.
        message
            A convergence message.
        explain
            Appended only with 'cos' or 'sin' weighting and infinite
            integration limits, it contains an explanation of the codes in
            infodict['ierlst']
        
        Other Parameters
        ----------------
        epsabs : float or int, optional
            Absolute error tolerance.
        epsrel : float or int, optional
            Relative error tolerance.
        limit : float or int, optional
            An upper bound on the number of subintervals used in the adaptive
            algorithm.
        points : (sequence of floats,ints), optional
            A sequence of break points in the bounded integration interval
            where local difficulties of the integrand may occur (e.g.,
            singularities, discontinuities). The sequence does not have
            to be sorted.
        weight : float or int, optional
            String indicating weighting function. Full explanation for this
            and the remaining arguments can be found below.
        wvar : optional
            Variables for use with weighting functions.
        wopts : optional
            Optional input for reusing Chebyshev moments.
        maxp1 : float or int, optional
            An upper bound on the number of Chebyshev moments.
        limlst : int, optional
            Upper bound on the number of cycles (>=3) for use with a sinusoidal
            weighting and an infinite end-point.
        
        See Also
        --------
        dblquad : double integral
        tplquad : triple integral
        nquad : n-dimensional integrals (uses `quad` recursively)
        fixed_quad : fixed-order Gaussian quadrature
        quadrature : adaptive Gaussian quadrature
        odeint : ODE integrator
        ode : ODE integrator
        simps : integrator for sampled data
        romb : integrator for sampled data
        scipy.special : for coefficients and roots of orthogonal polynomials
        
        Notes
        -----
        
        **Extra information for quad() inputs and outputs**
        
        If full_output is non-zero, then the third output argument
        (infodict) is a dictionary with entries as tabulated below.  For
        infinite limits, the range is transformed to (0,1) and the
        optional outputs are given with respect to this transformed range.
        Let M be the input argument limit and let K be infodict['last'].
        The entries are:
        
        'neval'
            The number of function evaluations.
        'last'
            The number, K, of subintervals produced in the subdivision process.
        'alist'
            A rank-1 array of length M, the first K elements of which are the
            left end points of the subintervals in the partition of the
            integration range.
        'blist'
            A rank-1 array of length M, the first K elements of which are the
            right end points of the subintervals.
        'rlist'
            A rank-1 array of length M, the first K elements of which are the
            integral approximations on the subintervals.
        'elist'
            A rank-1 array of length M, the first K elements of which are the
            moduli of the absolute error estimates on the subintervals.
        'iord'
            A rank-1 integer array of length M, the first L elements of
            which are pointers to the error estimates over the subintervals
            with ``L=K`` if ``K<=M/2+2`` or ``L=M+1-K`` otherwise. Let I be the
            sequence ``infodict['iord']`` and let E be the sequence
            ``infodict['elist']``.  Then ``E[I[1]], ..., E[I[L]]`` forms a
            decreasing sequence.
        
        If the input argument points is provided (i.e. it is not None),
        the following additional outputs are placed in the output
        dictionary.  Assume the points sequence is of length P.
        
        'pts'
            A rank-1 array of length P+2 containing the integration limits
            and the break points of the intervals in ascending order.
            This is an array giving the subintervals over which integration
            will occur.
        'level'
            A rank-1 integer array of length M (=limit), containing the
            subdivision levels of the subintervals, i.e., if (aa,bb) is a
            subinterval of ``(pts[1], pts[2])`` where ``pts[0]`` and ``pts[2]``
            are adjacent elements of ``infodict['pts']``, then (aa,bb) has level l
            if ``|bb-aa| = |pts[2]-pts[1]| * 2**(-l)``.
        'ndin'
            A rank-1 integer array of length P+2.  After the first integration
            over the intervals (pts[1], pts[2]), the error estimates over some
            of the intervals may have been increased artificially in order to
            put their subdivision forward.  This array has ones in slots
            corresponding to the subintervals for which this happens.
        
        **Weighting the integrand**
        
        The input variables, *weight* and *wvar*, are used to weight the
        integrand by a select list of functions.  Different integration
        methods are used to compute the integral with these weighting
        functions.  The possible values of weight and the corresponding
        weighting functions are.
        
        ==========  ===================================   =====================
        ``weight``  Weight function used                  ``wvar``
        ==========  ===================================   =====================
        'cos'       cos(w*x)                              wvar = w
        'sin'       sin(w*x)                              wvar = w
        'alg'       g(x) = ((x-a)**alpha)*((b-x)**beta)   wvar = (alpha, beta)
        'alg-loga'  g(x)*log(x-a)                         wvar = (alpha, beta)
        'alg-logb'  g(x)*log(b-x)                         wvar = (alpha, beta)
        'alg-log'   g(x)*log(x-a)*log(b-x)                wvar = (alpha, beta)
        'cauchy'    1/(x-c)                               wvar = c
        ==========  ===================================   =====================
        
        wvar holds the parameter w, (alpha, beta), or c depending on the weight
        selected.  In these expressions, a and b are the integration limits.
        
        For the 'cos' and 'sin' weighting, additional inputs and outputs are
        available.
        
        For finite integration limits, the integration is performed using a
        Clenshaw-Curtis method which uses Chebyshev moments.  For repeated
        calculations, these moments are saved in the output dictionary:
        
        'momcom'
            The maximum level of Chebyshev moments that have been computed,
            i.e., if ``M_c`` is ``infodict['momcom']`` then the moments have been
            computed for intervals of length ``|b-a| * 2**(-l)``,
            ``l=0,1,...,M_c``.
        'nnlog'
            A rank-1 integer array of length M(=limit), containing the
            subdivision levels of the subintervals, i.e., an element of this
            array is equal to l if the corresponding subinterval is
            ``|b-a|* 2**(-l)``.
        'chebmo'
            A rank-2 array of shape (25, maxp1) containing the computed
            Chebyshev moments.  These can be passed on to an integration
            over the same interval by passing this array as the second
            element of the sequence wopts and passing infodict['momcom'] as
            the first element.
        
        If one of the integration limits is infinite, then a Fourier integral is
        computed (assuming w neq 0).  If full_output is 1 and a numerical error
        is encountered, besides the error message attached to the output tuple,
        a dictionary is also appended to the output tuple which translates the
        error codes in the array ``info['ierlst']`` to English messages.  The
        output information dictionary contains the following entries instead of
        'last', 'alist', 'blist', 'rlist', and 'elist':
        
        'lst'
            The number of subintervals needed for the integration (call it ``K_f``).
        'rslst'
            A rank-1 array of length M_f=limlst, whose first ``K_f`` elements
            contain the integral contribution over the interval
            ``(a+(k-1)c, a+kc)`` where ``c = (2*floor(|w|) + 1) * pi / |w|``
            and ``k=1,2,...,K_f``.
        'erlst'
            A rank-1 array of length ``M_f`` containing the error estimate
            corresponding to the interval in the same position in
            ``infodict['rslist']``.
        'ierlst'
            A rank-1 integer array of length ``M_f`` containing an error flag
            corresponding to the interval in the same position in
            ``infodict['rslist']``.  See the explanation dictionary (last entry
            in the output tuple) for the meaning of the codes.
        
        Examples
        --------
        Calculate :math:`\int^4_0 x^2 dx` and compare with an analytic result
        
        >>> from scipy import integrate
        >>> x2 = lambda x: x**2
        >>> integrate.quad(x2, 0, 4)
        (21.333333333333332, 2.3684757858670003e-13)
        >>> print(4**3 / 3.)  # analytical result
        21.3333333333
        
        Calculate :math:`\int^\infty_0 e^{-x} dx`
        
        >>> invexp = lambda x: np.exp(-x)
        >>> integrate.quad(invexp, 0, np.inf)
        (1.0, 5.842605999138044e-11)
        
        >>> f = lambda x,a : a*x
        >>> y, err = integrate.quad(f, 0, 1, args=(1,))
        >>> y
        0.5
        >>> y, err = integrate.quad(f, 0, 1, args=(3,))
        >>> y
        1.5
        
        Calculate :math:`\int^1_0 x^2 + y^2 dx` with ctypes, holding
        y parameter as 1::
        
            testlib.c =>
                double func(int n, double args[n]){
                    return args[0]*args[0] + args[1]*args[1];}
            compile to library testlib.*
        
        ::
        
           from scipy import integrate
           import ctypes
           lib = ctypes.CDLL('/home/.../testlib.*') #use absolute path
           lib.func.restype = ctypes.c_double
           lib.func.argtypes = (ctypes.c_int,ctypes.c_double)
           integrate.quad(lib.func,0,1,(1))
           #(1.3333333333333333, 1.4802973661668752e-14)
           print((1.0**3/3.0 + 1.0) - (0.0**3/3.0 + 0.0)) #Analytic result
           # 1.3333333333333333
    
    quad_explain(output=<ipykernel.iostream.OutStream object at 0x7f427da94208>)
        Print extra information about integrate.quad() parameters and returns.
        
        Parameters
        ----------
        output : instance with "write" method, optional
            Information about `quad` is passed to ``output.write()``.
            Default is ``sys.stdout``.
        
        Returns
        -------
        None
    
    quadrature(func, a, b, args=(), tol=1.49e-08, rtol=1.49e-08, maxiter=50, vec_func=True, miniter=1)
        Compute a definite integral using fixed-tolerance Gaussian quadrature.
        
        Integrate `func` from `a` to `b` using Gaussian quadrature
        with absolute tolerance `tol`.
        
        Parameters
        ----------
        func : function
            A Python function or method to integrate.
        a : float
            Lower limit of integration.
        b : float
            Upper limit of integration.
        args : tuple, optional
            Extra arguments to pass to function.
        tol, rtol : float, optional
            Iteration stops when error between last two iterates is less than
            `tol` OR the relative change is less than `rtol`.
        maxiter : int, optional
            Maximum order of Gaussian quadrature.
        vec_func : bool, optional
            True or False if func handles arrays as arguments (is
            a "vector" function). Default is True.
        miniter : int, optional
            Minimum order of Gaussian quadrature.
        
        Returns
        -------
        val : float
            Gaussian quadrature approximation (within tolerance) to integral.
        err : float
            Difference between last two estimates of the integral.
        
        See also
        --------
        romberg: adaptive Romberg quadrature
        fixed_quad: fixed-order Gaussian quadrature
        quad: adaptive quadrature using QUADPACK
        dblquad: double integrals
        tplquad: triple integrals
        romb: integrator for sampled data
        simps: integrator for sampled data
        cumtrapz: cumulative integration for sampled data
        ode: ODE integrator
        odeint: ODE integrator
    
    romb(y, dx=1.0, axis=-1, show=False)
        Romberg integration using samples of a function.
        
        Parameters
        ----------
        y : array_like
            A vector of ``2**k + 1`` equally-spaced samples of a function.
        dx : float, optional
            The sample spacing. Default is 1.
        axis : int, optional
            The axis along which to integrate. Default is -1 (last axis).
        show : bool, optional
            When `y` is a single 1-D array, then if this argument is True
            print the table showing Richardson extrapolation from the
            samples. Default is False.
        
        Returns
        -------
        romb : ndarray
            The integrated result for `axis`.
        
        See also
        --------
        quad : adaptive quadrature using QUADPACK
        romberg : adaptive Romberg quadrature
        quadrature : adaptive Gaussian quadrature
        fixed_quad : fixed-order Gaussian quadrature
        dblquad : double integrals
        tplquad : triple integrals
        simps : integrators for sampled data
        cumtrapz : cumulative integration for sampled data
        ode : ODE integrators
        odeint : ODE integrators
    
    romberg(function, a, b, args=(), tol=1.48e-08, rtol=1.48e-08, show=False, divmax=10, vec_func=False)
        Romberg integration of a callable function or method.
        
        Returns the integral of `function` (a function of one variable)
        over the interval (`a`, `b`).
        
        If `show` is 1, the triangular array of the intermediate results
        will be printed.  If `vec_func` is True (default is False), then
        `function` is assumed to support vector arguments.
        
        Parameters
        ----------
        function : callable
            Function to be integrated.
        a : float
            Lower limit of integration.
        b : float
            Upper limit of integration.
        
        Returns
        -------
        results  : float
            Result of the integration.
        
        Other Parameters
        ----------------
        args : tuple, optional
            Extra arguments to pass to function. Each element of `args` will
            be passed as a single argument to `func`. Default is to pass no
            extra arguments.
        tol, rtol : float, optional
            The desired absolute and relative tolerances. Defaults are 1.48e-8.
        show : bool, optional
            Whether to print the results. Default is False.
        divmax : int, optional
            Maximum order of extrapolation. Default is 10.
        vec_func : bool, optional
            Whether `func` handles arrays as arguments (i.e whether it is a
            "vector" function). Default is False.
        
        See Also
        --------
        fixed_quad : Fixed-order Gaussian quadrature.
        quad : Adaptive quadrature using QUADPACK.
        dblquad : Double integrals.
        tplquad : Triple integrals.
        romb : Integrators for sampled data.
        simps : Integrators for sampled data.
        cumtrapz : Cumulative integration for sampled data.
        ode : ODE integrator.
        odeint : ODE integrator.
        
        References
        ----------
        .. [1] 'Romberg's method' http://en.wikipedia.org/wiki/Romberg%27s_method
        
        Examples
        --------
        Integrate a gaussian from 0 to 1 and compare to the error function.
        
        >>> from scipy import integrate
        >>> from scipy.special import erf
        >>> gaussian = lambda x: 1/np.sqrt(np.pi) * np.exp(-x**2)
        >>> result = integrate.romberg(gaussian, 0, 1, show=True)
        Romberg integration of <function vfunc at ...> from [0, 1]
        
        ::
        
           Steps  StepSize  Results
               1  1.000000  0.385872
               2  0.500000  0.412631  0.421551
               4  0.250000  0.419184  0.421368  0.421356
               8  0.125000  0.420810  0.421352  0.421350  0.421350
              16  0.062500  0.421215  0.421350  0.421350  0.421350  0.421350
              32  0.031250  0.421317  0.421350  0.421350  0.421350  0.421350  0.421350
        
        The final result is 0.421350396475 after 33 function evaluations.
        
        >>> print("%g %g" % (2*result, erf(1)))
        0.842701 0.842701
    
    simps(y, x=None, dx=1, axis=-1, even='avg')
        Integrate y(x) using samples along the given axis and the composite
        Simpson's rule.  If x is None, spacing of dx is assumed.
        
        If there are an even number of samples, N, then there are an odd
        number of intervals (N-1), but Simpson's rule requires an even number
        of intervals.  The parameter 'even' controls how this is handled.
        
        Parameters
        ----------
        y : array_like
            Array to be integrated.
        x : array_like, optional
            If given, the points at which `y` is sampled.
        dx : int, optional
            Spacing of integration points along axis of `y`. Only used when
            `x` is None. Default is 1.
        axis : int, optional
            Axis along which to integrate. Default is the last axis.
        even : {'avg', 'first', 'str'}, optional
            'avg' : Average two results:1) use the first N-2 intervals with
                      a trapezoidal rule on the last interval and 2) use the last
                      N-2 intervals with a trapezoidal rule on the first interval.
        
            'first' : Use Simpson's rule for the first N-2 intervals with
                    a trapezoidal rule on the last interval.
        
            'last' : Use Simpson's rule for the last N-2 intervals with a
                   trapezoidal rule on the first interval.
        
        See Also
        --------
        quad: adaptive quadrature using QUADPACK
        romberg: adaptive Romberg quadrature
        quadrature: adaptive Gaussian quadrature
        fixed_quad: fixed-order Gaussian quadrature
        dblquad: double integrals
        tplquad: triple integrals
        romb: integrators for sampled data
        cumtrapz: cumulative integration for sampled data
        ode: ODE integrators
        odeint: ODE integrators
        
        Notes
        -----
        For an odd number of samples that are equally spaced the result is
        exact if the function is a polynomial of order 3 or less.  If
        the samples are not equally spaced, then the result is exact only
        if the function is a polynomial of order 2 or less.
    
    solve_bvp(fun, bc, x, y, p=None, S=None, fun_jac=None, bc_jac=None, tol=0.001, max_nodes=1000, verbose=0)
        Solve a boundary-value problem for a system of ODEs.
        
        This function numerically solves a first order system of ODEs subject to
        two-point boundary conditions::
        
            dy / dx = f(x, y, p) + S * y / (x - a), a <= x <= b
            bc(y(a), y(b), p) = 0
        
        Here x is a 1-dimensional independent variable, y(x) is a n-dimensional
        vector-valued function and p is a k-dimensional vector of unknown
        parameters which is to be found along with y(x). For the problem to be
        determined there must be n + k boundary conditions, i.e. bc must be
        (n + k)-dimensional function.
        
        The last singular term in the right-hand side of the system is optional.
        It is defined by an n-by-n matrix S, such that the solution must satisfy
        S y(a) = 0. This condition will be forced during iterations, so it must not
        contradict boundary conditions. See [2]_ for the explanation how this term
        is handled when solving BVPs numerically.
        
        Problems in a complex domain can be solved as well. In this case y and p
        are considered to be complex, and f and bc are assumed to be complex-valued
        functions, but x stays real. Note that f and bc must be complex
        differentiable (satisfy Cauchy-Riemann equations [4]_), otherwise you
        should rewrite your problem for real and imaginary parts separately. To
        solve a problem in a complex domain, pass an initial guess for y with a
        complex data type (see below).
        
        Parameters
        ----------
        fun : callable
            Right-hand side of the system. The calling signature is ``fun(x, y)``,
            or ``fun(x, y, p)`` if parameters are present. All arguments are
            ndarray: ``x`` with shape (m,), ``y`` with shape (n, m), meaning that
            ``y[:, i]`` corresponds to ``x[i]``, and ``p`` with shape (k,). The
            return value must be an array with shape (n, m) and with the same
            layout as ``y``.
        bc : callable
            Function evaluating residuals of the boundary conditions. The calling
            signature is ``bc(ya, yb)``, or ``bc(ya, yb, p)`` if parameters are
            present. All arguments are ndarray: ``ya`` and ``yb`` with shape (n,),
            and ``p`` with shape (k,). The return value must be an array with
            shape (n + k,).
        x : array_like, shape (m,)
            Initial mesh. Must be a strictly increasing sequence of real numbers
            with ``x[0]=a`` and ``x[-1]=b``.
        y : array_like, shape (n, m)
            Initial guess for the function values at the mesh nodes, i-th column
            corresponds to ``x[i]``. For problems in a complex domain pass `y`
            with a complex data type (even if the initial guess is purely real).
        p : array_like with shape (k,) or None, optional
            Initial guess for the unknown parameters. If None (default), it is
            assumed that the problem doesn't depend on any parameters.
        S : array_like with shape (n, n) or None
            Matrix defining the singular term. If None (default), the problem is
            solved without the singular term.
        fun_jac : callable or None, optional
            Function computing derivatives of f with respect to y and p. The
            calling signature is ``fun_jac(x, y)``, or ``fun_jac(x, y, p)`` if
            parameters are present. The return must contain 1 or 2 elements in the
            following order:
        
                * df_dy : array_like with shape (n, n, m) where an element
                  (i, j, q) equals to d f_i(x_q, y_q, p) / d (y_q)_j.
                * df_dp : array_like with shape (n, k, m) where an element
                  (i, j, q) equals to d f_i(x_q, y_q, p) / d p_j.
        
            Here q numbers nodes at which x and y are defined, whereas i and j
            number vector components. If the problem is solved without unknown
            parameters df_dp should not be returned.
        
            If `fun_jac` is None (default), the derivatives will be estimated
            by the forward finite differences.
        bc_jac : callable or None, optional
            Function computing derivatives of bc with respect to ya, yb and p.
            The calling signature is ``bc_jac(ya, yb)``, or ``bc_jac(ya, yb, p)``
            if parameters are present. The return must contain 2 or 3 elements in
            the following order:
        
                * dbc_dya : array_like with shape (n, n) where an element (i, j)
                  equals to d bc_i(ya, yb, p) / d ya_j.
                * dbc_dyb : array_like with shape (n, n) where an element (i, j)
                  equals to d bc_i(ya, yb, p) / d yb_j.
                * dbc_dp : array_like with shape (n, k) where an element (i, j)
                  equals to d bc_i(ya, yb, p) / d p_j.
        
            If the problem is solved without unknown parameters dbc_dp should not
            be returned.
        
            If `bc_jac` is None (default), the derivatives will be estimated by
            the forward finite differences.
        tol : float, optional
            Desired tolerance of the solution. If we define ``r = y' - f(x, y)``
            where y is the found solution, then the solver tries to achieve on each
            mesh interval ``norm(r / (1 + abs(f)) < tol``, where ``norm`` is
            estimated in a root mean squared sense (using a numerical quadrature
            formula). Default is 1e-3.
        max_nodes : int, optional
            Maximum allowed number of the mesh nodes. If exceeded, the algorithm
            terminates. Default is 1000.
        verbose : {0, 1, 2}, optional
            Level of algorithm's verbosity:
        
                * 0 (default) : work silently.
                * 1 : display a termination report.
                * 2 : display progress during iterations.
        
        Returns
        -------
        Bunch object with the following fields defined:
        sol : PPoly
            Found solution for y as `scipy.interpolate.PPoly` instance, a C1
            continuous cubic spline.
        p : ndarray or None, shape (k,)
            Found parameters. None, if the parameters were not present in the
            problem.
        x : ndarray, shape (m,)
            Nodes of the final mesh.
        y : ndarray, shape (n, m)
            Solution values at the mesh nodes.
        yp : ndarray, shape (n, m)
            Solution derivatives at the mesh nodes.
        rms_residuals : ndarray, shape (m - 1,)
            RMS values of the relative residuals over each mesh interval (see the
            description of `tol` parameter).
        niter : int
            Number of completed iterations.
        status : int
            Reason for algorithm termination:
        
                * 0: The algorithm converged to the desired accuracy.
                * 1: The maximum number of mesh nodes is exceeded.
                * 2: A singular Jacobian encountered when solving the collocation
                  system.
        
        message : string
            Verbal description of the termination reason.
        success : bool
            True if the algorithm converged to the desired accuracy (``status=0``).
        
        Notes
        -----
        This function implements a 4-th order collocation algorithm with the
        control of residuals similar to [1]_. A collocation system is solved
        by a damped Newton method with an affine-invariant criterion function as
        described in [3]_.
        
        Note that in [1]_  integral residuals are defined without normalization
        by interval lengths. So their definition is different by a multiplier of
        h**0.5 (h is an interval length) from the definition used here.
        
        .. versionadded:: 0.18.0
        
        References
        ----------
        .. [1] J. Kierzenka, L. F. Shampine, "A BVP Solver Based on Residual
               Control and the Maltab PSE", ACM Trans. Math. Softw., Vol. 27,
               Number 3, pp. 299-316, 2001.
        .. [2] L.F. Shampine, P. H. Muir and H. Xu, "A User-Friendly Fortran BVP
               Solver".
        .. [3] U. Ascher, R. Mattheij and R. Russell "Numerical Solution of
               Boundary Value Problems for Ordinary Differential Equations".
        .. [4] `Cauchy-Riemann equations
                <https://en.wikipedia.org/wiki/Cauchy-Riemann_equations>`_ on
                Wikipedia.
        
        Examples
        --------
        In the first example we solve Bratu's problem::
        
            y'' + k * exp(y) = 0
            y(0) = y(1) = 0
        
        for k = 1.
        
        We rewrite the equation as a first order system and implement its
        right-hand side evaluation::
        
            y1' = y2
            y2' = -exp(y1)
        
        >>> def fun(x, y):
        ...     return np.vstack((y[1], -np.exp(y[0])))
        
        Implement evaluation of the boundary condition residuals:
        
        >>> def bc(ya, yb):
        ...     return np.array([ya[0], yb[0]])
        
        Define the initial mesh with 5 nodes:
        
        >>> x = np.linspace(0, 1, 5)
        
        This problem is known to have two solutions. To obtain both of them we
        use two different initial guesses for y. We denote them by subscripts
        a and b.
        
        >>> y_a = np.zeros((2, x.size))
        >>> y_b = np.zeros((2, x.size))
        >>> y_b[0] = 3
        
        Now we are ready to run the solver.
        
        >>> from scipy.integrate import solve_bvp
        >>> res_a = solve_bvp(fun, bc, x, y_a)
        >>> res_b = solve_bvp(fun, bc, x, y_b)
        
        Let's plot the two found solutions. We take an advantage of having the
        solution in a spline form to produce a smooth plot.
        
        >>> x_plot = np.linspace(0, 1, 100)
        >>> y_plot_a = res_a.sol(x_plot)[0]
        >>> y_plot_b = res_b.sol(x_plot)[0]
        >>> import matplotlib.pyplot as plt
        >>> plt.plot(x_plot, y_plot_a, label='y_a')
        >>> plt.plot(x_plot, y_plot_b, label='y_b')
        >>> plt.legend()
        >>> plt.xlabel("x")
        >>> plt.ylabel("y")
        >>> plt.show()
        
        We see that the two solutions have similar shape, but differ in scale
        significantly.
        
        In the second example we solve a simple Sturm-Liouville problem::
        
            y'' + k**2 * y = 0
            y(0) = y(1) = 0
        
        It is known that a non-trivial solution y = A * sin(k * x) is possible for
        k = pi * n, where n is an integer. To establish the normalization constant
        A = 1 we add a boundary condition::
        
            y'(0) = k
        
        Again we rewrite our equation as a first order system and implement its
        right-hand side evaluation::
        
            y1' = y2
            y2' = -k**2 * y1
        
        >>> def fun(x, y, p):
        ...     k = p[0]
        ...     return np.vstack((y[1], -k**2 * y[0]))
        
        Note that parameters p are passed as a vector (with one element in our
        case).
        
        Implement the boundary conditions:
        
        >>> def bc(ya, yb, p):
        ...     k = p[0]
        ...     return np.array([ya[0], yb[0], ya[1] - k])
        
        Setup the initial mesh and guess for y. We aim to find the solution for
        k = 2 * pi, to achieve that we set values of y to approximately follow
        sin(2 * pi * x):
        
        >>> x = np.linspace(0, 1, 5)
        >>> y = np.zeros((2, x.size))
        >>> y[0, 1] = 1
        >>> y[0, 3] = -1
        
        Run the solver with 6 as an initial guess for k.
        
        >>> sol = solve_bvp(fun, bc, x, y, p=[6])
        
        We see that the found k is approximately correct:
        
        >>> sol.p[0]
        6.28329460046
        
        And finally plot the solution to see the anticipated sinusoid:
        
        >>> x_plot = np.linspace(0, 1, 100)
        >>> y_plot = sol.sol(x_plot)[0]
        >>> plt.plot(x_plot, y_plot)
        >>> plt.xlabel("x")
        >>> plt.ylabel("y")
        >>> plt.show()
    
    tplquad(func, a, b, gfun, hfun, qfun, rfun, args=(), epsabs=1.49e-08, epsrel=1.49e-08)
        Compute a triple (definite) integral.
        
        Return the triple integral of ``func(z, y, x)`` from ``x = a..b``,
        ``y = gfun(x)..hfun(x)``, and ``z = qfun(x,y)..rfun(x,y)``.
        
        Parameters
        ----------
        func : function
            A Python function or method of at least three variables in the
            order (z, y, x).
        a, b : float
            The limits of integration in x: `a` < `b`
        gfun : function
            The lower boundary curve in y which is a function taking a single
            floating point argument (x) and returning a floating point result:
            a lambda function can be useful here.
        hfun : function
            The upper boundary curve in y (same requirements as `gfun`).
        qfun : function
            The lower boundary surface in z.  It must be a function that takes
            two floats in the order (x, y) and returns a float.
        rfun : function
            The upper boundary surface in z. (Same requirements as `qfun`.)
        args : tuple, optional
            Extra arguments to pass to `func`.
        epsabs : float, optional
            Absolute tolerance passed directly to the innermost 1-D quadrature
            integration. Default is 1.49e-8.
        epsrel : float, optional
            Relative tolerance of the innermost 1-D integrals. Default is 1.49e-8.
        
        Returns
        -------
        y : float
            The resultant integral.
        abserr : float
            An estimate of the error.
        
        See Also
        --------
        quad: Adaptive quadrature using QUADPACK
        quadrature: Adaptive Gaussian quadrature
        fixed_quad: Fixed-order Gaussian quadrature
        dblquad: Double integrals
        nquad : N-dimensional integrals
        romb: Integrators for sampled data
        simps: Integrators for sampled data
        ode: ODE integrators
        odeint: ODE integrators
        scipy.special: For coefficients and roots of orthogonal polynomials
    
    trapz(y, x=None, dx=1.0, axis=-1)
        Integrate along the given axis using the composite trapezoidal rule.
        
        Integrate `y` (`x`) along given axis.
        
        Parameters
        ----------
        y : array_like
            Input array to integrate.
        x : array_like, optional
            The sample points corresponding to the `y` values. If `x` is None,
            the sample points are assumed to be evenly spaced `dx` apart. The
            default is None.
        dx : scalar, optional
            The spacing between sample points when `x` is None. The default is 1.
        axis : int, optional
            The axis along which to integrate.
        
        Returns
        -------
        trapz : float
            Definite integral as approximated by trapezoidal rule.
        
        See Also
        --------
        sum, cumsum
        
        Notes
        -----
        Image [2]_ illustrates trapezoidal rule -- y-axis locations of points
        will be taken from `y` array, by default x-axis distances between
        points will be 1.0, alternatively they can be provided with `x` array
        or with `dx` scalar.  Return value will be equal to combined area under
        the red lines.
        
        
        References
        ----------
        .. [1] Wikipedia page: http://en.wikipedia.org/wiki/Trapezoidal_rule
        
        .. [2] Illustration image:
               http://en.wikipedia.org/wiki/File:Composite_trapezoidal_rule_illustration.png
        
        Examples
        --------
        >>> np.trapz([1,2,3])
        4.0
        >>> np.trapz([1,2,3], x=[4,6,8])
        8.0
        >>> np.trapz([1,2,3], dx=2)
        8.0
        >>> a = np.arange(6).reshape(2, 3)
        >>> a
        array([[0, 1, 2],
               [3, 4, 5]])
        >>> np.trapz(a, axis=0)
        array([ 1.5,  2.5,  3.5])
        >>> np.trapz(a, axis=1)
        array([ 2.,  8.])

DATA
    __all__ = ['IntegrationWarning', 'absolute_import', 'complex_ode', 'cu...
    absolute_import = _Feature((2, 5, 0, 'alpha', 1), (3, 0, 0, 'alpha', 0...
    division = _Feature((2, 2, 0, 'alpha', 2), (3, 0, 0, 'alpha', 0), 8192...
    print_function = _Feature((2, 6, 0, 'alpha', 2), (3, 0, 0, 'alpha', 0)...

FILE
    /home/zingale/.local/lib/python3.5/site-packages/scipy/integrate/__init__.py


quad is the basic integrator for a general (not sampled) function. It uses a general-interface from the Fortran package QUADPACK (QAGS or QAGI). It will return the integral in an interval and an estimate of the error in the approximation


In [4]:
def f(x):
    return np.sin(x)**2

In [5]:
I, err = integrate.quad(f, 0.0, 2.0*np.pi, epsabs=1.e-14)
print(I)
print(err)


3.141592653589793
2.3058791671639882e-09

In [6]:
help(integrate.quad)


Help on function quad in module scipy.integrate.quadpack:

quad(func, a, b, args=(), full_output=0, epsabs=1.49e-08, epsrel=1.49e-08, limit=50, points=None, weight=None, wvar=None, wopts=None, maxp1=50, limlst=50)
    Compute a definite integral.
    
    Integrate func from `a` to `b` (possibly infinite interval) using a
    technique from the Fortran library QUADPACK.
    
    Parameters
    ----------
    func : function
        A Python function or method to integrate.  If `func` takes many
        arguments, it is integrated along the axis corresponding to the
        first argument.
        If the user desires improved integration performance, then f may
        instead be a ``ctypes`` function of the form:
    
            f(int n, double args[n]),
    
        where ``args`` is an array of function arguments and ``n`` is the
        length of ``args``. ``f.argtypes`` should be set to
        ``(c_int, c_double)``, and ``f.restype`` should be ``(c_double,)``.
    a : float
        Lower limit of integration (use -numpy.inf for -infinity).
    b : float
        Upper limit of integration (use numpy.inf for +infinity).
    args : tuple, optional
        Extra arguments to pass to `func`.
    full_output : int, optional
        Non-zero to return a dictionary of integration information.
        If non-zero, warning messages are also suppressed and the
        message is appended to the output tuple.
    
    Returns
    -------
    y : float
        The integral of func from `a` to `b`.
    abserr : float
        An estimate of the absolute error in the result.
    infodict : dict
        A dictionary containing additional information.
        Run scipy.integrate.quad_explain() for more information.
    message
        A convergence message.
    explain
        Appended only with 'cos' or 'sin' weighting and infinite
        integration limits, it contains an explanation of the codes in
        infodict['ierlst']
    
    Other Parameters
    ----------------
    epsabs : float or int, optional
        Absolute error tolerance.
    epsrel : float or int, optional
        Relative error tolerance.
    limit : float or int, optional
        An upper bound on the number of subintervals used in the adaptive
        algorithm.
    points : (sequence of floats,ints), optional
        A sequence of break points in the bounded integration interval
        where local difficulties of the integrand may occur (e.g.,
        singularities, discontinuities). The sequence does not have
        to be sorted.
    weight : float or int, optional
        String indicating weighting function. Full explanation for this
        and the remaining arguments can be found below.
    wvar : optional
        Variables for use with weighting functions.
    wopts : optional
        Optional input for reusing Chebyshev moments.
    maxp1 : float or int, optional
        An upper bound on the number of Chebyshev moments.
    limlst : int, optional
        Upper bound on the number of cycles (>=3) for use with a sinusoidal
        weighting and an infinite end-point.
    
    See Also
    --------
    dblquad : double integral
    tplquad : triple integral
    nquad : n-dimensional integrals (uses `quad` recursively)
    fixed_quad : fixed-order Gaussian quadrature
    quadrature : adaptive Gaussian quadrature
    odeint : ODE integrator
    ode : ODE integrator
    simps : integrator for sampled data
    romb : integrator for sampled data
    scipy.special : for coefficients and roots of orthogonal polynomials
    
    Notes
    -----
    
    **Extra information for quad() inputs and outputs**
    
    If full_output is non-zero, then the third output argument
    (infodict) is a dictionary with entries as tabulated below.  For
    infinite limits, the range is transformed to (0,1) and the
    optional outputs are given with respect to this transformed range.
    Let M be the input argument limit and let K be infodict['last'].
    The entries are:
    
    'neval'
        The number of function evaluations.
    'last'
        The number, K, of subintervals produced in the subdivision process.
    'alist'
        A rank-1 array of length M, the first K elements of which are the
        left end points of the subintervals in the partition of the
        integration range.
    'blist'
        A rank-1 array of length M, the first K elements of which are the
        right end points of the subintervals.
    'rlist'
        A rank-1 array of length M, the first K elements of which are the
        integral approximations on the subintervals.
    'elist'
        A rank-1 array of length M, the first K elements of which are the
        moduli of the absolute error estimates on the subintervals.
    'iord'
        A rank-1 integer array of length M, the first L elements of
        which are pointers to the error estimates over the subintervals
        with ``L=K`` if ``K<=M/2+2`` or ``L=M+1-K`` otherwise. Let I be the
        sequence ``infodict['iord']`` and let E be the sequence
        ``infodict['elist']``.  Then ``E[I[1]], ..., E[I[L]]`` forms a
        decreasing sequence.
    
    If the input argument points is provided (i.e. it is not None),
    the following additional outputs are placed in the output
    dictionary.  Assume the points sequence is of length P.
    
    'pts'
        A rank-1 array of length P+2 containing the integration limits
        and the break points of the intervals in ascending order.
        This is an array giving the subintervals over which integration
        will occur.
    'level'
        A rank-1 integer array of length M (=limit), containing the
        subdivision levels of the subintervals, i.e., if (aa,bb) is a
        subinterval of ``(pts[1], pts[2])`` where ``pts[0]`` and ``pts[2]``
        are adjacent elements of ``infodict['pts']``, then (aa,bb) has level l
        if ``|bb-aa| = |pts[2]-pts[1]| * 2**(-l)``.
    'ndin'
        A rank-1 integer array of length P+2.  After the first integration
        over the intervals (pts[1], pts[2]), the error estimates over some
        of the intervals may have been increased artificially in order to
        put their subdivision forward.  This array has ones in slots
        corresponding to the subintervals for which this happens.
    
    **Weighting the integrand**
    
    The input variables, *weight* and *wvar*, are used to weight the
    integrand by a select list of functions.  Different integration
    methods are used to compute the integral with these weighting
    functions.  The possible values of weight and the corresponding
    weighting functions are.
    
    ==========  ===================================   =====================
    ``weight``  Weight function used                  ``wvar``
    ==========  ===================================   =====================
    'cos'       cos(w*x)                              wvar = w
    'sin'       sin(w*x)                              wvar = w
    'alg'       g(x) = ((x-a)**alpha)*((b-x)**beta)   wvar = (alpha, beta)
    'alg-loga'  g(x)*log(x-a)                         wvar = (alpha, beta)
    'alg-logb'  g(x)*log(b-x)                         wvar = (alpha, beta)
    'alg-log'   g(x)*log(x-a)*log(b-x)                wvar = (alpha, beta)
    'cauchy'    1/(x-c)                               wvar = c
    ==========  ===================================   =====================
    
    wvar holds the parameter w, (alpha, beta), or c depending on the weight
    selected.  In these expressions, a and b are the integration limits.
    
    For the 'cos' and 'sin' weighting, additional inputs and outputs are
    available.
    
    For finite integration limits, the integration is performed using a
    Clenshaw-Curtis method which uses Chebyshev moments.  For repeated
    calculations, these moments are saved in the output dictionary:
    
    'momcom'
        The maximum level of Chebyshev moments that have been computed,
        i.e., if ``M_c`` is ``infodict['momcom']`` then the moments have been
        computed for intervals of length ``|b-a| * 2**(-l)``,
        ``l=0,1,...,M_c``.
    'nnlog'
        A rank-1 integer array of length M(=limit), containing the
        subdivision levels of the subintervals, i.e., an element of this
        array is equal to l if the corresponding subinterval is
        ``|b-a|* 2**(-l)``.
    'chebmo'
        A rank-2 array of shape (25, maxp1) containing the computed
        Chebyshev moments.  These can be passed on to an integration
        over the same interval by passing this array as the second
        element of the sequence wopts and passing infodict['momcom'] as
        the first element.
    
    If one of the integration limits is infinite, then a Fourier integral is
    computed (assuming w neq 0).  If full_output is 1 and a numerical error
    is encountered, besides the error message attached to the output tuple,
    a dictionary is also appended to the output tuple which translates the
    error codes in the array ``info['ierlst']`` to English messages.  The
    output information dictionary contains the following entries instead of
    'last', 'alist', 'blist', 'rlist', and 'elist':
    
    'lst'
        The number of subintervals needed for the integration (call it ``K_f``).
    'rslst'
        A rank-1 array of length M_f=limlst, whose first ``K_f`` elements
        contain the integral contribution over the interval
        ``(a+(k-1)c, a+kc)`` where ``c = (2*floor(|w|) + 1) * pi / |w|``
        and ``k=1,2,...,K_f``.
    'erlst'
        A rank-1 array of length ``M_f`` containing the error estimate
        corresponding to the interval in the same position in
        ``infodict['rslist']``.
    'ierlst'
        A rank-1 integer array of length ``M_f`` containing an error flag
        corresponding to the interval in the same position in
        ``infodict['rslist']``.  See the explanation dictionary (last entry
        in the output tuple) for the meaning of the codes.
    
    Examples
    --------
    Calculate :math:`\int^4_0 x^2 dx` and compare with an analytic result
    
    >>> from scipy import integrate
    >>> x2 = lambda x: x**2
    >>> integrate.quad(x2, 0, 4)
    (21.333333333333332, 2.3684757858670003e-13)
    >>> print(4**3 / 3.)  # analytical result
    21.3333333333
    
    Calculate :math:`\int^\infty_0 e^{-x} dx`
    
    >>> invexp = lambda x: np.exp(-x)
    >>> integrate.quad(invexp, 0, np.inf)
    (1.0, 5.842605999138044e-11)
    
    >>> f = lambda x,a : a*x
    >>> y, err = integrate.quad(f, 0, 1, args=(1,))
    >>> y
    0.5
    >>> y, err = integrate.quad(f, 0, 1, args=(3,))
    >>> y
    1.5
    
    Calculate :math:`\int^1_0 x^2 + y^2 dx` with ctypes, holding
    y parameter as 1::
    
        testlib.c =>
            double func(int n, double args[n]){
                return args[0]*args[0] + args[1]*args[1];}
        compile to library testlib.*
    
    ::
    
       from scipy import integrate
       import ctypes
       lib = ctypes.CDLL('/home/.../testlib.*') #use absolute path
       lib.func.restype = ctypes.c_double
       lib.func.argtypes = (ctypes.c_int,ctypes.c_double)
       integrate.quad(lib.func,0,1,(1))
       #(1.3333333333333333, 1.4802973661668752e-14)
       print((1.0**3/3.0 + 1.0) - (0.0**3/3.0 + 0.0)) #Analytic result
       # 1.3333333333333333

sometimes our integrand function takes optional arguments


In [7]:
def g(x, A, sigma):
    return A*np.exp(-x**2/sigma**2)

In [8]:
I, err = integrate.quad(g, -1.0, 1.0, args=(1.0, 2.0))
print(I, err)


1.8451240256511698 2.0484991765669867e-14

numpy defines the inf quantity which can be used in the integration limits. We can integrate a Gaussian (we know the answer is sqrt(pi)

Note: behind the scenes, what the integration function does is do a variable transform like: $t = 1/x$. This works when one limit is $\infty$, giving

$$\int_a^b f(x) dx = \int_{1/b}^{1/a} \frac{1}{t^2} f\left (\frac{1}{t}\right) dt$$

In [9]:
I, err = integrate.quad(g, -np.inf, np.inf, args=(1.0, 1.0))
print(I, err)


1.7724538509055159 1.4202636780944923e-08

integration of a sampled function

here we integrate a function that is defined only at a sequece of points. Recall that Simpson's rule will use piecewise parabola data. Let's compute

$$I = \int_0^{2\pi} f(x_i) dx$$

with $x_i = 0, \ldots, 2\pi$ defined at $N$ points


In [10]:
N = 17
x = np.linspace(0.0, 2.0*np.pi, N, endpoint=True)
y = np.sin(x)**2

I = integrate.simps(y, x)
print(I)


3.14159265359

Romberg integration is specific to equally-spaced samples, where $N = 2^k + 1$ and can be more converge faster (it uses extrapolation of coarser integration results to achieve higher accuracy)


In [11]:
N = 17
x = np.linspace(0.0, 2.0*np.pi, N, endpoint=True)
y = np.sin(x)**2

I = integrate.romb(y, dx=x[1]-x[0])
print(I)


3.14306583533

Root Finding

Often we need to find a value of a variable that zeros a function -- this is root finding. Sometimes, this is a multidimensional problem.

The brentq() routine offers a very robust method for find roots from a scalar function. You do need to provide an interval that bounds the root.

$f(x) = \frac{x e^x}{e^x - 1} - 5$


In [12]:
import scipy.optimize as optimize

def f(x):
    # this is the non-linear equation that comes up in deriving Wien's law for radiation
    return (x*np.exp(x)/(np.exp(x) - 1.0) - 5.0)

root, r = optimize.brentq(f, 0.1, 10.0, full_output=True)

print(root)
print(r.converged)


4.965114231744287
True

In [13]:
x = np.linspace(0.1, 10.0, 1000)
plt.plot(x, f(x))
plt.plot(np.array([root]), np.array([f(root)]), 'ro')


Out[13]:
[<matplotlib.lines.Line2D at 0x7f4280aa9400>]

ODEs

Many methods exist for integrating ordinary differential equations. Most will want you to write your ODEs as a system of first order equations.

This system of ODEs is the Lorenz system:

$$\frac{dx}{dt} = \sigma (y - x)$$$$\frac{dy}{dt} = rx - y - xz$$$$\frac{dz}{dt} = xy - bz$$

the steady states of this system correspond to:

$${\bf f}({\bf x}) = \left ( \sigma (y -x), rx - y -xz, xy - bz \right )^\intercal = 0$$

In [14]:
# system parameters
sigma = 10.0
b = 8./3.
r = 28.0

def rhs(t, x):
    xdot = sigma*(x[1] - x[0])
    ydot = r*x[0] - x[1] - x[0]*x[2]
    zdot = x[0]*x[1] - b*x[2]

    return np.array([xdot, ydot, zdot])


def jac(t, x):

    return np.array(
        [ [-sigma, sigma, 0.0], 
          [r - x[2], -1.0, -x[0]],
          [x[1], x[0], -b] ])


def f(x):
    return rhs(0.,x), jac(0.,x)

This class stores the integrated solution in a simple datatype


In [15]:
class IntHistory(object):
    """ a simple container to store the integrated history """
    
    def __init__(self, t=None, x=None, y=None, z=None):
        self.t = np.array(t)
        self.x = np.array(x)
        self.y = np.array(y)
        self.z = np.array(z)

In [16]:
def ode_integrate(X0, dt, tmax):
    """ integrate using the VODE method, storing the solution each dt """

    r = integrate.ode(rhs, jac).set_integrator("vode", method="adams",
                                               with_jacobian=True,
                                               atol=1.e-10, rtol=1.e-10,
                                               nsteps = 15000, order=12)

    t = 0.0
    r.set_initial_value(X0, t)

    tout = [t]
    x1out = [X0[0]]
    x2out = [X0[1]]
    x3out = [X0[2]]

    while r.successful() and r.t < tmax:
        r.integrate(r.t+dt)

        tout.append(r.t)
        x1out.append(r.y[0])
        x2out.append(r.y[1])
        x3out.append(r.y[2])


    return IntHistory(np.array(tout), np.array(x1out), np.array(x2out), np.array(x3out))

In [17]:
H1 = ode_integrate([1.0, 1.0, 20.0], 0.02, 30)

from mpl_toolkits.mplot3d import Axes3D

fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot(H1.x, H1.y, H1.z)
fig.set_size_inches(8.0,6.0)


Multi-variate root find

we can find the steady points in this system by doing a multi-variate root find on the RHS vector


In [18]:
sol1 = optimize.root(f, [1., 1., 1.], jac=True)
print(sol1.x)

sol2 = optimize.root(f, [10., 10., 10.], jac=True)
print(sol2.x)

sol3 = optimize.root(f, [-10., -10., -10.], jac=True)
print(sol3.x)


[ 0.  0.  0.]
[  8.48528137   8.48528137  27.        ]
[ -8.48528137  -8.48528137  27.        ]

In [19]:
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot(H1.x, H1.y, H1.z)

ax.scatter(sol1.x[0], sol1.x[1], sol1.x[2], marker="x", color="r")
ax.scatter(sol2.x[0], sol2.x[1], sol2.x[2], marker="x", color="r")
ax.scatter(sol3.x[0], sol3.x[1], sol3.x[2], marker="x", color="r")

ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("z")


Out[19]:
<matplotlib.text.Text at 0x7f42486ebe10>

Fitting

Fitting is used to match a model to experimental data. E.g. we have N points of $(x_i, y_i)$ with associated errors, $\sigma_i$, and we want to find a simply function that best represents the data.

Usually this means that we will need to define a metric, often called the residual, for how well our function matches the data, and then minimize this residual. Least-squares fitting is a popular formulation.

We want to fit our data to a function $Y(x, \{a_j\})$, where $a_j$ are model parameters we can adjust. We want to find the optimal $a_j$ to minimize the distance of $Y$ from our data: $$\Delta_i = Y(x_i, \{a_j\}) - y_i$$

Least-squares minimizes $\chi^2$: $$\chi^2(\{a_j\}) = \sum_{i=1}^N \left ( \frac{\Delta_i}{\sigma_i} \right )^2$$

general linear least squares

First we'll make some experimental data (a quadratic with random fashion). We use the randn() function to provide Gaussian normalized errors.


In [20]:
def y_experiment2(a1, a2, a3, sigma, x):
    """ return the experimental data in a quadratic + random fashion,                              
        with a1, a2, a3 the coefficients of the quadratic and sigma is                             
        the error.  This will be poorly matched to a linear fit for                                
        a3 != 0 """

    N = len(x)

    # randn gives samples from the "standard normal" distribution                                  
    r = np.random.randn(N)

    y = a1 + a2*x + a3*x*x + sigma*r

    return y

N = 40
sigma = 5.0*np.ones(N)

x = np.linspace(0, 100.0, N)
y = y_experiment2(2.0, 1.50, -0.02, sigma, x)

plt.scatter(x,y)
plt.errorbar(x, y, yerr=sigma, fmt=None)


/home/zingale/.local/lib/python3.5/site-packages/matplotlib/axes/_axes.py:2818: MatplotlibDeprecationWarning: Use of None object as fmt keyword argument to suppress plotting of data values is deprecated since 1.4; use the string "none" instead.
  warnings.warn(msg, mplDeprecation, stacklevel=1)
Out[20]:
<Container object of 3 artists>

In [21]:
def resid(avec, x, y, sigma):
    """ the residual function -- this is what will be minimized by the
        scipy.optimize.leastsq() routine.  avec is the parameters we
        are optimizing -- they are packed in here, so we unpack to
        begin.  (x, y) are the data points 

        scipy.optimize.leastsq() minimizes:

           x = arg min(sum(func(y)**2,axis=0))
                    y

        so this should just be the distance from a point to the curve,
        and it will square it and sum over the points
        """

    a0, a1, a2 = avec
    
    return (y - (a0 + a1*x + a2*x**2))/sigma


# initial guesses
a0, a1, a2 = 1, 1, 1

afit, flag = optimize.leastsq(resid, [a0, a1, a2], args=(x, y, sigma))

print(afit)

plt.plot(x, afit[0] + afit[1]*x + afit[2]*x*x )
plt.scatter(x,y)
plt.errorbar(x, y, yerr=sigma, fmt=None)


[ 4.16669496  1.40563242 -0.01925657]
/home/zingale/.local/lib/python3.5/site-packages/matplotlib/axes/_axes.py:2818: MatplotlibDeprecationWarning: Use of None object as fmt keyword argument to suppress plotting of data values is deprecated since 1.4; use the string "none" instead.
  warnings.warn(msg, mplDeprecation, stacklevel=1)
Out[21]:
<Container object of 3 artists>

$\chi^2$


In [22]:
chisq = sum(np.power(resid(afit, x, y, sigma),2))
normalization = len(x)-len(afit)
print(chisq/normalization)


0.937643967456

a nonlinear example

our experiemental data -- an exponential


In [23]:
a0 = 2.5
a1 = 2./3.
sigma = 5.0

a0_orig, a1_orig = a0, a1

x = np.linspace(0.0, 4.0, 25)
y = a0*np.exp(a1*x) + sigma*np.random.randn(len(x))

plt.scatter(x,y)
plt.errorbar(x, y, yerr=sigma, fmt=None, label="_nolegend_")


/home/zingale/.local/lib/python3.5/site-packages/matplotlib/axes/_axes.py:2818: MatplotlibDeprecationWarning: Use of None object as fmt keyword argument to suppress plotting of data values is deprecated since 1.4; use the string "none" instead.
  warnings.warn(msg, mplDeprecation, stacklevel=1)
Out[23]:
<Container object of 3 artists>

our function to minimize


In [24]:
def resid(avec, x, y):
    """ the residual function -- this is what will be minimized by the                             
        scipy.optimize.leastsq() routine.  avec is the parameters we                               
        are optimizing -- they are packed in here, so we unpack to                                 
        begin.  (x, y) are the data points                                                         
                                                                                                   
        scipy.optimize.leastsq() minimizes:                                                        
                                                                                                   
           x = arg min(sum(func(y)**2,axis=0))                                                     
                    y                                                                              
                                                                                                   
        so this should just be the distance from a point to the curve,                             
        and it will square it and sum over the points                                              
        """

    a0, a1 = avec

    # note: if we wanted to deal with error bars, we would weight each                             
    # residual accordingly                                                                         
    return y - a0*np.exp(a1*x)

In [25]:
a0, a1 = 1, 1
afit, flag = optimize.leastsq(resid, [a0, a1], args=(x, y))

print(flag)
print(afit)


1
[ 2.21564232  0.72038716]

In [26]:
plt.plot(x, afit[0]*np.exp(afit[1]*x),
           label=r"$a_0 = $ %f; $a_1 = $ %f" % (afit[0], afit[1]))
plt.plot(x, a0_orig*np.exp(a1_orig*x), ":", label="original function")
plt.legend(numpoints=1, frameon=False)
plt.scatter(x,y, c="k")
plt.errorbar(x, y, yerr=sigma, fmt=None, label="_nolegend_")


/home/zingale/.local/lib/python3.5/site-packages/matplotlib/axes/_axes.py:2818: MatplotlibDeprecationWarning: Use of None object as fmt keyword argument to suppress plotting of data values is deprecated since 1.4; use the string "none" instead.
  warnings.warn(msg, mplDeprecation, stacklevel=1)
Out[26]:
<Container object of 3 artists>

FFTs

Fourier transforms convert a physical-space (or time series) representation of a function into frequency space. This provides an equivalent representation of the data with a new view.

The FFT and its inverse in NumPy use: $$F_k = \sum_{n=0}^{N-1} f_n e^{-2\pi i nk/N}$$

$$f_n = \frac{1}{N} \sum_{k=0}^{N-1} F_k e^{2\pi i n k/N}$$

Both NumPy and SciPy have FFT routines that are similar. However, the NumPy version returns the data in a more convenient form.

It's always best to start with something you understand -- let's do a simple sine wave. Since our function is real, we can use the rfft routines in NumPy -- the understand that we are working with real data and they don't return the negative frequency components.

One important caveat -- FFTs assume that you are periodic. If you include both endpoints of the domain in the points that comprise your sample then you will not match this assumption. Here we use endpoint=False with linspace()


In [27]:
def single_freq_sine(npts):

    # a pure sine with no phase shift will result in pure imaginary                                         
    # signal                                                                                                
    f_0 = 0.2

    xmax = 10.0/f_0
    
    xx = np.linspace(0.0, xmax, npts, endpoint=False)

    f = np.sin(2.0*np.pi*f_0*xx)

    return xx, f

To make our life easier, we'll define a function that plots all the stages of the FFT process


In [43]:
def plot_FFT(xx, f):

    npts = len(xx)

    # Forward transform: f(x) -> F(k)                                                                       
    fk = np.fft.rfft(f)

    # Normalization -- the '2' here comes from the fact that we are                                         
    # neglecting the negative portion of the frequency space, since                                         
    # the FFT of a real function contains redundant information, so                                         
    # we are only dealing with 1/2 of the frequency space.                                                  
    #                                                                                                       
    # technically, we should only scale the 0 bin by N, since k=0 is                                        
    # not duplicated -- we won't worry about that for these plots                                           
    norm = 2.0/npts

    fk = fk*norm

    fk_r = fk.real
    fk_i = fk.imag

    # the fftfreq returns the postive and negative (and 0) frequencies                                      
    # the newer versions of numpy (>=1.8) have an rfftfreq() function                                       
    # that really does what we want -- we'll use that here.                                                 
    k = np.fft.rfftfreq(npts)

    # to make these dimensional, we need to divide by dx.  Note that                                        
    # max(xx) is not the true length, since we didn't have a point                                          
    # at the endpoint of the domain.                                                                        
    kfreq = k*npts/(max(xx) + xx[1])

    # Inverse transform: F(k) -> f(x) -- without the normalization                                          
    fkinv = np.fft.irfft(fk/norm)

    # plots
    plt.subplot(411)

    plt.plot(xx, f)
    plt.xlabel("x")
    plt.ylabel("f(x)")

    plt.subplot(412)

    plt.plot(kfreq, fk_r, label=r"Re($\mathcal{F}$)")
    plt.plot(kfreq, fk_i, ls=":", label=r"Im($\mathcal{F}$)")
    plt.xlabel(r"$\nu_k$")
    plt.ylabel("F(k)")

    plt.legend(fontsize="small", frameon=False)

    plt.subplot(413)

    plt.plot(kfreq, np.abs(fk))
    plt.xlabel(r"$\nu_k$")
    plt.ylabel(r"|F(k)|")

    plt.subplot(414)

    plt.plot(xx, fkinv.real)
    plt.xlabel(r"$\nu_k$")
    plt.ylabel(r"inverse F(k)")

    f = plt.gcf()
    
    f.set_size_inches(10,8)
    plt.tight_layout()

In [44]:
npts = 128
xx, f = single_freq_sine(npts)
plot_FFT(xx, f)


A cosine is shifted in phase by pi/2


In [45]:
def single_freq_cosine(npts):

    # a pure cosine with no phase shift will result in pure real                                            
    # signal                                                                                                
    f_0 = 0.2

    xmax = 10.0/f_0

    xx = np.linspace(0.0, xmax, npts, endpoint=False)

    f = np.cos(2.0*np.pi*f_0*xx)

    return xx, f

In [46]:
xx, f = single_freq_cosine(npts)
plot_FFT(xx, f)


Now let's look at a sine with a pi/4 phase shift


In [47]:
def single_freq_sine_plus_shift(npts):

    # a pure sine with no phase shift will result in pure imaginary                                         
    # signal                                                                                                
    f_0 = 0.2

    xmax = 10.0/f_0

    xx = np.linspace(0.0, xmax, npts, endpoint=False)

    f = np.sin(2.0*np.pi*f_0*xx + np.pi/4)

    return xx, f

In [48]:
xx, f = single_freq_sine_plus_shift(npts)
plot_FFT(xx, f)


A frequency filter

we'll setup a simple two-frequency sine wave and filter a component


In [49]:
def two_freq_sine(npts):

    # a pure sine with no phase shift will result in pure imaginary             
    # signal                                                                    
    f_0 = 0.2
    f_1 = 0.5

    xmax = 10.0/f_0

    # we call with endpoint=False -- if we include the endpoint, then for       
    # a periodic function, the first and last point are identical -- this       
    # shows up as a signal in the FFT.                                          
    xx = np.linspace(0.0, xmax, npts, endpoint=False)

    f = 0.5*(np.sin(2.0*np.pi*f_0*xx) + np.sin(2.0*np.pi*f_1*xx))

    return xx, f

In [50]:
npts = 256

xx, f = two_freq_sine(npts)

plt.plot(xx, f)


Out[50]:
[<matplotlib.lines.Line2D at 0x7f0604e97208>]

we'll take the transform: f(x) -> F(k)


In [51]:
# normalization factor: the 2 here comes from the fact that we neglect          
# the negative portion of frequency space because our input function            
# is real                                                                       
norm = 2.0/npts
fk = norm*np.fft.rfft(f)

ofk_r = fk.real.copy()
ofk_i = fk.imag.copy()

# get the frequencies
k = np.fft.rfftfreq(len(xx))

# since we don't include the endpoint in xx, to normalize things, we need       
# max(xx) + dx to get the true length of the domain
#
# This makes the frequencies essentially multiples of 1/dx
kfreq = k*npts/(max(xx) + xx[1])


plt.plot(kfreq, fk.real, label="real")
plt.plot(kfreq, fk.imag, ":", label="imaginary")
plt.legend(frameon=False)


Out[51]:
<matplotlib.legend.Legend at 0x7f0604ec8be0>

Filter out the higher frequencies


In [52]:
fk[kfreq > 0.4] = 0.0

# element 0 of fk is the DC component                                           
fk_r = fk.real
fk_i = fk.imag


# Inverse transform: F(k) -> f(x)                                               
fkinv = np.fft.irfft(fk/norm)

plt.plot(xx, fkinv.real)


Out[52]:
[<matplotlib.lines.Line2D at 0x7f0604e85630>]

Linear Algebra

general manipulations of matrices

you can use regular NumPy arrays or you can use a special matrix class that offers some short cuts


In [53]:
a = np.array([[1.0, 2.0], [3.0, 4.0]])

In [54]:
print(a)
print(a.transpose())
print(a.T)


[[ 1.  2.]
 [ 3.  4.]]
[[ 1.  3.]
 [ 2.  4.]]
[[ 1.  3.]
 [ 2.  4.]]

In [55]:
ainv = np.linalg.inv(a)
print(ainv)


[[-2.   1. ]
 [ 1.5 -0.5]]

In [56]:
print(np.dot(a, ainv))


[[  1.00000000e+00   1.11022302e-16]
 [  0.00000000e+00   1.00000000e+00]]

the eye() function will generate an identity matrix (as will the identity())


In [57]:
print(np.eye(2))
print(np.identity(2))


[[ 1.  0.]
 [ 0.  1.]]
[[ 1.  0.]
 [ 0.  1.]]

we can solve Ax = b


In [58]:
b = np.array([5, 7])
x = np.linalg.solve(a, b)
print(x)


[-3.  4.]

The matrix class


In [59]:
A = np.matrix('1.0 2.0; 3.0 4.0')
print(A)
print(A.T)


[[ 1.  2.]
 [ 3.  4.]]
[[ 1.  3.]
 [ 2.  4.]]

In [60]:
X = np.matrix('5.0 7.0')
Y = X.T

print(A*Y)


[[ 19.]
 [ 43.]]

In [61]:
print(A.I*Y)


[[-3.]
 [ 4.]]