14. scipy-basics

In [1]:

import math
import numpy as np
%pylab inline

%pylab is deprecated, use %matplotlib inline and import the required libraries.
Populating the interactive namespace from numpy and matplotlib

Integration¶

We'll do some integrals of the form $$I = \int_a^b f(x) dx$$

In [2]:

from scipy import integrate
integrate?

Type:        module
String form: <module 'scipy.integrate' from 'c:\\Users\\WIN10HX64\\AppData\\Local\\Programs\\Python\\Python310\\lib\\site-packages\\scipy\\integrate\\__init__.py'>
File:        c:\users\win10hx64\appdata\local\programs\python\python310\lib\site-packages\scipy\integrate\__init__.py
Docstring:  
=============================================
Integration and ODEs (:mod:`scipy.integrate`)
=============================================

.. currentmodule:: scipy.integrate

Integrating functions, given function object
============================================

.. autosummary::
   :toctree: generated/

   quad          -- General purpose integration
   quad_vec      -- General purpose integration of vector-valued functions
   dblquad       -- General purpose double integration
   tplquad       -- General purpose triple integration
   nquad         -- General purpose N-D 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
   newton_cotes  -- Weights and error coefficient for Newton-Cotes integration
   IntegrationWarning -- Warning on issues during integration
   AccuracyWarning  -- Warning on issues during quadrature integration

Integrating functions, given fixed samples
==========================================

.. autosummary::
   :toctree: generated/

   trapezoid            -- Use trapezoidal rule to compute integral.
   cumulative_trapezoid -- Use trapezoidal rule to cumulatively compute integral.
   simpson              -- 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.

Solving initial value problems for ODE systems
==============================================

The solvers are implemented as individual classes, which can be used directly
(low-level usage) or through a convenience function.

.. autosummary::
   :toctree: generated/

   solve_ivp     -- Convenient function for ODE integration.
   RK23          -- Explicit Runge-Kutta solver of order 3(2).
   RK45          -- Explicit Runge-Kutta solver of order 5(4).
   DOP853        -- Explicit Runge-Kutta solver of order 8.
   Radau         -- Implicit Runge-Kutta solver of order 5.
   BDF           -- Implicit multi-step variable order (1 to 5) solver.
   LSODA         -- LSODA solver from ODEPACK Fortran package.
   OdeSolver     -- Base class for ODE solvers.
   DenseOutput   -- Local interpolant for computing a dense output.
   OdeSolution   -- Class which represents a continuous ODE solution.


Old API
-------

These are the routines developed earlier for SciPy. They wrap older solvers
implemented in Fortran (mostly ODEPACK). While the interface to them is not
particularly convenient and certain features are missing compared to the new
API, the solvers themselves are of good quality and work fast as compiled
Fortran code. In some cases, it might be worth using this old API.

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


Solving boundary value problems for ODE systems
===============================================

.. autosummary::
   :toctree: generated/

   solve_bvp     -- Solve a boundary value problem for a system of ODEs.

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 [3]:

def f(x):
    return x/sqrt(x**4 + 10*x**2 - 96*x - 71)

In [4]:

I, err = integrate.quad(f, -2, -1, epsabs=1.e-14)
print(I)
print(err)

-0.15198691759900226
1.0458679567595425e-12

In [5]:

integrate.quad?

Signature:
integrate.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,
)
Docstring:
Compute a definite integral.

Integrate func from `a` to `b` (possibly infinite interval) using a
technique from the Fortran library QUADPACK.

Parameters
----------
func : {function, scipy.LowLevelCallable}
    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
    be a `scipy.LowLevelCallable` with one of the signatures::

        double func(double x)
        double func(double x, void *user_data)
        double func(int n, double *xx)
        double func(int n, double *xx, void *user_data)

    The ``user_data`` is the data contained in the `scipy.LowLevelCallable`.
    In the call forms with ``xx``,  ``n`` is the length of the ``xx``
    array which contains ``xx[0] == x`` and the rest of the items are
    numbers contained in the ``args`` argument of quad.

    In addition, certain ctypes call signatures are supported for
    backward compatibility, but those should not be used in new code.
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.
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. Default is 1.49e-8. `quad` tries to obtain
    an accuracy of ``abs(i-result) <= max(epsabs, epsrel*abs(i))``
    where ``i`` = integral of `func` from `a` to `b`, and ``result`` is the
    numerical approximation. See `epsrel` below.
epsrel : float or int, optional
    Relative error tolerance. Default is 1.49e-8.
    If ``epsabs <= 0``, `epsrel` must be greater than both 5e-29
    and ``50 * (machine epsilon)``. See `epsabs` above.
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. Note that this option cannot be used in conjunction
    with ``weight``.
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
simpson : 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, and these do not support specifying break points. 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.


**Details of QUADPACK level routines**

`quad` calls routines from the FORTRAN library QUADPACK. This section
provides details on the conditions for each routine to be called and a
short description of each routine. The routine called depends on
`weight`, `points` and the integration limits `a` and `b`.

================  ==============  ==========  =====================
QUADPACK routine  `weight`        `points`    infinite bounds
================  ==============  ==========  =====================
qagse             None            No          No
qagie             None            No          Yes
qagpe             None            Yes         No
qawoe             'sin', 'cos'    No          No
qawfe             'sin', 'cos'    No          either `a` or `b`
qawse             'alg*'          No          No
qawce             'cauchy'        No          No
================  ==============  ==========  =====================

The following provides a short desciption from [1]_ for each
routine.

qagse
    is an integrator based on globally adaptive interval
    subdivision in connection with extrapolation, which will
    eliminate the effects of integrand singularities of
    several types.
qagie
    handles integration over infinite intervals. The infinite range is
    mapped onto a finite interval and subsequently the same strategy as
    in ``QAGS`` is applied.
qagpe
    serves the same purposes as QAGS, but also allows the
    user to provide explicit information about the location
    and type of trouble-spots i.e. the abscissae of internal
    singularities, discontinuities and other difficulties of
    the integrand function.
qawoe
    is an integrator for the evaluation of
    :math:`\int^b_a \cos(\omega x)f(x)dx` or
    :math:`\int^b_a \sin(\omega x)f(x)dx`
    over a finite interval [a,b], where :math:`\omega` and :math:`f`
    are specified by the user. The rule evaluation component is based
    on the modified Clenshaw-Curtis technique

    An adaptive subdivision scheme is used in connection
    with an extrapolation procedure, which is a modification
    of that in ``QAGS`` and allows the algorithm to deal with
    singularities in :math:`f(x)`.
qawfe
    calculates the Fourier transform
    :math:`\int^\infty_a \cos(\omega x)f(x)dx` or
    :math:`\int^\infty_a \sin(\omega x)f(x)dx`
    for user-provided :math:`\omega` and :math:`f`. The procedure of
    ``QAWO`` is applied on successive finite intervals, and convergence
    acceleration by means of the :math:`\varepsilon`-algorithm is applied
    to the series of integral approximations.
qawse
    approximate :math:`\int^b_a w(x)f(x)dx`, with :math:`a < b` where
    :math:`w(x) = (x-a)^{\alpha}(b-x)^{\beta}v(x)` with
    :math:`\alpha,\beta > -1`, where :math:`v(x)` may be one of the
    following functions: :math:`1`, :math:`\log(x-a)`, :math:`\log(b-x)`,
    :math:`\log(x-a)\log(b-x)`.

    The user specifies :math:`\alpha`, :math:`\beta` and the type of the
    function :math:`v`. A globally adaptive subdivision strategy is
    applied, with modified Clenshaw-Curtis integration on those
    subintervals which contain `a` or `b`.
qawce
    compute :math:`\int^b_a f(x) / (x-c)dx` where the integral must be
    interpreted as a Cauchy principal value integral, for user specified
    :math:`c` and :math:`f`. The strategy is globally adaptive. Modified
    Clenshaw-Curtis integration is used on those intervals containing the
    point :math:`x = c`.

References
----------

.. [1] Piessens, Robert; de Doncker-Kapenga, Elise;
       Überhuber, Christoph W.; Kahaner, David (1983).
       QUADPACK: A subroutine package for automatic integration.
       Springer-Verlag.
       ISBN 978-3-540-12553-2.

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)

Calculate :math:`\int^1_0 a x \,dx` for :math:`a = 1, 3`

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

Be aware that pulse shapes and other sharp features as compared to the
size of the integration interval may not be integrated correctly using
this method. A simplified example of this limitation is integrating a
y-axis reflected step function with many zero values within the integrals
bounds.

>>> y = lambda x: 1 if x<=0 else 0
>>> integrate.quad(y, -1, 1)
(1.0, 1.1102230246251565e-14)
>>> integrate.quad(y, -1, 100)
(1.0000000002199108, 1.0189464580163188e-08)
>>> integrate.quad(y, -1, 10000)
(0.0, 0.0)
File:      c:\users\win10hx64\appdata\local\programs\python\python310\lib\site-packages\scipy\integrate\_quadpack_py.py
Type:      function

Sometimes our integrand function takes optional arguments.

In [6]:

def g(x, A, sigma):
    return A*math.exp(-x**2/sigma**2)

In [7]:

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).

In [8]:

I, err = integrate.quad(g, -np.inf, np.inf, args=(1.0, 1.0))
print(I, err)

1.7724538509055159 1.4202636756659625e-08

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

Multidimensional Integrals¶

Multidimensional integration can be done with successive calls to quad(), but there are wrappers that help.

Let's compute

$$I = \int_{y=0}^{1/2} \int_{x=0}^{1-2y} xy dxdy = \frac{1}{96}$$

(this example comes from the SciPy tutorial)

Notice that the limits of integration in x depend on y.

Note the form of the function:

dblquad(f, a, b, ylo, yhi)

where f = f(y, x) -- the y argument is first.

The integral will be from: x = [a,b], and y = [ylo(x), yhi(x)].

In [9]:

def integrand(x,y):
    return x*y

def x_lower_lim(y):
    return 0
    
def x_upper_lim(y):
    return 1-2*y

# we change the definitions of x and y in this call
I, err = integrate.dblquad(integrand, 0.0, 0.5, x_lower_lim, x_upper_lim)
print(I, 1.0/I)

0.010416666666666668 95.99999999999999

If you remember the python lambda functions (one expression functions), you can do this more compactly:

In [10]:

I, err = integrate.dblquad(lambda x, y: x*y, 0.0, 0.5, lambda y: 0, lambda y: 1-2*y)
print(I)

0.010416666666666668

In [11]:

integrate.dblquad?

Signature:
integrate.dblquad(
    func,
    a,
    b,
    gfun,
    hfun,
    args=(),
    epsabs=1.49e-08,
    epsrel=1.49e-08,
)
Docstring:
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 or float
    The lower boundary curve in y which is a function taking a single
    floating point argument (x) and returning a floating point result
    or a float indicating a constant boundary curve.
hfun : callable or float
    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. ``dblquad`` tries to obtain
    an accuracy of ``abs(i-result) <= max(epsabs, epsrel*abs(i))``
    where ``i`` = inner integral of ``func(y, x)`` from ``gfun(x)``
    to ``hfun(x)``, and ``result`` is the numerical approximation.
    See `epsrel` below.
epsrel : float, optional
    Relative tolerance of the inner 1-D integrals. Default is 1.49e-8.
    If ``epsabs <= 0``, `epsrel` must be greater than both 5e-29
    and ``50 * (machine epsilon)``. See `epsabs` above.

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
simpson : integrator for sampled data
romb : integrator for sampled data
scipy.special : for coefficients and roots of orthogonal polynomials


Notes
-----

**Details of QUADPACK level routines**

`quad` calls routines from the FORTRAN library QUADPACK. This section
provides details on the conditions for each routine to be called and a
short description of each routine. For each level of integration, ``qagse``
is used for finite limits or ``qagie`` is used if either limit (or both!)
are infinite. The following provides a short description from [1]_ for each
routine.

qagse
    is an integrator based on globally adaptive interval
    subdivision in connection with extrapolation, which will
    eliminate the effects of integrand singularities of
    several types.
qagie
    handles integration over infinite intervals. The infinite range is
    mapped onto a finite interval and subsequently the same strategy as
    in ``QAGS`` is applied.

References
----------

.. [1] Piessens, Robert; de Doncker-Kapenga, Elise;
       Überhuber, Christoph W.; Kahaner, David (1983).
       QUADPACK: A subroutine package for automatic integration.
       Springer-Verlag.
       ISBN 978-3-540-12553-2.

Examples
--------
Compute the double integral of ``x * y**2`` over the box
``x`` ranging from 0 to 2 and ``y`` ranging from 0 to 1.
That is, :math:`\int^{x=2}_{x=0} \int^{y=1}_{y=0} x y^2 \,dy \,dx`.

>>> from scipy import integrate
>>> f = lambda y, x: x*y**2
>>> integrate.dblquad(f, 0, 2, 0, 1)
    (0.6666666666666667, 7.401486830834377e-15)

Calculate :math:`\int^{x=\pi/4}_{x=0} \int^{y=\cos(x)}_{y=\sin(x)} 1
\,dy \,dx`.

>>> f = lambda y, x: 1
>>> integrate.dblquad(f, 0, np.pi/4, np.sin, np.cos)
    (0.41421356237309503, 1.1083280054755938e-14)

Calculate :math:`\int^{x=1}_{x=0} \int^{y=x}_{y=2-x} a x y \,dy \,dx`
for :math:`a=1, 3`.

>>> f = lambda y, x, a: a*x*y
>>> integrate.dblquad(f, 0, 1, lambda x: x, lambda x: 2-x, args=(1,))
    (0.33333333333333337, 5.551115123125783e-15)
>>> integrate.dblquad(f, 0, 1, lambda x: x, lambda x: 2-x, args=(3,))
    (0.9999999999999999, 1.6653345369377348e-14)
File:      c:\users\win10hx64\appdata\local\programs\python\python310\lib\site-packages\scipy\integrate\_quadpack_py.py
Type:      function

Integration of a Smpled 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 [12]:

N = 17
x = np.linspace(0.0, 2.0*math.pi, N, endpoint=True)
y = np.sin(x)**2

plot(x, y, 'o')

Out[12]:

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

In [13]:

integrate.simpson?
I = integrate.simpson(y, x)
print(I)

3.141592653589793
Signature: integrate.simpson(y, x=None, dx=1.0, axis=-1, even='avg')
Docstring:
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 : float, optional
    Spacing of integration points along axis of `x`. Only used when
    `x` is None. Default is 1.
axis : int, optional
    Axis along which to integrate. Default is the last axis.
even : str {'avg', 'first', 'last'}, 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
cumulative_trapezoid : 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.

Examples
--------
>>> from scipy import integrate
>>> x = np.arange(0, 10)
>>> y = np.arange(0, 10)

>>> integrate.simpson(y, x)
40.5

>>> y = np.power(x, 3)
>>> integrate.simpson(y, x)
1642.5
>>> integrate.quad(lambda x: x**3, 0, 9)[0]
1640.25

>>> integrate.simpson(y, x, even='first')
1644.5
File:      c:\users\win10hx64\appdata\local\programs\python\python310\lib\site-packages\scipy\integrate\_quadrature.py
Type:      function

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 [14]:

N = 17
x = np.linspace(0.0, 2.0*math.pi, N, endpoint=True)
y = np.sin(x)**2

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

3.1430658353300385

Interpolation¶

The interp1d function allows for a variety of 1-d interpolation methods. It returns an object that acts as a function, which can be evaluated at any point.

In [15]:

import scipy.interpolate as interpolate

In [16]:

interpolate.interp1d?

Init signature:
interpolate.interp1d(
    x,
    y,
    kind='linear',
    axis=-1,
    copy=True,
    bounds_error=None,
    fill_value=nan,
    assume_sorted=False,
)
Docstring:     
Interpolate a 1-D function.

`x` and `y` are arrays of values used to approximate some function f:
``y = f(x)``. This class returns a function whose call method uses
interpolation to find the value of new points.

Parameters
----------
x : (N,) array_like
    A 1-D array of real values.
y : (...,N,...) array_like
    A N-D array of real values. The length of `y` along the interpolation
    axis must be equal to the length of `x`.
kind : str or int, optional
    Specifies the kind of interpolation as a string or as an integer
    specifying the order of the spline interpolator to use.
    The string has to be one of 'linear', 'nearest', 'nearest-up', 'zero',
    'slinear', 'quadratic', 'cubic', 'previous', or 'next'. 'zero',
    'slinear', 'quadratic' and 'cubic' refer to a spline interpolation of
    zeroth, first, second or third order; 'previous' and 'next' simply
    return the previous or next value of the point; 'nearest-up' and
    'nearest' differ when interpolating half-integers (e.g. 0.5, 1.5)
    in that 'nearest-up' rounds up and 'nearest' rounds down. Default
    is 'linear'.
axis : int, optional
    Specifies the axis of `y` along which to interpolate.
    Interpolation defaults to the last axis of `y`.
copy : bool, optional
    If True, the class makes internal copies of x and y.
    If False, references to `x` and `y` are used. The default is to copy.
bounds_error : bool, optional
    If True, a ValueError is raised any time interpolation is attempted on
    a value outside of the range of x (where extrapolation is
    necessary). If False, out of bounds values are assigned `fill_value`.
    By default, an error is raised unless ``fill_value="extrapolate"``.
fill_value : array-like or (array-like, array_like) or "extrapolate", optional
    - if a ndarray (or float), this value will be used to fill in for
      requested points outside of the data range. If not provided, then
      the default is NaN. The array-like must broadcast properly to the
      dimensions of the non-interpolation axes.
    - If a two-element tuple, then the first element is used as a
      fill value for ``x_new < x[0]`` and the second element is used for
      ``x_new > x[-1]``. Anything that is not a 2-element tuple (e.g.,
      list or ndarray, regardless of shape) is taken to be a single
      array-like argument meant to be used for both bounds as
      ``below, above = fill_value, fill_value``. Using a two-element tuple
      or ndarray requires ``bounds_error=False``.

      .. versionadded:: 0.17.0
    - If "extrapolate", then points outside the data range will be
      extrapolated.

      .. versionadded:: 0.17.0
assume_sorted : bool, optional
    If False, values of `x` can be in any order and they are sorted first.
    If True, `x` has to be an array of monotonically increasing values.

Attributes
----------
fill_value

Methods
-------
__call__

See Also
--------
splrep, splev
    Spline interpolation/smoothing based on FITPACK.
UnivariateSpline : An object-oriented wrapper of the FITPACK routines.
interp2d : 2-D interpolation

Notes
-----
Calling `interp1d` with NaNs present in input values results in
undefined behaviour.

Input values `x` and `y` must be convertible to `float` values like
`int` or `float`.

If the values in `x` are not unique, the resulting behavior is
undefined and specific to the choice of `kind`, i.e., changing
`kind` will change the behavior for duplicates.


Examples
--------
>>> import matplotlib.pyplot as plt
>>> from scipy import interpolate
>>> x = np.arange(0, 10)
>>> y = np.exp(-x/3.0)
>>> f = interpolate.interp1d(x, y)

>>> xnew = np.arange(0, 9, 0.1)
>>> ynew = f(xnew)   # use interpolation function returned by `interp1d`
>>> plt.plot(x, y, 'o', xnew, ynew, '-')
>>> plt.show()
Init docstring: Initialize a 1-D linear interpolation class.
File:           c:\users\win10hx64\appdata\local\programs\python\python310\lib\site-packages\scipy\interpolate\_interpolate.py
Type:           type
Subclasses:     

In [40]:

def fExact(x):
    return sin(x)*x

N = 10
x = np.linspace(0, 20, N)

y = fExact(x)

fInterp = interpolate.interp1d(x, y, kind=3)

In [41]:

# use finer points when we plot
xplot = np.linspace(0, 20, 10)

pylab.plot(x, y, "ro", label="known points")
pylab.plot(xplot, fExact(xplot), "b:", label="exact function")
pylab.plot(xplot, fInterp(xplot), "r-", label="interpolant")

pylab.legend(frameon=False, loc="best")

Out[41]:

<matplotlib.legend.Legend at 0x263a620bee0>

Multi-d Interpolation¶

Here's an example of mult-d interpolation from the official tutorial.

First we define the "answer" -- this is the true function that we will sample at a number of points and then try to use interpolation to recover.

In [19]:

def func(x, y):
    return x*(1-x)*np.cos(4*np.pi*x) * np.sin(4*np.pi*y**2)**2

Here we will use mgrid to create the grid of (x,y) where we know func exactly -- this will be for plotting. Note the fun trick here, this is not really a function, but rather something that can magically look like an array, and we index it with the start:stop:stride. If we set stride to an imaginary number, then it is interpreted as the number of points to put between the start and stop.

In [20]:

grid_x, grid_y = np.mgrid[0:1:100j, 0:1:200j]

In [21]:

print(grid_x.shape)
print(grid_y.shape)

(100, 200)
(100, 200)

Here is what the exact function looks like -- note that our function is defined in x,y, but imshow is meant for plotting an array, so the first index is the row. We take the transpose when plotting.

In [22]:

pylab.imshow(func(grid_x, grid_y).T, extent=(0,1,0,1), origin="lower")

Out[22]:

<matplotlib.image.AxesImage at 0x263a2fcf850>

Now we'll define 1000 random points where we'll sample the function.

In [23]:

points = np.random.rand(1000, 2)
values = func(points[:,0], points[:,1])

Here's what those points look like:

In [24]:

pylab.scatter(points[:,0], points[:,1], s=1)
pylab.xlim(0,1)
pylab.ylim(0,1)

Out[24]:

(0.0, 1.0)

The interpolate.griddata() function provides many ways to interpolate a collection of points into a uniform grid. There are many different interpolation methods within this function

In [25]:

interpolate.griddata?

Signature:
interpolate.griddata(
    points,
    values,
    xi,
    method='linear',
    fill_value=nan,
    rescale=False,
)
Docstring:
Interpolate unstructured D-D data.

Parameters
----------
points : 2-D ndarray of floats with shape (n, D), or length D tuple of 1-D ndarrays with shape (n,).
    Data point coordinates.
values : ndarray of float or complex, shape (n,)
    Data values.
xi : 2-D ndarray of floats with shape (m, D), or length D tuple of ndarrays broadcastable to the same shape.
    Points at which to interpolate data.
method : {'linear', 'nearest', 'cubic'}, optional
    Method of interpolation. One of

    ``nearest``
      return the value at the data point closest to
      the point of interpolation. See `NearestNDInterpolator` for
      more details.

    ``linear``
      tessellate the input point set to N-D
      simplices, and interpolate linearly on each simplex. See
      `LinearNDInterpolator` for more details.

    ``cubic`` (1-D)
      return the value determined from a cubic
      spline.

    ``cubic`` (2-D)
      return the value determined from a
      piecewise cubic, continuously differentiable (C1), and
      approximately curvature-minimizing polynomial surface. See
      `CloughTocher2DInterpolator` for more details.
fill_value : float, optional
    Value used to fill in for requested points outside of the
    convex hull of the input points. If not provided, then the
    default is ``nan``. This option has no effect for the
    'nearest' method.
rescale : bool, optional
    Rescale points to unit cube before performing interpolation.
    This is useful if some of the input dimensions have
    incommensurable units and differ by many orders of magnitude.

    .. versionadded:: 0.14.0

Returns
-------
ndarray
    Array of interpolated values.

Notes
-----

.. versionadded:: 0.9

Examples
--------

Suppose we want to interpolate the 2-D function

>>> def func(x, y):
...     return x*(1-x)*np.cos(4*np.pi*x) * np.sin(4*np.pi*y**2)**2

on a grid in [0, 1]x[0, 1]

>>> grid_x, grid_y = np.mgrid[0:1:100j, 0:1:200j]

but we only know its values at 1000 data points:

>>> rng = np.random.default_rng()
>>> points = rng.random((1000, 2))
>>> values = func(points[:,0], points[:,1])

This can be done with `griddata` -- below we try out all of the
interpolation methods:

>>> from scipy.interpolate import griddata
>>> grid_z0 = griddata(points, values, (grid_x, grid_y), method='nearest')
>>> grid_z1 = griddata(points, values, (grid_x, grid_y), method='linear')
>>> grid_z2 = griddata(points, values, (grid_x, grid_y), method='cubic')

One can see that the exact result is reproduced by all of the
methods to some degree, but for this smooth function the piecewise
cubic interpolant gives the best results:

>>> import matplotlib.pyplot as plt
>>> plt.subplot(221)
>>> plt.imshow(func(grid_x, grid_y).T, extent=(0,1,0,1), origin='lower')
>>> plt.plot(points[:,0], points[:,1], 'k.', ms=1)
>>> plt.title('Original')
>>> plt.subplot(222)
>>> plt.imshow(grid_z0.T, extent=(0,1,0,1), origin='lower')
>>> plt.title('Nearest')
>>> plt.subplot(223)
>>> plt.imshow(grid_z1.T, extent=(0,1,0,1), origin='lower')
>>> plt.title('Linear')
>>> plt.subplot(224)
>>> plt.imshow(grid_z2.T, extent=(0,1,0,1), origin='lower')
>>> plt.title('Cubic')
>>> plt.gcf().set_size_inches(6, 6)
>>> plt.show()

See Also
--------
LinearNDInterpolator :
    Piecewise linear interpolant in N dimensions.
NearestNDInterpolator :
    Nearest-neighbor interpolation in N dimensions.
CloughTocher2DInterpolator :
    Piecewise cubic, C1 smooth, curvature-minimizing interpolant in 2D.
File:      c:\users\win10hx64\appdata\local\programs\python\python310\lib\site-packages\scipy\interpolate\_ndgriddata.py
Type:      function

In [26]:

grid_z0 = interpolate.griddata(points, values, (grid_x, grid_y), method='nearest')
pylab.imshow(grid_z0.T, extent=(0,1,0,1), origin="lower")

Out[26]:

<matplotlib.image.AxesImage at 0x263a50cdd20>

In [27]:

grid_z0 = interpolate.griddata(points, values, (grid_x, grid_y), method='linear')
pylab.imshow(grid_z0.T, extent=(0,1,0,1), origin="lower")

Out[27]:

<matplotlib.image.AxesImage at 0x263a6148c40>

In [28]:

grid_z0 = interpolate.griddata(points, values, (grid_x, grid_y), method='cubic')
pylab.imshow(grid_z0.T, extent=(0,1,0,1), origin="lower")

Out[28]:

<matplotlib.image.AxesImage at 0x263a618f910>

Root Finding¶

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 [42]:

import scipy.optimize as optimize

In [44]:

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 [31]:

optimize.brentq?

Signature:
optimize.brentq(
    f,
    a,
    b,
    args=(),
    xtol=2e-12,
    rtol=8.881784197001252e-16,
    maxiter=100,
    full_output=False,
    disp=True,
)
Docstring:
Find a root of a function in a bracketing interval using Brent's method.

Uses the classic Brent's method to find a zero of the function `f` on
the sign changing interval [a , b]. Generally considered the best of the
rootfinding routines here. It is a safe version of the secant method that
uses inverse quadratic extrapolation. Brent's method combines root
bracketing, interval bisection, and inverse quadratic interpolation. It is
sometimes known as the van Wijngaarden-Dekker-Brent method. Brent (1973)
claims convergence is guaranteed for functions computable within [a,b].

[Brent1973]_ provides the classic description of the algorithm. Another
description can be found in a recent edition of Numerical Recipes, including
[PressEtal1992]_. A third description is at
http://mathworld.wolfram.com/BrentsMethod.html. It should be easy to
understand the algorithm just by reading our code. Our code diverges a bit
from standard presentations: we choose a different formula for the
extrapolation step.

Parameters
----------
f : function
    Python function returning a number. The function :math:`f`
    must be continuous, and :math:`f(a)` and :math:`f(b)` must
    have opposite signs.
a : scalar
    One end of the bracketing interval :math:`[a, b]`.
b : scalar
    The other end of the bracketing interval :math:`[a, b]`.
xtol : number, optional
    The computed root ``x0`` will satisfy ``np.allclose(x, x0,
    atol=xtol, rtol=rtol)``, where ``x`` is the exact root. The
    parameter must be nonnegative. For nice functions, Brent's
    method will often satisfy the above condition with ``xtol/2``
    and ``rtol/2``. [Brent1973]_
rtol : number, optional
    The computed root ``x0`` will satisfy ``np.allclose(x, x0,
    atol=xtol, rtol=rtol)``, where ``x`` is the exact root. The
    parameter cannot be smaller than its default value of
    ``4*np.finfo(float).eps``. For nice functions, Brent's
    method will often satisfy the above condition with ``xtol/2``
    and ``rtol/2``. [Brent1973]_
maxiter : int, optional
    If convergence is not achieved in `maxiter` iterations, an error is
    raised. Must be >= 0.
args : tuple, optional
    Containing extra arguments for the function `f`.
    `f` is called by ``apply(f, (x)+args)``.
full_output : bool, optional
    If `full_output` is False, the root is returned. If `full_output` is
    True, the return value is ``(x, r)``, where `x` is the root, and `r` is
    a `RootResults` object.
disp : bool, optional
    If True, raise RuntimeError if the algorithm didn't converge.
    Otherwise, the convergence status is recorded in any `RootResults`
    return object.

Returns
-------
x0 : float
    Zero of `f` between `a` and `b`.
r : `RootResults` (present if ``full_output = True``)
    Object containing information about the convergence. In particular,
    ``r.converged`` is True if the routine converged.

Notes
-----
`f` must be continuous.  f(a) and f(b) must have opposite signs.

Related functions fall into several classes:

multivariate local optimizers
  `fmin`, `fmin_powell`, `fmin_cg`, `fmin_bfgs`, `fmin_ncg`
nonlinear least squares minimizer
  `leastsq`
constrained multivariate optimizers
  `fmin_l_bfgs_b`, `fmin_tnc`, `fmin_cobyla`
global optimizers
  `basinhopping`, `brute`, `differential_evolution`
local scalar minimizers
  `fminbound`, `brent`, `golden`, `bracket`
N-D root-finding
  `fsolve`
1-D root-finding
  `brenth`, `ridder`, `bisect`, `newton`
scalar fixed-point finder
  `fixed_point`

References
----------
.. [Brent1973]
   Brent, R. P.,
   *Algorithms for Minimization Without Derivatives*.
   Englewood Cliffs, NJ: Prentice-Hall, 1973. Ch. 3-4.

.. [PressEtal1992]
   Press, W. H.; Flannery, B. P.; Teukolsky, S. A.; and Vetterling, W. T.
   *Numerical Recipes in FORTRAN: The Art of Scientific Computing*, 2nd ed.
   Cambridge, England: Cambridge University Press, pp. 352-355, 1992.
   Section 9.3:  "Van Wijngaarden-Dekker-Brent Method."

Examples
--------
>>> def f(x):
...     return (x**2 - 1)

>>> from scipy import optimize

>>> root = optimize.brentq(f, -2, 0)
>>> root
-1.0

>>> root = optimize.brentq(f, 0, 2)
>>> root
1.0
File:      c:\users\win10hx64\appdata\local\programs\python\python310\lib\site-packages\scipy\optimize\_zeros_py.py
Type:      function

In [45]:

x = np.linspace(0.1, 10.0, 1000)
pylab.plot(x, f(x))
pylab.plot(np.array([root]), np.array([f(root)]), 'ro')

Out[45]:

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

Fitting¶

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 [33]:

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)

pylab.scatter(x,y)
pylab.errorbar(x, y, yerr=sigma, fmt='none')

Out[33]:

<ErrorbarContainer object of 3 artists>

In [34]:

import scipy.optimize as optimize

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)

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

[ 1.75895379  1.44739172 -0.01915192]

Out[34]:

<ErrorbarContainer object of 3 artists>

$\chi^2$

In [35]:

chisq = sum(power(resid(afit, x, y, sigma),2))
normalization = len(x)-len(afit)
print(chisq/normalization)

1.2422349570213052

In [36]:

optimize.leastsq?

Signature:
optimize.leastsq(
    func,
    x0,
    args=(),
    Dfun=None,
    full_output=0,
    col_deriv=0,
    ftol=1.49012e-08,
    xtol=1.49012e-08,
    gtol=0.0,
    maxfev=0,
    epsfcn=None,
    factor=100,
    diag=None,
)
Docstring:
Minimize the sum of squares of a set of equations.

::

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

Parameters
----------
func : callable
    Should take at least one (possibly length ``N`` vector) argument and
    returns ``M`` floating point numbers. It must not return NaNs or
    fitting might fail. ``M`` must be greater than or equal to ``N``.
x0 : ndarray
    The starting estimate for the minimization.
args : tuple, optional
    Any extra arguments to func are placed in this tuple.
Dfun : callable, optional
    A function or method to compute the Jacobian of func with derivatives
    across the rows. If this is None, the Jacobian will be estimated.
full_output : bool, optional
    non-zero to return all optional outputs.
col_deriv : bool, optional
    non-zero to specify that the Jacobian function computes derivatives
    down the columns (faster, because there is no transpose operation).
ftol : float, optional
    Relative error desired in the sum of squares.
xtol : float, optional
    Relative error desired in the approximate solution.
gtol : float, optional
    Orthogonality desired between the function vector and the columns of
    the Jacobian.
maxfev : int, optional
    The maximum number of calls to the function. If `Dfun` is provided,
    then the default `maxfev` is 100*(N+1) where N is the number of elements
    in x0, otherwise the default `maxfev` is 200*(N+1).
epsfcn : float, optional
    A variable used in determining a suitable step length for the forward-
    difference approximation of the Jacobian (for Dfun=None).
    Normally the actual step length will be sqrt(epsfcn)*x
    If epsfcn is less than the machine precision, it is assumed that the
    relative errors are of the order of the machine precision.
factor : float, optional
    A parameter determining the initial step bound
    (``factor * || diag * x||``). Should be in interval ``(0.1, 100)``.
diag : sequence, optional
    N positive entries that serve as a scale factors for the variables.

Returns
-------
x : ndarray
    The solution (or the result of the last iteration for an unsuccessful
    call).
cov_x : ndarray
    The inverse of the Hessian. `fjac` and `ipvt` are used to construct an
    estimate of the Hessian. A value of None indicates a singular matrix,
    which means the curvature in parameters `x` is numerically flat. To
    obtain the covariance matrix of the parameters `x`, `cov_x` must be
    multiplied by the variance of the residuals -- see curve_fit.
infodict : dict
    a dictionary of optional outputs with the keys:

    ``nfev``
        The number of function calls
    ``fvec``
        The function evaluated at the output
    ``fjac``
        A permutation of the R matrix of a QR
        factorization of the final approximate
        Jacobian matrix, stored column wise.
        Together with ipvt, the covariance of the
        estimate can be approximated.
    ``ipvt``
        An integer array of length N which defines
        a permutation matrix, p, such that
        fjac*p = q*r, where r is upper triangular
        with diagonal elements of nonincreasing
        magnitude. Column j of p is column ipvt(j)
        of the identity matrix.
    ``qtf``
        The vector (transpose(q) * fvec).

mesg : str
    A string message giving information about the cause of failure.
ier : int
    An integer flag. If it is equal to 1, 2, 3 or 4, the solution was
    found. Otherwise, the solution was not found. In either case, the
    optional output variable 'mesg' gives more information.

See Also
--------
least_squares : Newer interface to solve nonlinear least-squares problems
    with bounds on the variables. See ``method=='lm'`` in particular.

Notes
-----
"leastsq" is a wrapper around MINPACK's lmdif and lmder algorithms.

cov_x is a Jacobian approximation to the Hessian of the least squares
objective function.
This approximation assumes that the objective function is based on the
difference between some observed target data (ydata) and a (non-linear)
function of the parameters `f(xdata, params)` ::

       func(params) = ydata - f(xdata, params)

so that the objective function is ::

       min   sum((ydata - f(xdata, params))**2, axis=0)
     params

The solution, `x`, is always a 1-D array, regardless of the shape of `x0`,
or whether `x0` is a scalar.

Examples
--------
>>> from scipy.optimize import leastsq
>>> def func(x):
...     return 2*(x-3)**2+1
>>> leastsq(func, 0)
(array([2.99999999]), 1)
File:      c:\users\win10hx64\appdata\local\programs\python\python310\lib\site-packages\scipy\optimize\_minpack_py.py
Type:      function

In [ ]: