Cálculo de integrales con SciPy¶
SciPy es paquete que incluye una colección de algoritmos matemáticos y funciones construidas sobre el paquete NumPy. En esta clase nos vamos a centrar en el cálculo de integrales definidas.
Como siempre lo primero es lo primero, importemos lo paquetes que vamos a utilizar:
In [1]:
import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
scipy.integrate¶
Este subpaquete de SciPy proporciona algunas técnicas de integración tanto de funciones como de ecuaciones diferenciales. En primer lugar importémoslo y ejecutemos la ayuda para ver cuáles son estas funciones:
In [5]:
from scipy import integrate
In [6]:
from IPython.display import HTML
HTML('<iframe src="http://docs.scipy.org/doc/scipy/reference/integrate.html#module-scipy.integrate" width="800" height="600"></iframe>')
Out[6]:
In [7]:
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
Integrating functions, given fixed samples
==========================================
.. autosummary::
:toctree: generated/
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.
PACKAGE CONTENTS
_dop
_ode
_odepack
_quadpack
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)
| 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__(...)
|
| __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)
|
| 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``.
|
| 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)``.
| `f` should return a scalar, array or list (not a tuple).
| 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.
|
| 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
| - rband : None or int
| Jacobian band width, jac[i,j] != 0 for i-lband <= j <= i+rband.
| Setting these requires your jac routine to return the jacobian
| in packed format, jac_packed[i-j+lband, j] = jac[i,j].
| - method: 'adams' or 'bdf'
| Which solver to use, Adams (non-stiff) or BDF (stiff)
| - with_jacobian : bool
| Whether to use the 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
| - rband : None or int
| Jacobian band width, jac[i,j] != 0 for i-lband <= j <= i+rband.
| Setting these requires your jac routine to return the jacobian
| in packed format, jac_packed[i-j+lband, j] = jac[i,j].
| - with_jacobian : bool
| Whether to use the 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.
| - 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: [email protected], [email protected]
|
| 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', with_jacobian=True)
| >>> 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:
| >>> r.integrate(r.t+dt)
| >>> print("%g %g" % (r.t, r.y))
|
| 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)
|
| 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) : tuple
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
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
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 averages 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)
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.
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.
e.g. if ``func = f(x0, x1, x2)``, then ``ranges[0]`` may be defined as
either ``(a, b)`` or else as ``(a, b) = range0(x1, x2)``.
args : iterable object, optional
Additional arguments ``t0, ..., tn``, required by `func`.
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.quadare 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. 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
The ``full_output`` option from `quad` is unavailable, due to the
complexity of handling the large amount of data such an option would
return for this kind of nested integration. For more information on
these options, see `quad` and `quad_explain`.
Returns
-------
result : float
The result of the integration.
abserr : float
The maximum of the estimates of the absolute error in the various
integration results.
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,{},{},{}])
(1.5267454070738633, 2.9437360001402324e-14)
>>> 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.
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.
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.
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
quad_explain(output=<IPython.kernel.zmq.iostream.OutStream object at 0x7f3f287e60f0>)
Print extra information about integrate.quad() parameters and returns.
Parameters
----------
output : instance with "write" method
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, rol : 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 : array_like, 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.
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) : tuple
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 : Arguments
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
If `x` is None, then spacing between all `y` elements is `dx`.
dx : scalar, optional
If `x` is None, spacing given by `dx` is assumed. Default is 1.
axis : int, optional
Specify the axis.
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/alex/anaconda3/lib/python3.4/site-packages/scipy/integrate/__init__.py
quad¶
Como se puede ver en la ayuda, si queremos realizar una integración numérica de una función de una variable, debemos utilizar quad (aunque también podemos usar trapz, simps... La forma de acceder a ella tal y como hemos importado el paquete sería ejecutando integrate.quad. Sin emabrgo, sería más normal importar del siguiete modo:
In [8]:
from scipy.integrate import quad
De este modo se puede usar la función quad, simplemente como quad. Pero todavía no sabemos cómo funciona, ¿te atreves a investigarlo tú?
In [9]:
help(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.
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
Quizá esta ayuda te resulte más atractiva.
¿Qué es lo primero que necesitamos hacer para integrar una función? Pues sí, la función... definamos una:
$$f(x) = x \cdot sin(x)$$In [10]:
def fun(x):
return x * np.sin(x)
Ejercicio¶
Antes de integrarla genera esta gráfica:
In [11]:
x = np.linspace(0,10,100)
y = fun(x)
plt.title('$y = x sin(x)$', fontsize = 25)
plt.plot(x,y, linewidth = 2)
x_fill = np.linspace(2,9,100)
y_fill = fun(x_fill)
plt.fill_between(x_fill, y_fill, color='gray', alpha=0.5)
plt.grid()
Ahora integremos con quad¶
Integremos la función en el intervalo $[2, 9]$. Recuerda que esto te calcula la integral, no el área:
In [12]:
value, err = quad(fun, 2, 9)
print("El resultado es: ", value, "con un error de: ", err)
El resultado es: 6.870699742283883 con un error de: 2.864870105641461e-13
Simpson y trapecio¶
Según figura en la documentación a estos métodos hay que pasarles las coordenadas de los puntos (no la función). Esto puede ser útil si no disponemos de una función, sino de una serie da valores, que por ejemplo, provienen de un experimento.
In [13]:
x = np.linspace(2,9,100)
value = integrate.trapz(fun(x), x)
print("El resultado es: ", value)
El resultado es: 6.86742266171
In [14]:
x = np.linspace(2,9,100)
value = integrate.simps(fun(x), x)
print("El resultado es: ", value)
El resultado es: 6.8705759095
####
¡Síguenos en Twitter!
###### Follow @AeroPython
###### Este notebook ha sido realizado por: Juan Luis Cano, y Álex Sáez
##### 
Curso AeroPython por Juan Luis Cano Rodriguez y Alejandro Sáez Mollejo se distribuye bajo una Licencia Creative Commons Atribución 4.0 Internacional.
Las siguientes celdas contienen configuración del Notebook
Para visualizar y utlizar los enlaces a Twitter el notebook debe ejecutarse como seguro
File > Trusted Notebook
In [1]:
# Esta celda da el estilo al notebook
from IPython.core.display import HTML
css_file = '../styles/aeropython.css'
HTML(open(css_file, "r").read())
Out[1]:
/* This template is inspired in the one used by Lorena Barba
in the numerical-mooc repository: https://github.com/numerical-mooc/numerical-mooc
We thank her work and hope you also enjoy the look of the notobooks with this style */
El estilo se ha aplicado =)