In :
from IPython.core.display import HTML
css_file = './custom.css'

Out:
In :
import sys #only needed to determine Python version number
import matplotlib #only needed to determine Matplotlib version number
import scipy # idem
import sympy # idem

print('Python version ' + sys.version)
print('Matplotlib version ' + matplotlib.__version__ )
print('Scipy version ' + scipy.__version__ )
print('Sympy version ' + sympy.__version__ )

%reset -f

Python version 3.6.2 |Anaconda custom (64-bit)| (default, Jul 20 2017, 13:51:32)
[GCC 4.4.7 20120313 (Red Hat 4.4.7-1)]
Matplotlib version 2.1.0
Scipy version 1.0.0
Sympy version 1.1.1


# M62_CM2 Calcul approché (package SciPy) VS formel (package SymPy).¶

## Calcul approché avec le module SciPy¶

Cette librairie est un ensemble très complet de modules d’algèbre linéaire, statistiques et autres algorithmes numériques. Le site de la documentation en fournit la liste : http://docs.scipy.org/doc/scipy/reference

### fsolve¶

Si on ne peut pas calculer analytiquement la solution d'une équation, on peut l'approcher numériquement comme suit:

In :
from math import cos
from scipy.optimize import fsolve

sol=fsolve(lambda x: x-cos(x), 1)
print(sol)

0.7390851332151607

In :
%matplotlib inline
from matplotlib.pylab import *
x=linspace(-pi,pi,101)
y=[cos(xi) for xi in x ]
plot(x,x,'g-',x,y,'r-')
plot( [sol,sol,0],[0,sol,sol] ,'-.')
grid(True) In :
from math import exp
from matplotlib.pylab import *
from scipy.optimize import fsolve

sol=fsolve(lambda x: x**2-exp(-x*x), -1)
print(sol)

sol=fsolve(lambda x: x**2-exp(-x*x), 0.2)
print(sol)

x=linspace(-1.5,1.5,101)
plot(x,[xi**2 for xi in x],'r--',label=r'$x^2$')
plot(x,[exp(-xi*xi) for xi in x],label=r'$e^{-x^2}$')

# figure()
# plot(x,[xi**2-exp(-xi*xi) for xi in x],x,[0 for xi in x])

[-0.75308916]
[0.75308916] ### integrate.quad¶

Pour approcher la valeur numérique d'une intégrale on peut utiliser integrate.quad https://docs.scipy.org/doc/scipy/reference/tutorial/integrate.html

In :
from math import *
from scipy import integrate

a=0
b=1
f = lambda x:x**2
print("Integrale =",integr, " Erreur =",integr )

f = lambda x: exp(-x*x)
print(sqrt(pi)/2)

Integrale = 0.33333333333333337  Erreur = 3.700743415417189e-15
0.8862269254527579
0.8862269254527579

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

dblquad       -- General purpose double integration
tplquad       -- General purpose triple integration
nquad         -- General purpose n-dimensional integration
romberg       -- Integrate func using Romberg integration
newton_cotes  -- Weights and error coefficient for Newton-Cotes integration
IntegrationWarning -- Warning on issues during integration

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

.. autosummary::
:toctree: generated/

trapz         -- Use trapezoidal rule to compute integral.
cumtrapz      -- Use trapezoidal rule to cumulatively compute integral.
simps         -- Use Simpson's rule to compute integral from samples.
romb          -- Use Romberg Integration to compute integral from
-- (2**k + 1) evenly-spaced samples.

.. seealso::

:mod:scipy.special for orthogonal polynomials (special) for Gaussian
quadrature roots and weights for other weighting factors and regions.

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

PACKAGE CONTENTS
_bvp
_dop
_ivp (package)
_ode
_odepack
_test_multivariate
_test_odeint_banded
lsoda
odepack
setup
tests (package)
vode

CLASSES
builtins.UserWarning(builtins.Warning)
builtins.object
scipy.integrate._ivp.base.DenseOutput
scipy.integrate._ivp.base.OdeSolver
scipy.integrate._ivp.bdf.BDF
scipy.integrate._ivp.lsoda.LSODA
scipy.integrate._ivp.common.OdeSolution
scipy.integrate._ode.ode
scipy.integrate._ode.complex_ode
scipy.integrate._ivp.rk.RungeKutta(scipy.integrate._ivp.base.OdeSolver)
scipy.integrate._ivp.rk.RK23
scipy.integrate._ivp.rk.RK45

class BDF(scipy.integrate._ivp.base.OdeSolver)
|  Implicit method based on Backward Differentiation Formulas.
|
|  This is a variable order method with the order varying automatically from
|  1 to 5. The general framework of the BDF algorithm is described in _.
|  This class implements a quasi-constant step size approach as explained
|  in _. The error estimation strategy for the constant step BDF is derived
|  in _. An accuracy enhancement using modified formulas (NDF) _ is also
|  implemented.
|
|  Can be applied in a complex domain.
|
|  Parameters
|  ----------
|  fun : callable
|      Right-hand side of the system. The calling signature is fun(t, y).
|      Here t is a scalar and there are two options for ndarray y.
|      It can either have shape (n,), then fun must return array_like with
|      shape (n,). Or alternatively it can have shape (n, k), then fun
|      must return array_like with shape (n, k), i.e. each column
|      corresponds to a single column in y. The choice between the two
|      options is determined by vectorized argument (see below). The
|      vectorized implementation allows faster approximation of the Jacobian
|      by finite differences.
|  t0 : float
|      Initial time.
|  y0 : array_like, shape (n,)
|      Initial state.
|  t_bound : float
|      Boundary time --- the integration won't continue beyond it. It also
|      determines the direction of the integration.
|  max_step : float, optional
|      Maximum allowed step size. Default is np.inf, i.e. the step is not
|      bounded and determined solely by the solver.
|  rtol, atol : float and array_like, optional
|      Relative and absolute tolerances. The solver keeps the local error
|      estimates less than atol + rtol * abs(y). Here rtol controls a
|      relative accuracy (number of correct digits). But if a component of y
|      is approximately below atol then the error only needs to fall within
|      the same atol threshold, and the number of correct digits is not
|      guaranteed. If components of y have different scales, it might be
|      beneficial to set different atol values for different components by
|      passing array_like with shape (n,) for atol. Default values are
|      1e-3 for rtol and 1e-6 for atol.
|  jac : {None, array_like, sparse_matrix, callable}, optional
|      Jacobian matrix of the right-hand side of the system with respect to
|      y, required only by 'Radau' and 'BDF' methods. The Jacobian matrix
|      has shape (n, n) and its element (i, j) is equal to d f_i / d y_j.
|      There are 3 ways to define the Jacobian:
|
|          * If array_like or sparse_matrix, then the Jacobian is assumed to
|            be constant.
|          * If callable, then the Jacobian is assumed to depend on both
|            t and y, and will be called as jac(t, y) as necessary. The
|            return value might be a sparse matrix.
|          * If None (default), then the Jacobian will be approximated by
|            finite differences.
|
|      It is generally recommended to provide the Jacobian rather than
|      relying on a finite difference approximation.
|  jac_sparsity : {None, array_like, sparse matrix}, optional
|      Defines a sparsity structure of the Jacobian matrix for a finite
|      difference approximation, its shape must be (n, n). If the Jacobian has
|      only few non-zero elements in *each* row, providing the sparsity
|      structure will greatly speed up the computations _. A zero
|      entry means that a corresponding element in the Jacobian is identically
|      zero. If None (default), the Jacobian is assumed to be dense.
|  vectorized : bool, optional
|      Whether fun is implemented in a vectorized fashion. Default is False.
|
|  Attributes
|  ----------
|  n : int
|      Number of equations.
|  status : string
|      Current status of the solver: 'running', 'finished' or 'failed'.
|  t_bound : float
|      Boundary time.
|  direction : float
|      Integration direction: +1 or -1.
|  t : float
|      Current time.
|  y : ndarray
|      Current state.
|  t_old : float
|      Previous time. None if no steps were made yet.
|  step_size : float
|      Size of the last successful step. None if no steps were made yet.
|  nfev : int
|      Number of the system's rhs evaluations.
|  njev : int
|      Number of the Jacobian evaluations.
|  nlu : int
|      Number of LU decompositions.
|
|  References
|  ----------
|  ..  G. D. Byrne, A. C. Hindmarsh, "A Polyalgorithm for the Numerical
|         Solution of Ordinary Differential Equations", ACM Transactions on
|         Mathematical Software, Vol. 1, No. 1, pp. 71-96, March 1975.
|  ..  L. F. Shampine, M. W. Reichelt, "THE MATLAB ODE SUITE", SIAM J. SCI.
|         COMPUTE., Vol. 18, No. 1, pp. 1-22, January 1997.
|  ..  E. Hairer, G. Wanner, "Solving Ordinary Differential Equations I:
|         Nonstiff Problems", Sec. III.2.
|  ..  A. Curtis, M. J. D. Powell, and J. Reid, "On the estimation of
|         sparse Jacobian matrices", Journal of the Institute of Mathematics
|         and its Applications, 13, pp. 117-120, 1974.
|
|  Method resolution order:
|      BDF
|      scipy.integrate._ivp.base.OdeSolver
|      builtins.object
|
|  Methods defined here:
|
|  __init__(self, fun, t0, y0, t_bound, max_step=inf, rtol=0.001, atol=1e-06, jac=None, jac_sparsity=None, vectorized=False, **extraneous)
|      Initialize self.  See help(type(self)) for accurate signature.
|
|  ----------------------------------------------------------------------
|  Methods inherited from scipy.integrate._ivp.base.OdeSolver:
|
|  dense_output(self)
|      Compute a local interpolant over the last successful step.
|
|      Returns
|      -------
|      sol : DenseOutput
|          Local interpolant over the last successful step.
|
|  step(self)
|      Perform one integration step.
|
|      Returns
|      -------
|      message : string or None
|          Report from the solver. Typically a reason for a failure if
|          self.status is 'failed' after the step was taken or None
|          otherwise.
|
|  ----------------------------------------------------------------------
|  Data descriptors inherited from scipy.integrate._ivp.base.OdeSolver:
|
|  __dict__
|      dictionary for instance variables (if defined)
|
|  __weakref__
|      list of weak references to the object (if defined)
|
|  step_size
|
|  ----------------------------------------------------------------------
|  Data and other attributes inherited from scipy.integrate._ivp.base.OdeSolver:
|
|  TOO_SMALL_STEP = 'Required step size is less than spacing between numb...

class DenseOutput(builtins.object)
|  Base class for local interpolant over step made by an ODE solver.
|
|  It interpolates between t_min and t_max (see Attributes below).
|  Evaluation outside this interval is not forbidden, but the accuracy is not
|  guaranteed.
|
|  Attributes
|  ----------
|  t_min, t_max : float
|      Time range of the interpolation.
|
|  Methods defined here:
|
|  __call__(self, t)
|      Evaluate the interpolant.
|
|      Parameters
|      ----------
|      t : float or array_like with shape (n_points,)
|          Points to evaluate the solution at.
|
|      Returns
|      -------
|      y : ndarray, shape (n,) or (n, n_points)
|          Computed values. Shape depends on whether t was a scalar or a
|          1-d array.
|
|  __init__(self, t_old, t)
|      Initialize self.  See help(type(self)) for accurate signature.
|
|  ----------------------------------------------------------------------
|  Data descriptors defined here:
|
|  __dict__
|      dictionary for instance variables (if defined)
|
|  __weakref__
|      list of weak references to the object (if defined)

class IntegrationWarning(builtins.UserWarning)
|  Warning on issues during integration.
|
|  Method resolution order:
|      IntegrationWarning
|      builtins.UserWarning
|      builtins.Warning
|      builtins.Exception
|      builtins.BaseException
|      builtins.object
|
|  Data descriptors defined here:
|
|  __weakref__
|      list of weak references to the object (if defined)
|
|  ----------------------------------------------------------------------
|  Methods inherited from builtins.UserWarning:
|
|  __init__(self, /, *args, **kwargs)
|      Initialize self.  See help(type(self)) for accurate signature.
|
|  __new__(*args, **kwargs) from builtins.type
|      Create and return a new object.  See help(type) for accurate signature.
|
|  ----------------------------------------------------------------------
|  Methods inherited from builtins.BaseException:
|
|  __delattr__(self, name, /)
|      Implement delattr(self, name).
|
|  __getattribute__(self, name, /)
|      Return getattr(self, name).
|
|  __reduce__(...)
|      helper for pickle
|
|  __repr__(self, /)
|      Return repr(self).
|
|  __setattr__(self, name, value, /)
|      Implement setattr(self, name, value).
|
|  __setstate__(...)
|
|  __str__(self, /)
|      Return str(self).
|
|  with_traceback(...)
|      Exception.with_traceback(tb) --
|      set self.__traceback__ to tb and return self.
|
|  ----------------------------------------------------------------------
|  Data descriptors inherited from builtins.BaseException:
|
|  __cause__
|      exception cause
|
|  __context__
|      exception context
|
|  __dict__
|
|  __suppress_context__
|
|  __traceback__
|
|  args

class LSODA(scipy.integrate._ivp.base.OdeSolver)
|  Adams/BDF method with automatic stiffness detection and switching.
|
|  This is a wrapper to the Fortran solver from ODEPACK _. It switches
|  automatically between the nonstiff Adams method and the stiff BDF method.
|  The method was originally detailed in _.
|
|  Parameters
|  ----------
|  fun : callable
|      Right-hand side of the system. The calling signature is fun(t, y).
|      Here t is a scalar and there are two options for ndarray y.
|      It can either have shape (n,), then fun must return array_like with
|      shape (n,). Or alternatively it can have shape (n, k), then fun
|      must return array_like with shape (n, k), i.e. each column
|      corresponds to a single column in y. The choice between the two
|      options is determined by vectorized argument (see below). The
|      vectorized implementation allows faster approximation of the Jacobian
|      by finite differences.
|  t0 : float
|      Initial time.
|  y0 : array_like, shape (n,)
|      Initial state.
|  t_bound : float
|      Boundary time --- the integration won't continue beyond it. It also
|      determines the direction of the integration.
|  first_step : float or None, optional
|      Initial step size. Default is None which means that the algorithm
|      should choose.
|  min_step : float, optional
|      Minimum allowed step size. Default is 0.0, i.e. the step is not
|      bounded and determined solely by the solver.
|  max_step : float, optional
|      Maximum allowed step size. Default is np.inf, i.e. the step is not
|      bounded and determined solely by the solver.
|  rtol, atol : float and array_like, optional
|      Relative and absolute tolerances. The solver keeps the local error
|      estimates less than atol + rtol * abs(y). Here rtol controls a
|      relative accuracy (number of correct digits). But if a component of y
|      is approximately below atol then the error only needs to fall within
|      the same atol threshold, and the number of correct digits is not
|      guaranteed. If components of y have different scales, it might be
|      beneficial to set different atol values for different components by
|      passing array_like with shape (n,) for atol. Default values are
|      1e-3 for rtol and 1e-6 for atol.
|  jac : None or callable, optional
|      Jacobian matrix of the right-hand side of the system with respect to
|      y. The Jacobian matrix has shape (n, n) and its element (i, j) is
|      equal to d f_i / d y_j. The function will be called as
|      jac(t, y). If None (default), then the Jacobian will be
|      approximated by finite differences. It is generally recommended to
|      provide the Jacobian rather than relying on a finite difference
|      approximation.
|  lband, uband : int or None, optional
|      Jacobian band width:
|      jac[i, j] != 0 only for i - lband <= j <= i + uband. Setting these
|      requires your jac routine to return the Jacobian in the packed format:
|      the returned array must have n columns and uband + lband + 1
|      rows in which Jacobian diagonals are written. Specifically
|      jac_packed[uband + i - j , j] = jac[i, j]. The same format is used
|      in scipy.linalg.solve_banded (check for an illustration).
|      These parameters can be also used with jac=None to reduce the
|      number of Jacobian elements estimated by finite differences.
|  vectorized : bool, optional
|      Whether fun is implemented in a vectorized fashion. A vectorized
|      implementation offers no advantages for this solver. Default is False.
|
|  Attributes
|  ----------
|  n : int
|      Number of equations.
|  status : string
|      Current status of the solver: 'running', 'finished' or 'failed'.
|  t_bound : float
|      Boundary time.
|  direction : float
|      Integration direction: +1 or -1.
|  t : float
|      Current time.
|  y : ndarray
|      Current state.
|  t_old : float
|      Previous time. None if no steps were made yet.
|  nfev : int
|      Number of the system's rhs evaluations.
|  njev : int
|      Number of the Jacobian evaluations.
|
|  References
|  ----------
|  ..  A. C. Hindmarsh, "ODEPACK, A Systematized Collection of ODE
|         Solvers," IMACS Transactions on Scientific Computation, Vol 1.,
|         pp. 55-64, 1983.
|  ..  L. Petzold, "Automatic selection of methods for solving stiff and
|         nonstiff systems of ordinary differential equations", SIAM Journal
|         on Scientific and Statistical Computing, Vol. 4, No. 1, pp. 136-148,
|         1983.
|
|  Method resolution order:
|      LSODA
|      scipy.integrate._ivp.base.OdeSolver
|      builtins.object
|
|  Methods defined here:
|
|  __init__(self, fun, t0, y0, t_bound, first_step=None, min_step=0.0, max_step=inf, rtol=0.001, atol=1e-06, jac=None, lband=None, uband=None, vectorized=False, **extraneous)
|      Initialize self.  See help(type(self)) for accurate signature.
|
|  ----------------------------------------------------------------------
|  Methods inherited from scipy.integrate._ivp.base.OdeSolver:
|
|  dense_output(self)
|      Compute a local interpolant over the last successful step.
|
|      Returns
|      -------
|      sol : DenseOutput
|          Local interpolant over the last successful step.
|
|  step(self)
|      Perform one integration step.
|
|      Returns
|      -------
|      message : string or None
|          Report from the solver. Typically a reason for a failure if
|          self.status is 'failed' after the step was taken or None
|          otherwise.
|
|  ----------------------------------------------------------------------
|  Data descriptors inherited from scipy.integrate._ivp.base.OdeSolver:
|
|  __dict__
|      dictionary for instance variables (if defined)
|
|  __weakref__
|      list of weak references to the object (if defined)
|
|  step_size
|
|  ----------------------------------------------------------------------
|  Data and other attributes inherited from scipy.integrate._ivp.base.OdeSolver:
|
|  TOO_SMALL_STEP = 'Required step size is less than spacing between numb...

class OdeSolution(builtins.object)
|  Continuous ODE solution.
|
|  It is organized as a collection of DenseOutput objects which represent
|  local interpolants. It provides an algorithm to select a right interpolant
|  for each given point.
|
|  The interpolants cover the range between t_min and t_max (see
|  Attributes below). Evaluation outside this interval is not forbidden, but
|  the accuracy is not guaranteed.
|
|  When evaluating at a breakpoint (one of the values in ts) a segment with
|  the lower index is selected.
|
|  Parameters
|  ----------
|  ts : array_like, shape (n_segments + 1,)
|      Time instants between which local interpolants are defined. Must
|      be strictly increasing or decreasing (zero segment with two points is
|      also allowed).
|  interpolants : list of DenseOutput with n_segments elements
|      Local interpolants. An i-th interpolant is assumed to be defined
|      between ts[i] and ts[i + 1].
|
|  Attributes
|  ----------
|  t_min, t_max : float
|      Time range of the interpolation.
|
|  Methods defined here:
|
|  __call__(self, t)
|      Evaluate the solution.
|
|      Parameters
|      ----------
|      t : float or array_like with shape (n_points,)
|          Points to evaluate at.
|
|      Returns
|      -------
|      y : ndarray, shape (n_states,) or (n_states, n_points)
|          Computed values. Shape depends on whether t is a scalar or a
|          1-d array.
|
|  __init__(self, ts, interpolants)
|      Initialize self.  See help(type(self)) for accurate signature.
|
|  ----------------------------------------------------------------------
|  Data descriptors defined here:
|
|  __dict__
|      dictionary for instance variables (if defined)
|
|  __weakref__
|      list of weak references to the object (if defined)

class OdeSolver(builtins.object)
|  Base class for ODE solvers.
|
|  In order to implement a new solver you need to follow the guidelines:
|
|      1. A constructor must accept parameters presented in the base class
|         (listed below) along with any other parameters specific to a solver.
|      2. A constructor must accept arbitrary extraneous arguments
|         **extraneous, but warn that these arguments are irrelevant
|         using common.warn_extraneous function. Do not pass these
|         arguments to the base class.
|      3. A solver must implement a private method _step_impl(self) which
|         propagates a solver one step further. It must return tuple
|         (success, message), where success is a boolean indicating
|         whether a step was successful, and message is a string
|         containing description of a failure if a step failed or None
|         otherwise.
|      4. A solver must implement a private method _dense_output_impl(self)
|         which returns a DenseOutput object covering the last successful
|         step.
|      5. A solver must have attributes listed below in Attributes section.
|         Note that t_old and step_size are updated automatically.
|      6. Use fun(self, t, y) method for the system rhs evaluation, this
|         way the number of function evaluations (nfev) will be tracked
|         automatically.
|      7. For convenience a base class provides fun_single(self, t, y) and
|         fun_vectorized(self, t, y) for evaluating the rhs in
|         non-vectorized and vectorized fashions respectively (regardless of
|         how fun from the constructor is implemented). These calls don't
|         increment nfev.
|      8. If a solver uses a Jacobian matrix and LU decompositions, it should
|         track the number of Jacobian evaluations (njev) and the number of
|         LU decompositions (nlu).
|      9. By convention the function evaluations used to compute a finite
|         difference approximation of the Jacobian should not be counted in
|         nfev, thus use fun_single(self, t, y) or
|         fun_vectorized(self, t, y) when computing a finite difference
|         approximation of the Jacobian.
|
|  Parameters
|  ----------
|  fun : callable
|      Right-hand side of the system. The calling signature is fun(t, y).
|      Here t is a scalar and there are two options for ndarray y.
|      It can either have shape (n,), then fun must return array_like with
|      shape (n,). Or alternatively it can have shape (n, n_points), then
|      fun must return array_like with shape (n, n_points) (each column
|      corresponds to a single column in y). The choice between the two
|      options is determined by vectorized argument (see below).
|  t0 : float
|      Initial time.
|  y0 : array_like, shape (n,)
|      Initial state.
|  t_bound : float
|      Boundary time --- the integration won't continue beyond it. It also
|      determines the direction of the integration.
|  vectorized : bool
|      Whether fun is implemented in a vectorized fashion.
|  support_complex : bool, optional
|      Whether integration in a complex domain should be supported.
|      Generally determined by a derived solver class capabilities.
|      Default is False.
|
|  Attributes
|  ----------
|  n : int
|      Number of equations.
|  status : string
|      Current status of the solver: 'running', 'finished' or 'failed'.
|  t_bound : float
|      Boundary time.
|  direction : float
|      Integration direction: +1 or -1.
|  t : float
|      Current time.
|  y : ndarray
|      Current state.
|  t_old : float
|      Previous time. None if no steps were made yet.
|  step_size : float
|      Size of the last successful step. None if no steps were made yet.
|  nfev : int
|      Number of the system's rhs evaluations.
|  njev : int
|      Number of the Jacobian evaluations.
|  nlu : int
|      Number of LU decompositions.
|
|  Methods defined here:
|
|  __init__(self, fun, t0, y0, t_bound, vectorized, support_complex=False)
|      Initialize self.  See help(type(self)) for accurate signature.
|
|  dense_output(self)
|      Compute a local interpolant over the last successful step.
|
|      Returns
|      -------
|      sol : DenseOutput
|          Local interpolant over the last successful step.
|
|  step(self)
|      Perform one integration step.
|
|      Returns
|      -------
|      message : string or None
|          Report from the solver. Typically a reason for a failure if
|          self.status is 'failed' after the step was taken or None
|          otherwise.
|
|  ----------------------------------------------------------------------
|  Data descriptors defined here:
|
|  __dict__
|      dictionary for instance variables (if defined)
|
|  __weakref__
|      list of weak references to the object (if defined)
|
|  step_size
|
|  ----------------------------------------------------------------------
|  Data and other attributes defined here:
|
|  TOO_SMALL_STEP = 'Required step size is less than spacing between numb...

class RK23(RungeKutta)
|  Explicit Runge-Kutta method of order 3(2).
|
|  The Bogacki-Shamping pair of formulas is used _. The error is controlled
|  assuming 2nd order accuracy, but steps are taken using a 3rd oder accurate
|  formula (local extrapolation is done). A cubic Hermit polynomial is used
|  for the dense output.
|
|  Can be applied in a complex domain.
|
|  Parameters
|  ----------
|  fun : callable
|      Right-hand side of the system. The calling signature is fun(t, y).
|      Here t is a scalar and there are two options for ndarray y.
|      It can either have shape (n,), then fun must return array_like with
|      shape (n,). Or alternatively it can have shape (n, k), then fun
|      must return array_like with shape (n, k), i.e. each column
|      corresponds to a single column in y. The choice between the two
|      options is determined by vectorized argument (see below). The
|      vectorized implementation allows faster approximation of the Jacobian
|      by finite differences.
|  t0 : float
|      Initial time.
|  y0 : array_like, shape (n,)
|      Initial state.
|  t_bound : float
|      Boundary time --- the integration won't continue beyond it. It also
|      determines the direction of the integration.
|  max_step : float, optional
|      Maximum allowed step size. Default is np.inf, i.e. the step is not
|      bounded and determined solely by the solver.
|  rtol, atol : float and array_like, optional
|      Relative and absolute tolerances. The solver keeps the local error
|      estimates less than atol + rtol * abs(y). Here rtol controls a
|      relative accuracy (number of correct digits). But if a component of y
|      is approximately below atol then the error only needs to fall within
|      the same atol threshold, and the number of correct digits is not
|      guaranteed. If components of y have different scales, it might be
|      beneficial to set different atol values for different components by
|      passing array_like with shape (n,) for atol. Default values are
|      1e-3 for rtol and 1e-6 for atol.
|  vectorized : bool, optional
|      Whether fun is implemented in a vectorized fashion. Default is False.
|
|  Attributes
|  ----------
|  n : int
|      Number of equations.
|  status : string
|      Current status of the solver: 'running', 'finished' or 'failed'.
|  t_bound : float
|      Boundary time.
|  direction : float
|      Integration direction: +1 or -1.
|  t : float
|      Current time.
|  y : ndarray
|      Current state.
|  t_old : float
|      Previous time. None if no steps were made yet.
|  step_size : float
|      Size of the last successful step. None if no steps were made yet.
|  nfev : int
|      Number of the system's rhs evaluations.
|  njev : int
|      Number of the Jacobian evaluations.
|  nlu : int
|      Number of LU decompositions.
|
|  References
|  ----------
|  ..  P. Bogacki, L.F. Shampine, "A 3(2) Pair of Runge-Kutta Formulas",
|         Appl. Math. Lett. Vol. 2, No. 4. pp. 321-325, 1989.
|
|  Method resolution order:
|      RK23
|      RungeKutta
|      scipy.integrate._ivp.base.OdeSolver
|      builtins.object
|
|  Data and other attributes defined here:
|
|  A = [array([0.5]), array([0.  , 0.75])]
|
|  B = array([0.22222222, 0.33333333, 0.44444444])
|
|  C = array([0.5 , 0.75])
|
|  E = array([ 0.06944444, -0.08333333, -0.11111111,  0.125     ])
|
|  P = array([[ 1.        , -1.33333333,  0.55555556],
|  ...
|         [ 0.   ...
|
|  n_stages = 3
|
|  order = 2
|
|  ----------------------------------------------------------------------
|  Methods inherited from RungeKutta:
|
|  __init__(self, fun, t0, y0, t_bound, max_step=inf, rtol=0.001, atol=1e-06, vectorized=False, **extraneous)
|      Initialize self.  See help(type(self)) for accurate signature.
|
|  ----------------------------------------------------------------------
|  Methods inherited from scipy.integrate._ivp.base.OdeSolver:
|
|  dense_output(self)
|      Compute a local interpolant over the last successful step.
|
|      Returns
|      -------
|      sol : DenseOutput
|          Local interpolant over the last successful step.
|
|  step(self)
|      Perform one integration step.
|
|      Returns
|      -------
|      message : string or None
|          Report from the solver. Typically a reason for a failure if
|          self.status is 'failed' after the step was taken or None
|          otherwise.
|
|  ----------------------------------------------------------------------
|  Data descriptors inherited from scipy.integrate._ivp.base.OdeSolver:
|
|  __dict__
|      dictionary for instance variables (if defined)
|
|  __weakref__
|      list of weak references to the object (if defined)
|
|  step_size
|
|  ----------------------------------------------------------------------
|  Data and other attributes inherited from scipy.integrate._ivp.base.OdeSolver:
|
|  TOO_SMALL_STEP = 'Required step size is less than spacing between numb...

class RK45(RungeKutta)
|  Explicit Runge-Kutta method of order 5(4).
|
|  The Dormand-Prince pair of formulas is used _. The error is controlled
|  assuming 4th order accuracy, but steps are taken using a 5th
|  oder accurate formula (local extrapolation is done). A quartic
|  interpolation polynomial is used for the dense output _.
|
|  Can be applied in a complex domain.
|
|  Parameters
|  ----------
|  fun : callable
|      Right-hand side of the system. The calling signature is fun(t, y).
|      Here t is a scalar and there are two options for ndarray y.
|      It can either have shape (n,), then fun must return array_like with
|      shape (n,). Or alternatively it can have shape (n, k), then fun
|      must return array_like with shape (n, k), i.e. each column
|      corresponds to a single column in y. The choice between the two
|      options is determined by vectorized argument (see below). The
|      vectorized implementation allows faster approximation of the Jacobian
|      by finite differences.
|  t0 : float
|      Initial value of the independent variable.
|  y0 : array_like, shape (n,)
|      Initial values of the dependent variable.
|  t_bound : float
|      Boundary time --- the integration won't continue beyond it. It also
|      determines the direction of the integration.
|  max_step : float, optional
|      Maximum allowed step size. Default is np.inf, i.e. the step is not
|      bounded and determined solely by the solver.
|  rtol, atol : float and array_like, optional
|      Relative and absolute tolerances. The solver keeps the local error
|      estimates less than atol + rtol * abs(y). Here rtol controls a
|      relative accuracy (number of correct digits). But if a component of y
|      is approximately below atol then the error only needs to fall within
|      the same atol threshold, and the number of correct digits is not
|      guaranteed. If components of y have different scales, it might be
|      beneficial to set different atol values for different components by
|      passing array_like with shape (n,) for atol. Default values are
|      1e-3 for rtol and 1e-6 for atol.
|  vectorized : bool, optional
|      Whether fun is implemented in a vectorized fashion. Default is False.
|
|  Attributes
|  ----------
|  n : int
|      Number of equations.
|  status : string
|      Current status of the solver: 'running', 'finished' or 'failed'.
|  t_bound : float
|      Boundary time.
|  direction : float
|      Integration direction: +1 or -1.
|  t : float
|      Current time.
|  y : ndarray
|      Current state.
|  t_old : float
|      Previous time. None if no steps were made yet.
|  step_size : float
|      Size of the last successful step. None if no steps were made yet.
|  nfev : int
|      Number of the system's rhs evaluations.
|  njev : int
|      Number of the Jacobian evaluations.
|  nlu : int
|      Number of LU decompositions.
|
|  References
|  ----------
|  ..  J. R. Dormand, P. J. Prince, "A family of embedded Runge-Kutta
|         formulae", Journal of Computational and Applied Mathematics, Vol. 6,
|         No. 1, pp. 19-26, 1980.
|  ..  L. W. Shampine, "Some Practical Runge-Kutta Formulas", Mathematics
|         of Computation,, Vol. 46, No. 173, pp. 135-150, 1986.
|
|  Method resolution order:
|      RK45
|      RungeKutta
|      scipy.integrate._ivp.base.OdeSolver
|      builtins.object
|
|  Data and other attributes defined here:
|
|  A = [array([0.2]), array([0.075, 0.225]), array([ 0.97777778, -3.73333...
|
|  B = array([ 0.09114583,  0.        ,  0.4492363 ,  0.65104167, -0.3223...
|
|  C = array([0.2       , 0.3       , 0.8       , 0.88888889, 1.        ]...
|
|  E = array([-0.00123264,  0.        ,  0.00425277, -0...7,  0.0508638 ,...
|
|  P = array([[ 1.        , -2.85358007,  3.07174346, -...       ,  1.382...
|
|  n_stages = 6
|
|  order = 4
|
|  ----------------------------------------------------------------------
|  Methods inherited from RungeKutta:
|
|  __init__(self, fun, t0, y0, t_bound, max_step=inf, rtol=0.001, atol=1e-06, vectorized=False, **extraneous)
|      Initialize self.  See help(type(self)) for accurate signature.
|
|  ----------------------------------------------------------------------
|  Methods inherited from scipy.integrate._ivp.base.OdeSolver:
|
|  dense_output(self)
|      Compute a local interpolant over the last successful step.
|
|      Returns
|      -------
|      sol : DenseOutput
|          Local interpolant over the last successful step.
|
|  step(self)
|      Perform one integration step.
|
|      Returns
|      -------
|      message : string or None
|          Report from the solver. Typically a reason for a failure if
|          self.status is 'failed' after the step was taken or None
|          otherwise.
|
|  ----------------------------------------------------------------------
|  Data descriptors inherited from scipy.integrate._ivp.base.OdeSolver:
|
|  __dict__
|      dictionary for instance variables (if defined)
|
|  __weakref__
|      list of weak references to the object (if defined)
|
|  step_size
|
|  ----------------------------------------------------------------------
|  Data and other attributes inherited from scipy.integrate._ivp.base.OdeSolver:
|
|  TOO_SMALL_STEP = 'Required step size is less than spacing between numb...

|  Implicit Runge-Kutta method of Radau IIA family of order 5.
|
|  Implementation follows _. The error is controlled for a 3rd order
|  accurate embedded formula. A cubic polynomial which satisfies the
|  collocation conditions is used for the dense output.
|
|  Parameters
|  ----------
|  fun : callable
|      Right-hand side of the system. The calling signature is fun(t, y).
|      Here t is a scalar and there are two options for ndarray y.
|      It can either have shape (n,), then fun must return array_like with
|      shape (n,). Or alternatively it can have shape (n, k), then fun
|      must return array_like with shape (n, k), i.e. each column
|      corresponds to a single column in y. The choice between the two
|      options is determined by vectorized argument (see below). The
|      vectorized implementation allows faster approximation of the Jacobian
|      by finite differences.
|  t0 : float
|      Initial time.
|  y0 : array_like, shape (n,)
|      Initial state.
|  t_bound : float
|      Boundary time --- the integration won't continue beyond it. It also
|      determines the direction of the integration.
|  max_step : float, optional
|      Maximum allowed step size. Default is np.inf, i.e. the step is not
|      bounded and determined solely by the solver.
|  rtol, atol : float and array_like, optional
|      Relative and absolute tolerances. The solver keeps the local error
|      estimates less than atol + rtol * abs(y). Here rtol controls a
|      relative accuracy (number of correct digits). But if a component of y
|      is approximately below atol then the error only needs to fall within
|      the same atol threshold, and the number of correct digits is not
|      guaranteed. If components of y have different scales, it might be
|      beneficial to set different atol values for different components by
|      passing array_like with shape (n,) for atol. Default values are
|      1e-3 for rtol and 1e-6 for atol.
|  jac : {None, array_like, sparse_matrix, callable}, optional
|      Jacobian matrix of the right-hand side of the system with respect to
|      y, required only by 'Radau' and 'BDF' methods. The Jacobian matrix
|      has shape (n, n) and its element (i, j) is equal to d f_i / d y_j.
|      There are 3 ways to define the Jacobian:
|
|          * If array_like or sparse_matrix, then the Jacobian is assumed to
|            be constant.
|          * If callable, then the Jacobian is assumed to depend on both
|            t and y, and will be called as jac(t, y) as necessary. The
|            return value might be a sparse matrix.
|          * If None (default), then the Jacobian will be approximated by
|            finite differences.
|
|      It is generally recommended to provide the Jacobian rather than
|      relying on a finite difference approximation.
|  jac_sparsity : {None, array_like, sparse matrix}, optional
|      Defines a sparsity structure of the Jacobian matrix for a finite
|      difference approximation, its shape must be (n, n). If the Jacobian has
|      only few non-zero elements in *each* row, providing the sparsity
|      structure will greatly speed up the computations _. A zero
|      entry means that a corresponding element in the Jacobian is identically
|      zero. If None (default), the Jacobian is assumed to be dense.
|  vectorized : bool, optional
|      Whether fun is implemented in a vectorized fashion. Default is False.
|
|  Attributes
|  ----------
|  n : int
|      Number of equations.
|  status : string
|      Current status of the solver: 'running', 'finished' or 'failed'.
|  t_bound : float
|      Boundary time.
|  direction : float
|      Integration direction: +1 or -1.
|  t : float
|      Current time.
|  y : ndarray
|      Current state.
|  t_old : float
|      Previous time. None if no steps were made yet.
|  step_size : float
|      Size of the last successful step. None if no steps were made yet.
|  nfev : int
|      Number of the system's rhs evaluations.
|  njev : int
|      Number of the Jacobian evaluations.
|  nlu : int
|      Number of LU decompositions.
|
|  References
|  ----------
|  ..  E. Hairer, G. Wanner, "Solving Ordinary Differential Equations II:
|         Stiff and Differential-Algebraic Problems", Sec. IV.8.
|  ..  A. Curtis, M. J. D. Powell, and J. Reid, "On the estimation of
|         sparse Jacobian matrices", Journal of the Institute of Mathematics
|         and its Applications, 13, pp. 117-120, 1974.
|
|  Method resolution order:
|      scipy.integrate._ivp.base.OdeSolver
|      builtins.object
|
|  Methods defined here:
|
|  __init__(self, fun, t0, y0, t_bound, max_step=inf, rtol=0.001, atol=1e-06, jac=None, jac_sparsity=None, vectorized=False, **extraneous)
|      Initialize self.  See help(type(self)) for accurate signature.
|
|  ----------------------------------------------------------------------
|  Methods inherited from scipy.integrate._ivp.base.OdeSolver:
|
|  dense_output(self)
|      Compute a local interpolant over the last successful step.
|
|      Returns
|      -------
|      sol : DenseOutput
|          Local interpolant over the last successful step.
|
|  step(self)
|      Perform one integration step.
|
|      Returns
|      -------
|      message : string or None
|          Report from the solver. Typically a reason for a failure if
|          self.status is 'failed' after the step was taken or None
|          otherwise.
|
|  ----------------------------------------------------------------------
|  Data descriptors inherited from scipy.integrate._ivp.base.OdeSolver:
|
|  __dict__
|      dictionary for instance variables (if defined)
|
|  __weakref__
|      list of weak references to the object (if defined)
|
|  step_size
|
|  ----------------------------------------------------------------------
|  Data and other attributes inherited from scipy.integrate._ivp.base.OdeSolver:
|
|  TOO_SMALL_STEP = 'Required step size is less than spacing between numb...

class complex_ode(ode)
|  A wrapper of ode for complex systems.
|
|  This functions similarly as ode, but re-maps a complex-valued
|  equation system to a real-valued one before using the integrators.
|
|  Parameters
|  ----------
|  f : callable f(t, y, *f_args)
|      Rhs of the equation. t is a scalar, y.shape == (n,).
|      f_args is set by calling set_f_params(*args).
|  jac : callable jac(t, y, *jac_args)
|      Jacobian of the rhs, jac[i,j] = d f[i] / d y[j].
|      jac_args is set by calling set_f_params(*args).
|
|  Attributes
|  ----------
|  t : float
|      Current time.
|  y : ndarray
|      Current variable values.
|
|  Examples
|  --------
|  For usage examples, see ode.
|
|  Method resolution order:
|      complex_ode
|      ode
|      builtins.object
|
|  Methods defined here:
|
|  __init__(self, f, jac=None)
|      Initialize self.  See help(type(self)) for accurate signature.
|
|  integrate(self, t, step=False, relax=False)
|      Find y=y(t), set y as an initial condition, and return y.
|
|      Parameters
|      ----------
|      t : float
|          The endpoint of the integration step.
|      step : bool
|          If True, and if the integrator supports the step method,
|          then perform a single integration step and return.
|          This parameter is provided in order to expose internals of
|          the implementation, and should not be changed from its default
|          value in most cases.
|      relax : bool
|          If True and if the integrator supports the run_relax method,
|          then integrate until t_1 >= t and return. relax is not
|          referenced if step=True.
|          This parameter is provided in order to expose internals of
|          the implementation, and should not be changed from its default
|          value in most cases.
|
|      Returns
|      -------
|      y : float
|          The integrated value at t
|
|  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:
|
|  get_return_code(self)
|      Extracts the return code for the integration to enable better control
|      if the integration fails.
|
|  set_f_params(self, *args)
|      Set extra parameters for user-supplied function f.
|
|  set_jac_params(self, *args)
|      Set extra parameters for user-supplied function jac.
|
|  successful(self)
|      Check if integration was successful.
|
|  ----------------------------------------------------------------------
|  Data descriptors inherited from ode:
|
|  __dict__
|      dictionary for instance variables (if defined)
|
|  __weakref__
|      list of weak references to the object (if defined)

class ode(builtins.object)
|  A generic interface class to numeric integrators.
|
|  Solve an equation system :math:y'(t) = f(t,y) with (optional) jac = df/dy.
|
|  *Note*: The first two arguments of f(t, y, ...) are in the
|  opposite order of the arguments in the system definition function used
|  by scipy.integrate.odeint.
|
|  Parameters
|  ----------
|  f : callable f(t, y, *f_args)
|      Right-hand side of the differential equation. t is a scalar,
|      y.shape == (n,).
|      f_args is set by calling set_f_params(*args).
|      f should return a scalar, array or list (not a tuple).
|  jac : callable jac(t, y, *jac_args), optional
|      Jacobian of the right-hand side, jac[i,j] = d f[i] / d y[j].
|      jac_args is set by calling set_jac_params(*args).
|
|  Attributes
|  ----------
|  t : float
|      Current time.
|  y : ndarray
|      Current variable values.
|
|  --------
|  odeint : an integrator with a simpler interface based on lsoda from ODEPACK
|  quad : for finding the area under a curve
|
|  Notes
|  -----
|  Available integrators are listed below. They can be selected using
|  the set_integrator method.
|
|  "vode"
|
|      Real-valued Variable-coefficient Ordinary Differential Equation
|      solver, with fixed-leading-coefficient implementation. It provides
|      implicit Adams method (for non-stiff problems) and a method based on
|      backward differentiation formulas (BDF) (for stiff problems).
|
|      Source: http://www.netlib.org/ode/vode.f
|
|      .. warning::
|
|         This integrator is not re-entrant. You cannot have two ode
|         instances using the "vode" integrator at the same time.
|
|      This integrator accepts the following parameters in set_integrator
|      method of the ode class:
|
|      - atol : float or sequence
|        absolute tolerance for solution
|      - rtol : float or sequence
|        relative tolerance for solution
|      - lband : None or int
|      - uband : None or int
|        Jacobian band width, jac[i,j] != 0 for i-lband <= j <= i+uband.
|        Setting these requires your jac routine to return the jacobian
|        in packed format, jac_packed[i-j+uband, j] = jac[i,j]. The
|        dimension of the matrix must be (lband+uband+1, len(y)).
|      - method: 'adams' or 'bdf'
|        Which solver to use, Adams (non-stiff) or BDF (stiff)
|      - with_jacobian : bool
|        This option is only considered when the user has not supplied a
|        Jacobian function and has not indicated (by setting either band)
|        that the Jacobian is banded.  In this case, with_jacobian specifies
|        whether the iteration method of the ODE solver's correction step is
|        chord iteration with an internally generated full Jacobian or
|        functional iteration with no Jacobian.
|      - nsteps : int
|        Maximum number of (internally defined) steps allowed during one
|        call to the solver.
|      - first_step : float
|      - min_step : float
|      - max_step : float
|        Limits for the step sizes used by the integrator.
|      - order : int
|        Maximum order used by the integrator,
|        order <= 12 for Adams, <= 5 for BDF.
|
|  "zvode"
|
|      Complex-valued Variable-coefficient Ordinary Differential Equation
|      solver, with fixed-leading-coefficient implementation.  It provides
|      implicit Adams method (for non-stiff problems) and a method based on
|      backward differentiation formulas (BDF) (for stiff problems).
|
|      Source: http://www.netlib.org/ode/zvode.f
|
|      .. warning::
|
|         This integrator is not re-entrant. You cannot have two ode
|         instances using the "zvode" integrator at the same time.
|
|      This integrator accepts the same parameters in set_integrator
|      as the "vode" solver.
|
|      .. note::
|
|          When using ZVODE for a stiff system, it should only be used for
|          the case in which the function f is analytic, that is, when each f(i)
|          is an analytic function of each y(j).  Analyticity means that the
|          partial derivative df(i)/dy(j) is a unique complex number, and this
|          fact is critical in the way ZVODE solves the dense or banded linear
|          systems that arise in the stiff case.  For a complex stiff ODE system
|          in which f is not analytic, ZVODE is likely to have convergence
|          failures, and for this problem one should instead use DVODE on the
|          equivalent real system (in the real and imaginary parts of y).
|
|  "lsoda"
|
|      Real-valued Variable-coefficient Ordinary Differential Equation
|      solver, with fixed-leading-coefficient implementation. It provides
|      automatic method switching between implicit Adams method (for non-stiff
|      problems) and a method based on backward differentiation formulas (BDF)
|      (for stiff problems).
|
|      Source: http://www.netlib.org/odepack
|
|      .. warning::
|
|         This integrator is not re-entrant. You cannot have two ode
|         instances using the "lsoda" integrator at the same time.
|
|      This integrator accepts the following parameters in set_integrator
|      method of the ode class:
|
|      - atol : float or sequence
|        absolute tolerance for solution
|      - rtol : float or sequence
|        relative tolerance for solution
|      - lband : None or int
|      - uband : None or int
|        Jacobian band width, jac[i,j] != 0 for i-lband <= j <= i+uband.
|        Setting these requires your jac routine to return the jacobian
|        in packed format, jac_packed[i-j+uband, j] = jac[i,j].
|      - with_jacobian : bool
|        *Not used.*
|      - nsteps : int
|        Maximum number of (internally defined) steps allowed during one
|        call to the solver.
|      - first_step : float
|      - min_step : float
|      - max_step : float
|        Limits for the step sizes used by the integrator.
|      - max_order_ns : int
|        Maximum order used in the nonstiff case (default 12).
|      - max_order_s : int
|        Maximum order used in the stiff case (default 5).
|      - max_hnil : int
|        Maximum number of messages reporting too small step size (t + h = t)
|        (default 0)
|      - ixpr : int
|        Whether to generate extra printing at method switches (default False).
|
|  "dopri5"
|
|      This is an explicit runge-kutta method of order (4)5 due to Dormand &
|      Prince (with stepsize control and dense output).
|
|      Authors:
|
|          E. Hairer and G. Wanner
|          Universite de Geneve, Dept. de Mathematiques
|          CH-1211 Geneve 24, Switzerland
|          e-mail:  [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 + y, -arg1*y**2]
|  >>> def jac(t, y, arg1):
|  ...     return [[1j*arg1, 1], [0, -arg1*2*y]]
|
|  The integration:
|
|  >>> r = ode(f, jac).set_integrator('zvode', method='bdf')
|  >>> r.set_initial_value(y0, t0).set_f_params(2.0).set_jac_params(2.0)
|  >>> t1 = 10
|  >>> dt = 1
|  >>> while r.successful() and r.t < t1:
|  ...     print(r.t+dt, r.integrate(r.t+dt))
|  1 [-0.71038232+0.23749653j  0.40000271+0.j        ]
|  2.0 [ 0.19098503-0.52359246j  0.22222356+0.j        ]
|  3.0 [ 0.47153208+0.52701229j  0.15384681+0.j        ]
|  4.0 [-0.61905937+0.30726255j  0.11764744+0.j        ]
|  5.0 [ 0.02340997-0.61418799j  0.09523835+0.j        ]
|  6.0 [ 0.58643071+0.339819j  0.08000018+0.j      ]
|  7.0 [-0.52070105+0.44525141j  0.06896565+0.j        ]
|  8.0 [-0.15986733-0.61234476j  0.06060616+0.j        ]
|  9.0 [ 0.64850462+0.15048982j  0.05405414+0.j        ]
|  10.0 [-0.38404699+0.56382299j  0.04878055+0.j        ]
|
|  References
|  ----------
|  .. [HNW93] E. Hairer, S.P. Norsett and G. Wanner, Solving Ordinary
|      Differential Equations i. Nonstiff Problems. 2nd edition.
|      Springer Series in Computational Mathematics,
|      Springer-Verlag (1993)
|
|  Methods defined here:
|
|  __init__(self, f, jac=None)
|      Initialize self.  See help(type(self)) for accurate signature.
|
|  get_return_code(self)
|      Extracts the return code for the integration to enable better control
|      if the integration fails.
|
|  integrate(self, t, step=False, relax=False)
|      Find y=y(t), set y as an initial condition, and return y.
|
|      Parameters
|      ----------
|      t : float
|          The endpoint of the integration step.
|      step : bool
|          If True, and if the integrator supports the step method,
|          then perform a single integration step and return.
|          This parameter is provided in order to expose internals of
|          the implementation, and should not be changed from its default
|          value in most cases.
|      relax : bool
|          If True and if the integrator supports the run_relax method,
|          then integrate until t_1 >= t and return. relax is not
|          referenced if step=True.
|          This parameter is provided in order to expose internals of
|          the implementation, and should not be changed from its default
|          value in most cases.
|
|      Returns
|      -------
|      y : float
|          The integrated value at t
|
|  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 : float, 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 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.

--------
numpy.cumsum, numpy.cumprod
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.5 * x**2, 'b-')
>>> plt.show()

dblquad(func, a, b, gfun, hfun, args=(), epsabs=1.49e-08, epsrel=1.49e-08)
Compute a double integral.

Return the double (definite) integral of func(y, x) from x = a..b
and y = gfun(x)..hfun(x).

Parameters
----------
func : callable
A Python function or method of at least two variables: y must be the
first argument and x the second argument.
a, b : float
The limits of integration in x: a < b
gfun : callable
The lower boundary curve in y which is a function taking a single
floating point argument (x) and returning a floating point result: a
lambda function can be useful here.
hfun : callable
The upper boundary curve in y (same requirements as gfun).
args : sequence, optional
Extra arguments to pass to func.
epsabs : float, optional
Absolute tolerance passed directly to the inner 1-D quadrature
integration. Default is 1.49e-8.
epsrel : float, optional
Relative tolerance of the inner 1-D integrals. Default is 1.49e-8.

Returns
-------
y : float
The resultant integral.
abserr : float
An estimate of the error.

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

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).
If integrating a vector-valued function, the returned array must have
shape (..., len(x)).
a : float
Lower limit of integration.
b : float
Upper limit of integration.
args : tuple, optional
Extra arguments to pass to function, if any.
n : int, optional
Order of quadrature integration. Default is 5.

Returns
-------
val : float
Gaussian quadrature approximation to the integral
none : None
Statically returned value of None

--------
romb : integrators for sampled data
simps : integrators for sampled data
cumtrapz : cumulative integration for sampled data
ode : ODE integrator
odeint : ODE integrator

newton_cotes(rn, equal=0)
Return weights and error coefficient for Newton-Cotes integration.

Suppose we have (N+1) samples of f at the positions
x_0, x_1, ..., x_N.  Then an N-point Newton-Cotes formula for the
integral between x_0 and x_N is:

:math:\int_{x_0}^{x_N} f(x)dx = \Delta x \sum_{i=0}^{N} a_i f(x_i)
+ B_N (\Delta x)^{N+2} f^{N+1} (\xi)

where :math:\xi \in [x_0,x_N]
and :math:\Delta x = \frac{x_N-x_0}{N} is the average samples spacing.

If the samples are equally-spaced and N is even, then the error
term is :math:B_N (\Delta x)^{N+3} f^{N+2}(\xi).

Parameters
----------
rn : int
The integer order for equally-spaced data or the relative positions of
the samples with the first sample at 0 and the last at N, where N+1 is
the length of rn.  N is the order of the Newton-Cotes integration.
equal : int, optional
Set to 1 to enforce equally spaced data.

Returns
-------
an : ndarray
1-D array of weights to apply to the function at the provided sample
positions.
B : float
Error coefficient.

Notes
-----
Normally, the Newton-Cotes rules are used on smaller integration
regions and a composite rule is used to return the total integral.

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, scipy.LowLevelCallable}
The function to be integrated. Has arguments of x0, ... xn,
t0, tm, where integration is carried out over x0, ... xn, which
must be floats.  Function signature should be
func(x0, x1, ..., xn, t0, t1, ..., tm).  Integration is carried out
in order.  That is, integration over x0 is the innermost integral,
and xn is the outermost.

If the user desires improved integration performance, then f may
be a scipy.LowLevelCallable with one of the signatures::

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

where n is the number of extra parameters and args is an array
of doubles of the additional parameters, the xx array contains the
coordinates. The user_data is the data contained in the
scipy.LowLevelCallable.
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 corresponds to
integration over x0, and so on.  If an element of ranges is a callable,
then it will be called with all of the integration arguments available,
as well as any parametric arguments. e.g. if
func = f(x0, x1, x2, t0, t1), then ranges may be defined as
either (a, b) or else as (a, b) = range0(x1, x2, t0, t1).
args : iterable object, optional
Additional arguments t0, ..., tn, required by func, ranges, and
opts.
opts : iterable object or dict, optional
Options to be passed to quad.  May be empty, a dict, or
a sequence of dicts or functions that return a dict.  If empty, the
default options from scipy.integrate.quad are used.  If a dict, the same
options are used for all levels of integraion.  If a sequence, then each
element of the sequence corresponds to a particular integration. e.g.
opts corresponds to integration over x0, and so on. If a callable,
the signature must be the same as for ranges. The available
options together with their default values are:

- epsabs = 1.49e-08
- epsrel = 1.49e-08
- limit  = 50
- points = None
- weight = None
- wvar   = None
- wopts  = None

For more information on these options, see quad and quad_explain.

full_output : bool, optional
Partial implementation of full_output from scipy.integrate.quad.
The number of integrand function evaluations neval can be obtained
by setting full_output=True when calling nquad.

Returns
-------
result : float
The result of the integration.
abserr : float
The maximum of the estimates of the absolute error in the various
integration results.
out_dict : dict, optional
A dict containing additional information on the integration.

--------

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 + 0.5 + 0.25*args]}
>>> integrate.nquad(func, [[0,1], [-1,1], [.13,.8], [-.15,1]],
...                 opts=[opts0,{},{},{}], full_output=True)
(1.5267454070738633, 2.9437360001402324e-14, {'neval': 388962})

>>> scale = .1
>>> def func2(x0, x1, x2, x3, t0, t1):
...     return x0*x1*x3**2 + np.sin(x2) + 1 + (1 if x0+t1*x1-t0>0 else 0)
>>> def lim0(x1, x2, x3, t0, t1):
...     return [scale * (x1**2 + x2 + np.cos(x3)*t0*t1 + 1) - 1,
...             scale * (x1**2 + x2 + np.cos(x3)*t0*t1 + 1) + 1]
>>> def lim1(x2, x3, t0, t1):
...     return [scale * (t0*x2 + t1*x3) - 1,
...             scale * (t0*x2 + t1*x3) + 1]
>>> def lim2(x3, t0, t1):
...     return [scale * (x3 + t0**2*t1**3) - 1,
...             scale * (x3 + t0**2*t1**3) + 1]
>>> def lim3(t0, t1):
...     return [scale * (t0+t1) - 1, scale * (t0+t1) + 1]
>>> def opts0(x1, x2, x3, t0, t1):
...     return {'points' : [t0 - t1*x1]}
>>> def opts1(x2, x3, t0, t1):
...     return {}
>>> def opts2(x3, t0, t1):
...     return {}
>>> def opts3(t0, t1):
...     return {}
>>> integrate.nquad(func2, [lim0, lim1, lim2, lim3], args=(0,0),
...                 opts=[opts0, opts1, opts2, opts3])
(25.066666666666666, 2.7829590483937256e-13)

odeint(func, y0, t, args=(), Dfun=None, col_deriv=0, full_output=0, ml=None, mu=None, rtol=None, atol=None, tcrit=None, h0=0.0, hmax=0.0, hmin=0.0, ixpr=0, mxstep=0, mxhnil=0, mxordn=12, mxords=5, printmessg=0)
Integrate a system of ordinary differential equations.

Solve a system of ordinary differential equations using lsoda from the
FORTRAN library odepack.

Solves the initial value problem for stiff or non-stiff systems
of first order ode-s::

dy/dt = func(y, t0, ...)

where y can be a vector.

*Note*: The first two arguments of func(y, t0, ...) are in the
opposite order of the arguments in the system definition function used
by the scipy.integrate.ode class.

Parameters
----------
func : callable(y, t0, ...)
Computes the derivative of y at t0.
y0 : array
Initial condition on y (can be a vector).
t : array
A sequence of time points for which to solve for y.  The initial
value point should be the first element of this sequence.
args : tuple, optional
Extra arguments to pass to function.
Dfun : callable(y, t0, ...)
Gradient (Jacobian) of func.
col_deriv : bool, optional
True if Dfun defines derivatives down columns (faster),
otherwise Dfun should define derivatives across rows.
full_output : bool, optional
True if to return a dictionary of optional outputs as the second output
printmessg : bool, optional
Whether to print the convergence message

Returns
-------
y : array, shape (len(t), len(y0))
Array containing the value of y for each desired time in t,
with the initial value y0 in the first row.
infodict : dict, only returned if full_output == True

=======  ============================================================
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 ith equation with respect to the jth
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.

--------
ode : a more object-oriented integrator based on VODE.
quad : for finding the area under a curve.

Examples
--------
The second order differential equation for the angle theta of a
pendulum acted on by gravity with friction can be written::

theta''(t) + b*theta'(t) + c*sin(theta(t)) = 0

where b and c are positive constants, and a prime (') denotes a
derivative.  To solve this equation with odeint, we must first convert
it to a system of first order equations.  By defining the angular
velocity omega(t) = theta'(t), we obtain the system::

theta'(t) = omega(t)
omega'(t) = -b*omega(t) - c*sin(theta(t))

Let y be the vector [theta, omega].  We implement this system
in python as:

>>> def pend(y, t, b, c):
...     theta, omega = y
...     dydt = [omega, -b*omega - c*np.sin(theta)]
...     return dydt
...

We assume the constants are b = 0.25 and c = 5.0:

>>> b = 0.25
>>> c = 5.0

For initial conditions, we assume the pendulum is nearly vertical
with theta(0) = pi - 0.1, and it initially at rest, so
omega(0) = 0.  Then the vector of initial conditions is

>>> y0 = [np.pi - 0.1, 0.0]

We generate a solution 101 evenly spaced samples in the interval
0 <= t <= 10.  So our array of times is:

>>> t = np.linspace(0, 10, 101)

Call odeint to generate the solution.  To pass the parameters
b and c to pend, we give them to odeint using the args
argument.

>>> from scipy.integrate import odeint
>>> sol = odeint(pend, y0, t, args=(b, c))

The solution is an array with shape (101, 2).  The first column
is theta(t), and the second is omega(t).  The following code
plots both components.

>>> import matplotlib.pyplot as plt
>>> plt.plot(t, sol[:, 0], 'b', label='theta(t)')
>>> plt.plot(t, sol[:, 1], 'g', label='omega(t)')
>>> plt.legend(loc='best')
>>> plt.xlabel('t')
>>> plt.grid()
>>> plt.show()

quad(func, a, b, args=(), full_output=0, epsabs=1.49e-08, epsrel=1.49e-08, limit=50, points=None, weight=None, wvar=None, wopts=None, maxp1=50, limlst=50)
Compute a definite integral.

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

Parameters
----------
func : {function, 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 == 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
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.

--------
nquad : n-dimensional integrals (uses quad recursively)
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], ..., 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, pts) where pts and pts
are adjacent elements of infodict['pts'], then (aa,bb) has level l
if |bb-aa| = |pts-pts| * 2**(-l).
'ndin'
A rank-1 integer array of length P+2.  After the first integration
over the intervals (pts, pts), 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
(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)
(1.0, 5.842605999138044e-11)

>>> f = lambda x,a : a*x
>>> y, err = integrate.quad(f, 0, 1, args=(1,))
>>> y
0.5
>>> y, err = integrate.quad(f, 0, 1, args=(3,))
>>> y
1.5

Calculate :math:\int^1_0 x^2 + y^2 dx with ctypes, holding
y parameter as 1::

testlib.c =>
double func(int n, double args[n]){
return args*args + args*args;}
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)
#(1.3333333333333333, 1.4802973661668752e-14)
print((1.0**3/3.0 + 1.0) - (0.0**3/3.0 + 0.0)) #Analytic result
# 1.3333333333333333

Parameters
----------
output : instance with "write" method, optional
Information about quad is passed to output.write().
Default is sys.stdout.

Returns
-------
None

quadrature(func, a, b, args=(), tol=1.49e-08, rtol=1.49e-08, maxiter=50, vec_func=True, miniter=1)
Compute a definite integral using fixed-tolerance Gaussian quadrature.

Integrate func from a to b using Gaussian quadrature
with absolute tolerance tol.

Parameters
----------
func : function
A Python function or method to integrate.
a : float
Lower limit of integration.
b : float
Upper limit of integration.
args : tuple, optional
Extra arguments to pass to function.
tol, rtol : float, optional
Iteration stops when error between last two iterates is less than
tol OR the relative change is less than rtol.
maxiter : int, optional
vec_func : bool, optional
True or False if func handles arrays as arguments (is
a "vector" function). Default is True.
miniter : int, optional

Returns
-------
val : float
Gaussian quadrature approximation (within tolerance) to integral.
err : float
Difference between last two estimates of the integral.

--------
romb: integrator for sampled data
simps: integrator for sampled data
cumtrapz: cumulative integration for sampled data
ode: ODE integrator
odeint: ODE integrator

romb(y, dx=1.0, axis=-1, show=False)
Romberg integration using samples of a function.

Parameters
----------
y : array_like
A vector of 2**k + 1 equally-spaced samples of a function.
dx : float, optional
The sample spacing. Default is 1.
axis : int, optional
The axis along which to integrate. Default is -1 (last axis).
show : bool, optional
When y is a single 1-D array, then if this argument is True
print the table showing Richardson extrapolation from the
samples. Default is False.

Returns
-------
romb : ndarray
The integrated result for axis.

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

--------
romb : Integrators for sampled data.
simps : Integrators for sampled data.
cumtrapz : Cumulative integration for sampled data.
ode : ODE integrator.
odeint : ODE integrator.

References
----------
..  '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 : 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.

--------
romb: integrators for sampled data
cumtrapz: cumulative integration for sampled data
ode: ODE integrators
odeint: ODE integrators

Notes
-----
For an odd number of samples that are equally spaced the result is
exact if the function is a polynomial of order 3 or less.  If
the samples are not equally spaced, then the result is exact only
if the function is a polynomial of order 2 or less.

solve_bvp(fun, bc, x, y, p=None, S=None, fun_jac=None, bc_jac=None, tol=0.001, max_nodes=1000, verbose=0)
Solve a boundary-value problem for a system of ODEs.

This function numerically solves a first order system of ODEs subject to
two-point boundary conditions::

dy / dx = f(x, y, p) + S * y / (x - a), a <= x <= b
bc(y(a), y(b), p) = 0

Here x is a 1-dimensional independent variable, y(x) is a n-dimensional
vector-valued function and p is a k-dimensional vector of unknown
parameters which is to be found along with y(x). For the problem to be
determined there must be n + k boundary conditions, i.e. bc must be
(n + k)-dimensional function.

The last singular term in the right-hand side of the system is optional.
It is defined by an n-by-n matrix S, such that the solution must satisfy
S y(a) = 0. This condition will be forced during iterations, so it must not
contradict boundary conditions. See _ for the explanation how this term
is handled when solving BVPs numerically.

Problems in a complex domain can be solved as well. In this case y and p
are considered to be complex, and f and bc are assumed to be complex-valued
functions, but x stays real. Note that f and bc must be complex
differentiable (satisfy Cauchy-Riemann equations _), otherwise you
should rewrite your problem for real and imaginary parts separately. To
solve a problem in a complex domain, pass an initial guess for y with a
complex data type (see below).

Parameters
----------
fun : callable
Right-hand side of the system. The calling signature is fun(x, y),
or fun(x, y, p) if parameters are present. All arguments are
ndarray: x with shape (m,), y with shape (n, m), meaning that
y[:, i] corresponds to x[i], and p with shape (k,). The
return value must be an array with shape (n, m) and with the same
layout as y.
bc : callable
Function evaluating residuals of the boundary conditions. The calling
signature is bc(ya, yb), or bc(ya, yb, p) if parameters are
present. All arguments are ndarray: ya and yb with shape (n,),
and p with shape (k,). The return value must be an array with
shape (n + k,).
x : array_like, shape (m,)
Initial mesh. Must be a strictly increasing sequence of real numbers
with x=a and x[-1]=b.
y : array_like, shape (n, m)
Initial guess for the function values at the mesh nodes, i-th column
corresponds to x[i]. For problems in a complex domain pass y
with a complex data type (even if the initial guess is purely real).
p : array_like with shape (k,) or None, optional
Initial guess for the unknown parameters. If None (default), it is
assumed that the problem doesn't depend on any parameters.
S : array_like with shape (n, n) or None
Matrix defining the singular term. If None (default), the problem is
solved without the singular term.
fun_jac : callable or None, optional
Function computing derivatives of f with respect to y and p. The
calling signature is fun_jac(x, y), or fun_jac(x, y, p) if
parameters are present. The return must contain 1 or 2 elements in the
following order:

* df_dy : array_like with shape (n, n, m) where an element
(i, j, q) equals to d f_i(x_q, y_q, p) / d (y_q)_j.
* df_dp : array_like with shape (n, k, m) where an element
(i, j, q) equals to d f_i(x_q, y_q, p) / d p_j.

Here q numbers nodes at which x and y are defined, whereas i and j
number vector components. If the problem is solved without unknown
parameters df_dp should not be returned.

If fun_jac is None (default), the derivatives will be estimated
by the forward finite differences.
bc_jac : callable or None, optional
Function computing derivatives of bc with respect to ya, yb and p.
The calling signature is bc_jac(ya, yb), or bc_jac(ya, yb, p)
if parameters are present. The return must contain 2 or 3 elements in
the following order:

* dbc_dya : array_like with shape (n, n) where an element (i, j)
equals to d bc_i(ya, yb, p) / d ya_j.
* dbc_dyb : array_like with shape (n, n) where an element (i, j)
equals to d bc_i(ya, yb, p) / d yb_j.
* dbc_dp : array_like with shape (n, k) where an element (i, j)
equals to d bc_i(ya, yb, p) / d p_j.

If the problem is solved without unknown parameters dbc_dp should not
be returned.

If bc_jac is None (default), the derivatives will be estimated by
the forward finite differences.
tol : float, optional
Desired tolerance of the solution. If we define r = y' - f(x, y)
where y is the found solution, then the solver tries to achieve on each
mesh interval norm(r / (1 + abs(f)) < tol, where norm is
estimated in a root mean squared sense (using a numerical quadrature
formula). Default is 1e-3.
max_nodes : int, optional
Maximum allowed number of the mesh nodes. If exceeded, the algorithm
terminates. Default is 1000.
verbose : {0, 1, 2}, optional
Level of algorithm's verbosity:

* 0 (default) : work silently.
* 1 : display a termination report.
* 2 : display progress during iterations.

Returns
-------
Bunch object with the following fields defined:
sol : PPoly
Found solution for y as scipy.interpolate.PPoly instance, a C1
continuous cubic spline.
p : ndarray or None, shape (k,)
Found parameters. None, if the parameters were not present in the
problem.
x : ndarray, shape (m,)
Nodes of the final mesh.
y : ndarray, shape (n, m)
Solution values at the mesh nodes.
yp : ndarray, shape (n, m)
Solution derivatives at the mesh nodes.
rms_residuals : ndarray, shape (m - 1,)
RMS values of the relative residuals over each mesh interval (see the
description of tol parameter).
niter : int
Number of completed iterations.
status : int
Reason for algorithm termination:

* 0: The algorithm converged to the desired accuracy.
* 1: The maximum number of mesh nodes is exceeded.
* 2: A singular Jacobian encountered when solving the collocation
system.

message : string
Verbal description of the termination reason.
success : bool
True if the algorithm converged to the desired accuracy (status=0).

Notes
-----
This function implements a 4-th order collocation algorithm with the
control of residuals similar to _. A collocation system is solved
by a damped Newton method with an affine-invariant criterion function as
described in _.

Note that in _  integral residuals are defined without normalization
by interval lengths. So their definition is different by a multiplier of
h**0.5 (h is an interval length) from the definition used here.

References
----------
..  J. Kierzenka, L. F. Shampine, "A BVP Solver Based on Residual
Control and the Maltab PSE", ACM Trans. Math. Softw., Vol. 27,
Number 3, pp. 299-316, 2001.
..  L.F. Shampine, P. H. Muir and H. Xu, "A User-Friendly Fortran BVP
Solver".
..  U. Ascher, R. Mattheij and R. Russell "Numerical Solution of
Boundary Value Problems for Ordinary Differential Equations".
..  Cauchy-Riemann equations
<https://en.wikipedia.org/wiki/Cauchy-Riemann_equations>_ on
Wikipedia.

Examples
--------
In the first example we solve Bratu's problem::

y'' + k * exp(y) = 0
y(0) = y(1) = 0

for k = 1.

We rewrite the equation as a first order system and implement its
right-hand side evaluation::

y1' = y2
y2' = -exp(y1)

>>> def fun(x, y):
...     return np.vstack((y, -np.exp(y)))

Implement evaluation of the boundary condition residuals:

>>> def bc(ya, yb):
...     return np.array([ya, yb])

Define the initial mesh with 5 nodes:

>>> x = np.linspace(0, 1, 5)

This problem is known to have two solutions. To obtain both of them we
use two different initial guesses for y. We denote them by subscripts
a and b.

>>> y_a = np.zeros((2, x.size))
>>> y_b = np.zeros((2, x.size))
>>> y_b = 3

Now we are ready to run the solver.

>>> from scipy.integrate import solve_bvp
>>> res_a = solve_bvp(fun, bc, x, y_a)
>>> res_b = solve_bvp(fun, bc, x, y_b)

Let's plot the two found solutions. We take an advantage of having the
solution in a spline form to produce a smooth plot.

>>> x_plot = np.linspace(0, 1, 100)
>>> y_plot_a = res_a.sol(x_plot)
>>> y_plot_b = res_b.sol(x_plot)
>>> import matplotlib.pyplot as plt
>>> plt.plot(x_plot, y_plot_a, label='y_a')
>>> plt.plot(x_plot, y_plot_b, label='y_b')
>>> plt.legend()
>>> plt.xlabel("x")
>>> plt.ylabel("y")
>>> plt.show()

We see that the two solutions have similar shape, but differ in scale
significantly.

In the second example we solve a simple Sturm-Liouville problem::

y'' + k**2 * y = 0
y(0) = y(1) = 0

It is known that a non-trivial solution y = A * sin(k * x) is possible for
k = pi * n, where n is an integer. To establish the normalization constant
A = 1 we add a boundary condition::

y'(0) = k

Again we rewrite our equation as a first order system and implement its
right-hand side evaluation::

y1' = y2
y2' = -k**2 * y1

>>> def fun(x, y, p):
...     k = p
...     return np.vstack((y, -k**2 * y))

Note that parameters p are passed as a vector (with one element in our
case).

Implement the boundary conditions:

>>> def bc(ya, yb, p):
...     k = p
...     return np.array([ya, yb, ya - k])

Setup the initial mesh and guess for y. We aim to find the solution for
k = 2 * pi, to achieve that we set values of y to approximately follow
sin(2 * pi * x):

>>> x = np.linspace(0, 1, 5)
>>> y = np.zeros((2, x.size))
>>> y[0, 1] = 1
>>> y[0, 3] = -1

Run the solver with 6 as an initial guess for k.

>>> sol = solve_bvp(fun, bc, x, y, p=)

We see that the found k is approximately correct:

>>> sol.p
6.28329460046

And finally plot the solution to see the anticipated sinusoid:

>>> x_plot = np.linspace(0, 1, 100)
>>> y_plot = sol.sol(x_plot)
>>> plt.plot(x_plot, y_plot)
>>> plt.xlabel("x")
>>> plt.ylabel("y")
>>> plt.show()

solve_ivp(fun, t_span, y0, method='RK45', t_eval=None, dense_output=False, events=None, vectorized=False, **options)
Solve an initial value problem for a system of ODEs.

This function numerically integrates a system of ordinary differential
equations given an initial value::

dy / dt = f(t, y)
y(t0) = y0

Here t is a 1-dimensional independent variable (time), y(t) is an
n-dimensional vector-valued function (state) and an n-dimensional
vector-valued function f(t, y) determines the differential equations.
The goal is to find y(t) approximately satisfying the differential
equations, given an initial value y(t0)=y0.

Some of the solvers support integration in a complex domain, but note that
for stiff ODE solvers the right hand side must be complex differentiable
(satisfy Cauchy-Riemann equations _). To solve a problem in a complex
domain, pass y0 with a complex data type. Another option always available
is to rewrite your problem for real and imaginary parts separately.

Parameters
----------
fun : callable
Right-hand side of the system. The calling signature is fun(t, y).
Here t is a scalar and there are two options for ndarray y.
It can either have shape (n,), then fun must return array_like with
shape (n,). Or alternatively it can have shape (n, k), then fun
must return array_like with shape (n, k), i.e. each column
corresponds to a single column in y. The choice between the two
options is determined by vectorized argument (see below). The
vectorized implementation allows faster approximation of the Jacobian
by finite differences (required for stiff solvers).
t_span : 2-tuple of floats
Interval of integration (t0, tf). The solver starts with t=t0 and
integrates until it reaches t=tf.
y0 : array_like, shape (n,)
Initial state. For problems in a complex domain pass y0 with a
complex data type (even if the initial guess is purely real).
method : string or OdeSolver, optional
Integration method to use:

* 'RK45' (default): Explicit Runge-Kutta method of order 5(4) _.
The error is controlled assuming 4th order accuracy, but steps
are taken using a 5th oder accurate formula (local extrapolation
is done). A quartic interpolation polynomial is used for the
dense output _. Can be applied in a complex domain.
* 'RK23': Explicit Runge-Kutta method of order 3(2) _. The error
is controlled assuming 2nd order accuracy, but steps are taken
using a 3rd oder accurate formula (local extrapolation is done).
A cubic Hermit polynomial is used for the dense output.
Can be applied in a complex domain.
order 5 _. The error is controlled for a 3rd order accurate
embedded formula. A cubic polynomial which satisfies the
collocation conditions is used for the dense output.
* 'BDF': Implicit multi-step variable order (1 to 5) method based
on a Backward Differentiation Formulas for the derivative
approximation _. An implementation approach follows the one
described in _. A quasi-constant step scheme is used
and accuracy enhancement using NDF modification is also
implemented. Can be applied in a complex domain.
* 'LSODA': Adams/BDF method with automatic stiffness detection and
switching _, _. This is a wrapper of the Fortran solver
from ODEPACK.

You should use 'RK45' or 'RK23' methods for non-stiff problems and
'Radau' or 'BDF' for stiff problems _. If not sure, first try to run
'RK45' and if it does unusual many iterations or diverges then your
problem is likely to be stiff and you should use 'Radau' or 'BDF'.
'LSODA' can also be a good universal choice, but it might be somewhat
less  convenient to work with as it wraps an old Fortran code.

You can also pass an arbitrary class derived from OdeSolver which
implements the solver.
dense_output : bool, optional
Whether to compute a continuous solution. Default is False.
t_eval : array_like or None, optional
Times at which to store the computed solution, must be sorted and lie
within t_span. If None (default), use points selected by a solver.
events : callable, list of callables or None, optional
Events to track. Events are defined by functions which take
a zero value at a point of an event. Each function must have a
signature event(t, y) and return float, the solver will find an
accurate value of t at which event(t, y(t)) = 0 using a root
finding algorithm. Additionally each event function might have
attributes:

* terminal: bool, whether to terminate integration if this
event occurs. Implicitly False if not assigned.
* direction: float, direction of crossing a zero. If direction
is positive then event must go from negative to positive, and
vice-versa if direction is negative. If 0, then either way will
count. Implicitly 0 if not assigned.

You can assign attributes like event.terminal = True to any
function in Python. If None (default), events won't be tracked.
vectorized : bool, optional
Whether fun is implemented in a vectorized fashion. Default is False.
options
Options passed to a chosen solver constructor. All options available
for already implemented solvers are listed below.
max_step : float, optional
Maximum allowed step size. Default is np.inf, i.e. step is not
bounded and determined solely by the solver.
rtol, atol : float and array_like, optional
Relative and absolute tolerances. The solver keeps the local error
estimates less than atol + rtol * abs(y). Here rtol controls a
relative accuracy (number of correct digits). But if a component of y
is approximately below atol then the error only needs to fall within
the same atol threshold, and the number of correct digits is not
guaranteed. If components of y have different scales, it might be
beneficial to set different atol values for different components by
passing array_like with shape (n,) for atol. Default values are
1e-3 for rtol and 1e-6 for atol.
jac : {None, array_like, sparse_matrix, callable}, optional
Jacobian matrix of the right-hand side of the system with respect to
y, required by 'Radau', 'BDF' and 'LSODA' methods. The Jacobian matrix
has shape (n, n) and its element (i, j) is equal to d f_i / d y_j.
There are 3 ways to define the Jacobian:

* If array_like or sparse_matrix, then the Jacobian is assumed to
be constant. Not supported by 'LSODA'.
* If callable, then the Jacobian is assumed to depend on both
t and y, and will be called as jac(t, y) as necessary.
For 'Radau' and 'BDF' methods the return value might be a sparse
matrix.
* If None (default), then the Jacobian will be approximated by
finite differences.

It is generally recommended to provide the Jacobian rather than
relying on a finite difference approximation.
jac_sparsity : {None, array_like, sparse matrix}, optional
Defines a sparsity structure of the Jacobian matrix for a finite
difference approximation, its shape must be (n, n). If the Jacobian has
only few non-zero elements in *each* row, providing the sparsity
structure will greatly speed up the computations _. A zero
entry means that a corresponding element in the Jacobian is identically
zero. If None (default), the Jacobian is assumed to be dense.
Not supported by 'LSODA', see lband and uband instead.
lband, uband : int or None
Parameters defining the Jacobian matrix bandwidth for 'LSODA' method.
The Jacobian bandwidth means that
jac[i, j] != 0 only for i - lband <= j <= i + uband. Setting these
requires your jac routine to return the Jacobian in the packed format:
the returned array must have n columns and uband + lband + 1
rows in which Jacobian diagonals are written. Specifically
jac_packed[uband + i - j , j] = jac[i, j]. The same format is used
in scipy.linalg.solve_banded (check for an illustration).
These parameters can be also used with jac=None to reduce the
number of Jacobian elements estimated by finite differences.
min_step, first_step : float, optional
The minimum allowed step size and the initial step size respectively
for 'LSODA' method. By default min_step is zero and first_step is
selected automatically.

Returns
-------
Bunch object with the following fields defined:
t : ndarray, shape (n_points,)
Time points.
y : ndarray, shape (n, n_points)
Solution values at t.
sol : OdeSolution or None
Found solution as OdeSolution instance, None if dense_output was
set to False.
t_events : list of ndarray or None
Contains arrays with times at each a corresponding event was detected,
the length of the list equals to the number of events. None if events
was None.
nfev : int
Number of the system rhs evaluations.
njev : int
Number of the Jacobian evaluations.
nlu : int
Number of LU decompositions.
status : int
Reason for algorithm termination:

* -1: Integration step failed.
* 0: The solver successfully reached the interval end.
* 1: A termination event occurred.

message : string
Verbal description of the termination reason.
success : bool
True if the solver reached the interval end or a termination event
occurred (status >= 0).

References
----------
..  J. R. Dormand, P. J. Prince, "A family of embedded Runge-Kutta
formulae", Journal of Computational and Applied Mathematics, Vol. 6,
No. 1, pp. 19-26, 1980.
..  L. W. Shampine, "Some Practical Runge-Kutta Formulas", Mathematics
of Computation,, Vol. 46, No. 173, pp. 135-150, 1986.
..  P. Bogacki, L.F. Shampine, "A 3(2) Pair of Runge-Kutta Formulas",
Appl. Math. Lett. Vol. 2, No. 4. pp. 321-325, 1989.
..  E. Hairer, G. Wanner, "Solving Ordinary Differential Equations II:
Stiff and Differential-Algebraic Problems", Sec. IV.8.
..  Backward Differentiation Formula
<https://en.wikipedia.org/wiki/Backward_differentiation_formula>_
on Wikipedia.
..  L. F. Shampine, M. W. Reichelt, "THE MATLAB ODE SUITE", SIAM J. SCI.
COMPUTE., Vol. 18, No. 1, pp. 1-22, January 1997.
..  A. C. Hindmarsh, "ODEPACK, A Systematized Collection of ODE
Solvers," IMACS Transactions on Scientific Computation, Vol 1.,
pp. 55-64, 1983.
..  L. Petzold, "Automatic selection of methods for solving stiff and
nonstiff systems of ordinary differential equations", SIAM Journal
on Scientific and Statistical Computing, Vol. 4, No. 1, pp. 136-148,
1983.
..  Stiff equation <https://en.wikipedia.org/wiki/Stiff_equation>_ on
Wikipedia.
..  A. Curtis, M. J. D. Powell, and J. Reid, "On the estimation of
sparse Jacobian matrices", Journal of the Institute of Mathematics
and its Applications, 13, pp. 117-120, 1974.
..  Cauchy-Riemann equations
<https://en.wikipedia.org/wiki/Cauchy-Riemann_equations>_ on
Wikipedia.

Examples
--------
Basic exponential decay showing automatically chosen time points.

>>> from scipy.integrate import solve_ivp
>>> def exponential_decay(t, y): return -0.5 * y
>>> sol = solve_ivp(exponential_decay, [0, 10], [2, 4, 8])
>>> print(sol.t)
[  0.           0.11487653   1.26364188   3.06061781   4.85759374
6.65456967   8.4515456   10.        ]
>>> print(sol.y)
[[ 2.          1.88836035  1.06327177  0.43319312  0.17648948  0.0719045
0.02929499  0.01350938]
[ 4.          3.7767207   2.12654355  0.86638624  0.35297895  0.143809
0.05858998  0.02701876]
[ 8.          7.5534414   4.25308709  1.73277247  0.7059579   0.287618
0.11717996  0.05403753]]

Specifying points where the solution is desired.

>>> sol = solve_ivp(exponential_decay, [0, 10], [2, 4, 8],
...                 t_eval=[0, 1, 2, 4, 10])
>>> print(sol.t)
[ 0  1  2  4 10]
>>> print(sol.y)
[[ 2.          1.21305369  0.73534021  0.27066736  0.01350938]
[ 4.          2.42610739  1.47068043  0.54133472  0.02701876]
[ 8.          4.85221478  2.94136085  1.08266944  0.05403753]]

Cannon fired upward with terminal event upon impact. The terminal and
direction fields of an event are applied by monkey patching a function.
Here y is position and y is velocity. The projectile starts at
position 0 with velocity +10. Note that the integration never reaches t=100
because the event is terminal.

>>> def upward_cannon(t, y): return [y, -0.5]
>>> def hit_ground(t, y): return y
>>> hit_ground.terminal = True
>>> hit_ground.direction = -1
>>> sol = solve_ivp(upward_cannon, [0, 100], [0, 10], events=hit_ground)
>>> print(sol.t_events)
[array([ 20.])]
>>> print(sol.t)
[  0.00000000e+00   9.99900010e-05   1.09989001e-03   1.10988901e-02
1.11088891e-01   1.11098890e+00   1.11099890e+01   2.00000000e+01]

tplquad(func, a, b, gfun, hfun, qfun, rfun, args=(), epsabs=1.49e-08, epsrel=1.49e-08)
Compute a triple (definite) integral.

Return the triple integral of func(z, y, x) from x = a..b,
y = gfun(x)..hfun(x), and z = qfun(x,y)..rfun(x,y).

Parameters
----------
func : function
A Python function or method of at least three variables in the
order (z, y, x).
a, b : float
The limits of integration in x: a < b
gfun : function
The lower boundary curve in y which is a function taking a single
floating point argument (x) and returning a floating point result:
a lambda function can be useful here.
hfun : function
The upper boundary curve in y (same requirements as gfun).
qfun : function
The lower boundary surface in z.  It must be a function that takes
two floats in the order (x, y) and returns a float.
rfun : function
The upper boundary surface in z. (Same requirements as qfun.)
args : tuple, optional
Extra arguments to pass to func.
epsabs : float, optional
Absolute tolerance passed directly to the innermost 1-D quadrature
integration. Default is 1.49e-8.
epsrel : float, optional
Relative tolerance of the innermost 1-D integrals. Default is 1.49e-8.

Returns
-------
y : float
The resultant integral.
abserr : float
An estimate of the error.

--------
romb: Integrators for sampled data
simps: Integrators for sampled data
ode: ODE integrators
odeint: ODE integrators
scipy.special: For coefficients and roots of orthogonal polynomials

trapz(y, x=None, dx=1.0, axis=-1)
Integrate along the given axis using the composite trapezoidal rule.

Integrate y (x) along given axis.

Parameters
----------
y : array_like
Input array to integrate.
x : array_like, optional
The sample points corresponding to the y values. If x is None,
the sample points are assumed to be evenly spaced dx apart. The
default is None.
dx : scalar, optional
The spacing between sample points when x is None. The default is 1.
axis : int, optional
The axis along which to integrate.

Returns
-------
trapz : float
Definite integral as approximated by trapezoidal rule.

--------
sum, cumsum

Notes
-----
Image _ 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
----------

..  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__ = ['BDF', 'DenseOutput', 'IntegrationWarning', 'LSODA', 'OdeSo...
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/minnolina/anaconda3/lib/python3.6/site-packages/scipy/integrate/__init__.py



### odeint¶

L'ECUE M62 n'est qu'une introduction aux méthodes numériques d'approximation d'EDO en utilisant Python comme langage commun de programmation. Cependant une grande partie des méthodes numériques dont vous avez besoin en tant que scientifique (y compris l'intégration numérique d'équations différentielles) se trouve déjà dans le module SciPy.
Bien sûr vous devez toujours faire un peu de programmation avec Python et une compréhension des fondements de la méthode numérique que vous utilisez est toujours indispensables mais avouons que c'est commode quand quelque chose est déjà programmé.

Voyons sur deux exemples comment utiliser la fonction odeint du module SciPy pour approcher la solution d'abord d'une EDO ensuite d'un système d'EDO (ce qui inclut les équations différentielles d'ordre 2 ou plus).

In :
from matplotlib.pylab import *
from scipy.integrate import odeint


Le principe d’utilisation de odeint (pour intégrer numériquement des équations différentielles) est le suivant: pour avoir une estimation numérique de la solution du problème de $$\begin{cases} \mathbf{y}'(t)=\boldsymbol{\varphi}(t,\mathbf{y}(t)),\\ \mathbf{y}(t_0)=\mathbf{y}_0 \end{cases}$$ avec $\mathbf{y}(t)=(y_1(t),y_2(t),\dots,y_n(t))$ le vecteur des fonctions recherchées, dépendant de la variable $t$ et $\boldsymbol{\varphi}=(\varphi_1,\varphi_2,\dots,\varphi_n)$ une fonction de forme quelconque, on donne comme argument la fonction $\boldsymbol{\varphi}$ (qui doit avoir deux paramètres, même dans le cas autonome, avec $t$ comme deuxième paramètre), la condition initiale $\mathbf{y}_0$ et le domaine de temps qui nous intéresse (qui commence à $t_0$). Elle retourne un tableau (même si $t$ était une liste).

Notons que la résolution de $y''(t)=F(t,y(t),y'(t))$ passera par celle du système différentiel $$\begin{cases} y_1'(t)=y_2(t),\\ y_2'(t)=F(t,y_1(t),y_2(t)) \end{cases}$$ avec $y(t)=y_1(t)$, $y'(t)=y_1'(t)=y_2(t)$ et $\varphi_1(t,y_1(t),y_2(t))=y_2(t)$, $\varphi_2(t,y_1(t),y_2(t))=F(t,y_1(t),y_2(t))$.

#### Exemple de résolution approchée d’une équation différentielle¶

En guise d’exemple, considérons une équation logistique simple de la forme $$y'(t)=r y(t) \left( 1-\frac{y(t)}{K} \right).$$ Pour les simulations on prendra $r=1.5$, $K=6$ et $y_0=1$. On crée alors la fonction $\varphi$

In :
# def phi(y,t):
#     return 1.5*y*(1.-y/6.)
# ou
phi = lambda y,t : 1.5*y*(1.-y/6.)


La fonction odeint peut être appelée avec au minimum trois arguments: la fonction $\varphi$, la condition initiale $y(t_0)=y_0$ et le temps $t$ comme variable principale:

In :
y0  = 1.0
t   = linspace(0,5,201)
sol = odeint(phi,y0,t)


La solution peut alors être tracée simplement

In :
plot(t,sol)
title(r'$y\prime(t)=\frac{3}{2}y(t)\left(1-\frac{y(t)}{6}\right)$ et $y(0)=1$')
show() NB: le nombre de points en lesquels les résultats sont évalués n’est pas (du moins directement) relié à la précision des calculs internes (ne pas imaginer que cela fixe le pas de la méthode, en particulier).

Notons que la fonction quiver du module matplotlib permet de tracer un champ de vecteurs. Utilisons-la pour obtenir celui associé à notre équation différentielle. La courbe bleu est la solution déterminée par odeint.

In :
T,Y = np.meshgrid(linspace(0,5,31),linspace(0,8,31))
U = 1.
V = phi(Y,T)
r=sqrt(U**2+V**2)
quiver(T, Y, U/r, V/r)
plot(t,sol)
grid(True)
title(r'Champ des pentes'); #### Exemple de résolution approchée d’un système d'équations différentielles¶

Proposées indépendamment par Alfred James Lotka en 1925 et Vito Volterra en 1926, les équations dites de Lotka-Volterra ou modèle proie-prédateur, peuvent servir à reproduire les évolutions temporelles de deux populations animales, l’une étant le prédateur de l'autre.

Si ces deux populations sont représentées par des variables continues $y_1(t)$ et $y_2(t)$, le système différentiel en question est alors $$\begin{pmatrix} y_1\ y_2 # \end{pmatrix}'(t)¶ \begin{pmatrix} \varphi_1(t, y_1(t),y_2(t))\ \varphi_2(t, y_1(t),y_2(t)) # \end{pmatrix}¶ \begin{pmatrix} y_1(t)(a-by_2(t))\\ -y_2(t)(c-dy_1(t)) \end{pmatrix}.$$ Pour les simulations on prendra $a=2$, $b=1$, $c=1$ et $d=0.3$. On crée alors la fonction vectorielle $\boldsymbol{\varphi}=(\varphi_1,\varphi_2)$:

In :
# def phi(y,t):
#     return [ y*(2.-y) , -y*(1.-0.3*y) ]
# ou
phi = lambda y,t : [ y*(2.-y) , -y*(1.-0.3*y) ]


En faisant varier le temps sur l'intervalle $[0;10]$ et en prenant comme condition initiale le vecteur $\mathbf{y}_0=(2,1)$ on écrit

In :
y0  = [2.0,1.0]
t   = linspace(0,10,201)
sol = odeint(phi,y0,t)


Le tracé des évolutions de $y_1$ et $y_2$ en fonction du temps $t$ peut être obtenu par

In :
plot(t,sol)
legend(["y1","y2"]); Le tracé des évolutions de $y_2$ en fonction de $y_1$ peut être obtenu par

In :
plot(sol[:,0],sol[:,1]); #### Exemple de résolution approchée d’une équation différentielle d'ordre 2¶

Considérons l'EDO $y''(t)=-\sin(y(t))$ qui décrit le mouvement d'un pendule non amorti, équivalente au système différentiel $$\begin{cases} y_1'(t)=y_2(t),\\ y_2'(t)=-\sin(y_1(t)) \end{cases}$$ avec $y(0)=0$ et $y'(0)=1$.

In :
phi = lambda y,t : [ y, -math.sin(y)  ]
y0  = [0,1.0]
t   = linspace(0,10,1001)
sol = odeint(phi,y0,t)
figure(1, figsize=(10, 5))
subplot(1,2,1)
plot(t,sol)
subplot(1,2,2)
plot(sol[:,0],sol[:,1]); ## Calcul formel avec le module sympy¶

SymPy est une bibliothèque Python pour les mathématiques symboliques. Elle prévoit devenir un système complet de calcul formel (“CAS” en anglais : “Computer Algebra System”) tout en gardant le code aussi simple que possible afin qu’il soit compréhensible et facilement extensible.

Quelques commandes : https://www.sympygamma.com/

Pour tester en ligne: https://live.sympy.org/

In :
%reset -f
from sympy import *


Pour un meilleur aspect des résultats affichées on utilisera :

In :
init_printing()


Le symbole oo est utilisé pour une classe définissant l’infini mathématique :

In :
oo+1

Out:
$$\infty$$

Contrairement à d’autres systèmes de calcul formel, SymPy vous oblige à déclarer les variables symboliques explicitement :

In :
#  premiere methode pour declarer des variables
x = Symbol('x')
y = Symbol('y')
#  deuxieme methode pour declarer des variables
h,k = symbols('u,v')
#  troisieme methode pour declarer des variables
var('u,v');


### Manipulations algèbriques¶

In :
(x+y)+(x-y)

Out:
$$2 x$$

#### simplify simplifie des expressions¶

In :
expr=sin(x)**2 + cos(x)**2
expr

Out:
$$\sin^{2}{\left (x \right )} + \cos^{2}{\left (x \right )}$$
In :
simplify(expr)

Out:
$$1$$
In :
expr=(x**3 + x**2 - x - 1)/(x**2 + 2*x + 1)
expr

Out:
$$\frac{x^{3} + x^{2} - x - 1}{x^{2} + 2 x + 1}$$
In :
simplify(expr) # ou expr.simplify()

Out:
$$x - 1$$

#### cancel simplifie des fractions de polynômes¶

In :
expr=(x**3-y**3)/(x**2-y**2)
expr

Out:
$$\frac{x^{3} - y^{3}}{x^{2} - y^{2}}$$
In :
cancel(expr)

Out:
$$\frac{x^{2} + x y + y^{2}}{x + y}$$

#### expand dévéloppe une expression:¶

In :
expr=(x+y)**2
expr

Out:
$$\left(x + y\right)^{2}$$
In :
expand(expr) # ou ((x+y)**2).expand()

Out:
$$x^{2} + 2 x y + y^{2}$$

#### expr.subs(variable,valeur) remplacer dans expression la variable par valeur¶

In :
(expr).subs(x,1)

Out:
$$\left(y + 1\right)^{2}$$
In :
expr = x*x + x*y + x*y + y*y
expr.subs({x:1, y:2})

Out:
$$9$$
In :
expr.subs({x:1-y})
_.simplify()

Out:
$$9$$

#### apart(expr,x) décompose une fraction en éléments simples¶

In :
expr=(x**2+8*x+4)/(x**2-4)
expr

Out:
$$\frac{x^{2} + 8 x + 4}{x^{2} - 4}$$
In :
apart( expr ,x )

Out:
$$1 + \frac{2}{x + 2} + \frac{6}{x - 2}$$

#### Pour remettre tout ensemble, on utilise together(expr, x) :¶

In :
expr=-1/(x+2)+1/(x+1)
expr

Out:
$$- \frac{1}{x + 2} + \frac{1}{x + 1}$$
In :
together(expr ,x)

Out:
$$\frac{1}{\left(x + 1\right) \left(x + 2\right)}$$

#### factor factorise un polynôme¶

In :
expr=x**4/2+5*x**3/12-x**2/3
expr

Out:
$$\frac{x^{4}}{2} + \frac{5 x^{3}}{12} - \frac{x^{2}}{3}$$
In :
factor(expr)

Out:
$$\frac{x^{2}}{12} \left(2 x - 1\right) \left(3 x + 4\right)$$

### solve pour la résolution d'(in)équations et systèmes¶

#### Solution d'une équation¶

On utilise la commande solve pour résoudre une équation:

In :
eq=x**4-1
eq

Out:
$$x^{4} - 1$$
In :
solve(eq,x)

Out:
$$\left [ -1, \quad 1, \quad - i, \quad i\right ]$$
In :
eq=log(x)-x
eq

Out:
$$- x + \log{\left (x \right )}$$
In :
solve(eq)

Out:
$$\left [ - \operatorname{LambertW}{\left (-1 \right )}\right ]$$
In :
var('a,b,c')
expr = a*x**2 + b*x + c
sol=solve(expr, x)
sol

Out:
$$\left [ \frac{1}{2 a} \left(- b + \sqrt{- 4 a c + b^{2}}\right), \quad - \frac{1}{2 a} \left(b + \sqrt{- 4 a c + b^{2}}\right)\right ]$$

#### Solution d'une inéquation¶

Pour resoudre une inéquation on utilise solve_univariate_inequality

In :
solve_univariate_inequality(x**2 > 4, x)

Out:
$$\left(-\infty < x \wedge x < -2\right) \vee \left(2 < x \wedge x < \infty\right)$$

#### Solution d'un système¶

On peut résoudre des systèmes (linéaires ou non linéaires):

In :
x, y, z = symbols('x, y, z')
eqs=[3*x-z+2, y-3*z-5, x+y-z]
eqs

Out:
$$\left [ 3 x - z + 2, \quad y - 3 z - 5, \quad x + y - z\right ]$$
In :
solve(eqs,[x,y,z])

Out:
$$\left \{ x : - \frac{9}{7}, \quad y : - \frac{4}{7}, \quad z : - \frac{13}{7}\right \}$$
In :
eqs=[x**2 + y**2 - 1, x**2 - y**2 - S(1) / 2]
eqs

Out:
$$\left [ x^{2} + y^{2} - 1, \quad x^{2} - y^{2} - \frac{1}{2}\right ]$$
In :
solve(eqs, x, y)

Out:
$$\left [ \left ( - \frac{\sqrt{3}}{2}, \quad - \frac{1}{2}\right ), \quad \left ( - \frac{\sqrt{3}}{2}, \quad \frac{1}{2}\right ), \quad \left ( \frac{\sqrt{3}}{2}, \quad - \frac{1}{2}\right ), \quad \left ( \frac{\sqrt{3}}{2}, \quad \frac{1}{2}\right )\right ]$$

Notons l'utilisation de S() dans la definition de l'équation. En effet, si on écrit $\frac{1}{2}$ directement comme 1/2, Python traduira cette écriture comme $0.5$. Pour que Python considère vraiment la fraction, on peut par exemple forcer le numerateur a être un entier en utilisant la fonction S().

NB Il faut toujours avoir un expris critique même lors de calculs formels ! Voici un exemple. Soit le système linéaire (avec paramètre): $$\begin{cases} x-\alpha y=1,\\ \alpha x-y=1. \end{cases}$$ Par la méthode de Gauss ($L_2 \leftarrow L_2-\alpha L_1$) on obtient $$\begin{cases} x-\alpha y=1,\\ (\alpha^2-1)y= 1-\alpha. \end{cases}$$ Comme $\alpha^2-1=(\alpha-1)(\alpha+1)$ on conclut que

• si $\alpha=1$ (la dernière équation correspond à $0=0$) alors le système possède une infinité de solutions,
• si $\alpha=-1$ (la dernière équation correspond à $0=2$) alors le système ne possède aucune solution,
• si $\alpha\not\in\{-1;1\}$ alors le système possède une solution unique $x=\frac{1}{\alpha+1}$ et $y=-\frac{1}{\alpha+1}$.

Que nous dit sympy?

In :
x,y,alpha=symbols('x,y,alpha')
solve([x-alpha*y-1,alpha*x-y-1],[x,y])

Out:
$$\left \{ x : \frac{1}{\alpha + 1}, \quad y : - \frac{1}{\alpha + 1}\right \}$$
In :
expr1 = 2*x + 3*y - 6
expr2 = 3*x + 2*y - 12
solnliste=solve((expr1, expr2), dict=True)
solnliste
# (dict=True) specifies that we want the result to be returned as a list of Python dictionaries.

Out:
$$\left [ \left \{ x : \frac{24}{5}, \quad y : - \frac{6}{5}\right \}\right ]$$
In :
soln = solnliste # Le dictionnaire est dans une liste !
expr1.subs({x:soln[x], y:soln[y]})

Out:
$$0$$

### Algèbre linéaire¶

On définit une matrice

In :
M = Matrix(((1,1), (2,-1)))
M

Out:
$$\left[\begin{matrix}1 & 1\\2 & -1\end{matrix}\right]$$
In :
det(M)

Out:
$$-3$$

On calcule les valeurs propres et vecteurs propres

In :
P, D = M.diagonalize()

In :
P

Out:
$$\left[\begin{matrix}- \frac{\sqrt{3}}{2} + \frac{1}{2} & \frac{1}{2} + \frac{\sqrt{3}}{2}\\1 & 1\end{matrix}\right]$$
In :
D

Out:
$$\left[\begin{matrix}- \sqrt{3} & 0\\0 & \sqrt{3}\end{matrix}\right]$$
In :
simplify(P*D*P**-1)

Out:
$$\left[\begin{matrix}1 & 1\\2 & -1\end{matrix}\right]$$

Polynôme caractéristique

In :
var('t')
p=M.charpoly(t)
factor(p)

Out:
$$t^{2} - 3$$

### Calcul de limites¶

Pour calculer une limite, la syntaxe est limit(function, variable, point). Par exemple, pour calculer $\lim_{x\to0}f(x)$ il suffit d’entrer limit(f, x, 0) (ou f.limit(x,0)):

In :
limit(cos(x),x,0)

Out:
$$1$$
In :
limit(x,x,oo)

Out:
$$\infty$$
In :
limit(1/x,x,oo)

Out:
$$0$$
In :
limit(x**x,x,0)

Out:
$$1$$

On peut créer une limite sans l'évaluer avec Limit puis l'valuer avec doit():

In :
lim=Limit(1/x,x,oo)
lim

Out:
$$\lim_{x \to \infty} \frac{1}{x}$$
In :
lim.doit()

Out:
$$0$$

On peut aussi indiquer la direction (dir="-" pour la gauche, dir="+" pour la droite)

In :
limit(1/x,x,0,dir="-")

Out:
$$-\infty$$
In :
limit(1/x,x,0,dir="+")

Out:
$$\infty$$

### Graphe d'une fonction $x\mapsto f(x)$¶

In :
plot(sin(x)/x); In :
from sympy.plotting import plot_parametric
plot_parametric(cos(u), sin(u), (u, -5, 5)); In :
from sympy.plotting import plot3d
plot3d(x*y, (x, -5, 5), (y, -5, 5)); ### Calcul de dérivées¶

Pour calculer des dérivées, la syntaxe est diff(f,x) :

In :
diff(sin(x),x)

Out:
$$\cos{\left (x \right )}$$
In :
diff(sin(2*x),x)

Out:
$$2 \cos{\left (2 x \right )}$$

Vérifions ce résultat :

In :
display(Limit((sin(2*(x+h))-sin(2*x))/h,h,0))
limit((sin(2*(x+h))-sin(2*x))/h,h,0)

$$\lim_{u \to 0^+}\left(\frac{1}{u} \left(- \sin{\left (2 x \right )} + \sin{\left (2 u + 2 x \right )}\right)\right)$$
Out:
$$2 \cos{\left (2 x \right )}$$

Dérivée partielle

In :
diff(x*y**2,y)

Out:
$$2 x y$$

Dérivée d'ordre $n$:

In :
diff(x**4,x,3)

Out:
$$24 x$$

On peut créer un dérivée sans l'évaluer avec Derivative puis l'évaluer avec doit():

In :
deriv=Derivative(y*x**4,x,3)
display(deriv)
deriv.doit()

$$\frac{\partial^{3}}{\partial x^{3}} \left(x^{4} y\right)$$
Out:
$$24 x y$$

### Calcul de primitives et intégrales¶

Pour calculer une primitive on utilise la synatxe integrate(f,x):

In :
integrate(sin(x),x)

Out:
$$- \cos{\left (x \right )}$$
In :
integrate(exp(-x**2), x)

Out:
$$\frac{\sqrt{\pi}}{2} \operatorname{erf}{\left (x \right )}$$

On peut créer une primitive sans l'évaluer avec Integral puis l'évaluer avec doit():

In :
primit=Integral(sin(x),x)
display(primit)
primit.doit()

$$\int \sin{\left (x \right )}\, dx$$
Out:
$$- \cos{\left (x \right )}$$

Pour calculer l'intégrale définie $\int_a^b f(x) \mathrm{d}x$ on utilise la syntaxe integrate(f,(x,a,b))

In :
integrate(sin(x),(x,0,pi/2))

Out:
$$1$$

On peut aussi calculer une intégrale impropre :

In :
integrate(exp(-x),(x,0,oo))

Out:
$$1$$

### Développements en série¶

Pour calculer un dévéloppement en série la syntaxe est series(fonction,variable,point,ordre) :

In :
series(cos(x),x,0,5)

Out:
$$1 - \frac{x^{2}}{2} + \frac{x^{4}}{24} + \mathcal{O}\left(x^{5}\right)$$
In :
series(sqrt(1+x),x,0,5)

Out:
$$1 + \frac{x}{2} - \frac{x^{2}}{8} + \frac{x^{3}}{16} - \frac{5 x^{4}}{128} + \mathcal{O}\left(x^{5}\right)$$
In :
series(2*x+sqrt(x**2-1),x,oo,1)

Out:
$$3 x + \mathcal{O}\left(\frac{1}{x}; x\rightarrow \infty\right)$$

### Interpolation¶

Pour calculer le polynôme d'interpolation d'une liste de point on utilise la commande interpolate(liste,variable) :

In :
var('a,y_a,b,y_b,x')
pts=[(a,y_a),(b,y_b)]
interpolate(pts,x)

Out:
$$- \frac{a y_{b}}{- a + b} - \frac{b y_{a}}{a - b} + \frac{x y_{a}}{a - b} + \frac{x y_{b}}{- a + b}$$
In :
pts=[(0,1),(1,2),(2,0)]
interpolate(pts,x)

Out:
$$- \frac{3 x^{2}}{2} + \frac{5 x}{2} + 1$$
In :
pts=[(-a,0),(0,y_b),(a,0)]
interpolate(pts,x).collect(y_b)

Out:
$$y_{b} \left(1 - \frac{x^{2}}{a^{2}}\right)$$

### EDO¶

https://docs.sympy.org/latest/modules/solvers/ode.html

Considérons le problème de Cauchy \begin{cases} u'(x)=-3x^2u(x)+6x^2,\\ u(0)=2. \end{cases}

In :
x=Symbol('x')
C1=Symbol('C1')
u=Function('u')
f=6*x**2-3*x**2*u(x)
edo=Eq(diff(u(x),x),f)
edo

Out:
$$\frac{d}{d x} u{\left (x \right )} = - 3 x^{2} u{\left (x \right )} + 6 x^{2}$$
In :
classify_ode(edo)

Out:
('separable',
'1st_exact',
'1st_linear',
'Bernoulli',
'almost_linear',
'1st_power_series',
'lie_group',
'separable_Integral',
'1st_exact_Integral',
'1st_linear_Integral',
'Bernoulli_Integral',
'almost_linear_Integral')
In :
solgen=dsolve(edo,u(x))
solgen

Out:
$$u{\left (x \right )} = C_{1} e^{- x^{3}} + 2$$

Prise en compte des conditions initiales:

In :
x0=0
u0=4
consts = solve( [solgen.rhs.subs(x,x0)-u0 ], dict=True)
consts

Out:
$$\left \{ C_{1} : 2\right \}$$
In :
solpar=solgen.subs(consts)
solpar

Out:
$$u{\left (x \right )} = 2 + 2 e^{- x^{3}}$$

On transforme le résultat en une fonction:

In :
func = lambdify(x,solpar.rhs,'numpy')

In :
from matplotlib.pylab import *
xx=linspace(0,3,101)
yy=func(xx)
plot(xx,yy); ## Références¶

In [ ]: