class_docstr: 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 [1]_.
  This class implements a quasi-constant step size as explained in [2]_.
  The error estimation strategy for the constant-step BDF is derived in [3]_.
  An accuracy enhancement using modified formulas (NDF) [2]_ is also implemented.
  
  Can be applied in the 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 the ndarray ``y``:
      It can either have shape (n,); then ``fun`` must return array_like with
      shape (n,). Alternatively it can have shape (n, k); then ``fun``
      must return an 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 a faster approximation of the Jacobian
      by finite differences (required for this solver).
  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.
  max_step : float, optional
      Maximum allowed step size. Default is np.inf, i.e. the step size 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`, 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 this method. The Jacobian matrix has shape (n, n) and its
      element (i, j) is equal to ``d f_i / d y_j``.
      There are three ways to define the Jacobian:
  
          * If array_like or sparse_matrix, the Jacobian is assumed to
            be constant.
          * If callable, the Jacobian is assumed to depend on both
            t and y; it will be called as ``jac(t, y)`` as necessary.
            For the 'Radau' and 'BDF' methods, the return value might be a
            sparse matrix.
          * If None (default), 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). This argument
      is ignored if `jac` is not `None`. If the Jacobian has only few non-zero
      elements in *each* row, providing the sparsity structure will greatly
      speed up the computations [4]_. A zero entry means that a corresponding
      element in the Jacobian is always 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 evaluations of the right-hand side.
  njev : int
      Number of evaluations of the Jacobian.
  nlu : int
      Number of LU decompositions.
  
  References
  ----------
  .. [1] 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.
  .. [2] L. F. Shampine, M. W. Reichelt, "THE MATLAB ODE SUITE", SIAM J. SCI.
         COMPUTE., Vol. 18, No. 1, pp. 1-22, January 1997.
  .. [3] E. Hairer, G. Wanner, "Solving Ordinary Differential Equations I:
         Nonstiff Problems", Sec. III.2.
  .. [4] 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.
class_docstr_tok: ["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","[","1","]","_",".","This","class","implements","a","quasi-constant","step","size","as","explained","in","[","2","]","_",".","The","error","estimation","strategy","for","the","constant-step","BDF","is","derived","in","[","3","]","_",".","An","accuracy","enhancement","using","modified","formulas","(","NDF",")","[","2","]","_","is","also","implemented",".","Can","be","applied","in","the","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","the","ndarray","``","y","``",":","It","can","either","have","shape","(","n",",",")",";","then","``","fun","``","must","return","array_like","with","shape","(","n",",",")",".","Alternatively","it","can","have","shape","(","n",",","k",")",";","then","``","fun","``","must","return","an","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","a","faster","approximation","of","the","Jacobian","by","finite","differences","(","required","for","this","solver",")",".","t0",":","float","Initial","time",".","y0",":","array_like",",","shape","(","n",",",")","Initial","state",".","t_bound",":","float","Boundary","time","-","the","integration","wo","n'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",".","max_step",":","float",",","optional","Maximum","allowed","step","size",".","Default","is","np.inf",",","i.e",".","the","step","size","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","`",",","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","this","method",".","The","Jacobian","matrix","has","shape","(","n",",","n",")","and","its","element","(","i",",","j",")","is","equal","to","``","d","f_i","\/","d","y_j","``",".","There","are","three","ways","to","define","the","Jacobian",":","*","If","array_like","or","sparse_matrix",",","the","Jacobian","is","assumed","to","be","constant",".","*","If","callable",",","the","Jacobian","is","assumed","to","depend","on","both","t","and","y",";","it","will","be","called","as","``","jac","(","t",",","y",")","``","as","necessary",".","For","the","'Radau","'","and","'BDF","'","methods",",","the","return","value","might","be","a","sparse","matrix",".","*","If","None","(","default",")",",","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",")",".","This","argument","is","ignored","if","`","jac","`","is","not","`","None","`",".","If","the","Jacobian","has","only","few","non-zero","elements","in","*","each","*","row",",","providing","the","sparsity","structure","will","greatly","speed","up","the","computations","[","4","]","_",".","A","zero","entry","means","that","a","corresponding","element","in","the","Jacobian","is","always","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","evaluations","of","the","right-hand","side",".","njev",":","int","Number","of","evaluations","of","the","Jacobian",".","nlu",":","int","Number","of","LU","decompositions",".","References","--","--","--","--","--","..","[","1","]","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",".","..","[","2","]","L.","F.","Shampine",",","M.","W.","Reichelt",",","``","THE","MATLAB","ODE","SUITE","''",",","SIAM","J.","SCI",".","COMPUTE.",",","Vol",".","18",",","No",".","1",",","pp",".","1-22",",","January","1997",".","..","[","3","]","E.","Hairer",",","G.","Wanner",",","``","Solving","Ordinary","Differential","Equations","I",":","Nonstiff","Problems","''",",","Sec",".","III.2",".","..","[","4","]","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","."]
code: def _step_impl(self):
      t = self.t
      D = self.D
  
      max_step = self.max_step
      min_step = 10 * np.abs(
          np.nextafter(t, self.direction * np.inf) - t
      )
      if self.h_abs > max_step:
          h_abs = max_step
          change_D(D, self.order, max_step / self.h_abs)
          self.n_equal_steps = 0
      elif self.h_abs < min_step:
          h_abs = min_step
          change_D(D, self.order, min_step / self.h_abs)
          self.n_equal_steps = 0
      else:
          h_abs = self.h_abs
  
      atol = self.atol
      rtol = self.rtol
      order = self.order
  
      alpha = self.alpha
      gamma = self.gamma
      error_const = self.error_const
  
      J = self.J
      LU = self.LU
      current_jac = self.jac is None
  
      step_accepted = False
      while not step_accepted:
          if h_abs < min_step:
              return False, self.TOO_SMALL_STEP
  
          h = h_abs * self.direction
          t_new = t + h
  
          if self.direction * (t_new - self.t_bound) > 0:
              t_new = self.t_bound
              change_D(D, order, np.abs(t_new - t) / h_abs)
              self.n_equal_steps = 0
              LU = None
  
          h = t_new - t
          h_abs = np.abs(h)
  
          y_predict = np.sum(D[: order + 1], axis=0)
  
          scale = atol + rtol * np.abs(y_predict)
          psi = (
              np.dot(D[1 : order + 1].T, gamma[1 : order + 1])
              / alpha[order]
          )
  
          converged = False
          c = h / alpha[order]
          while not converged:
              if LU is None:
                  LU = self.lu(self.I - c * J)
  
              converged, n_iter, y_new, d = solve_bdf_system(
                  self.fun,
                  t_new,
                  y_predict,
                  c,
                  psi,
                  LU,
                  self.solve_lu,
                  scale,
                  self.newton_tol,
              )
  
              if not converged:
                  if current_jac:
                      break
                  J = self.jac(t_new, y_predict)
                  LU = None
                  current_jac = True
  
          if not converged:
              factor = 0.5
              h_abs *= factor
              change_D(D, order, factor)
              self.n_equal_steps = 0
              LU = None
              continue
  
          safety = (
              0.9
              * (2 * NEWTON_MAXITER + 1)
              / (2 * NEWTON_MAXITER + n_iter)
          )
  
          scale = atol + rtol * np.abs(y_new)
          error = error_const[order] * d
          error_norm = norm(error / scale)
  
          if error_norm > 1:
              factor = max(
                  MIN_FACTOR, safety * error_norm ** (-1 / (order + 1))
              )
              h_abs *= factor
              change_D(D, order, factor)
              self.n_equal_steps = 0
              # As we didn't have problems with convergence, we don't
              # reset LU here.
          else:
              step_accepted = True
  
      self.n_equal_steps += 1
  
      self.t = t_new
      self.y = y_new
  
      self.h_abs = h_abs
      self.J = J
      self.LU = LU
  
      # Update differences. The principal relation here is
      # D^{j + 1} y_n = D^{j} y_n - D^{j} y_{n - 1}. Keep in mind that D
      # contained difference for previous interpolating polynomial and
      # d = D^{k + 1} y_n. Thus this elegant code follows.
      D[order + 2] = d - D[order + 1]
      D[order + 1] = d
      for i in reversed(range(order + 1)):
          D[i] += D[i + 1]
  
      if self.n_equal_steps < order + 1:
          return True, None
  
      if order > 1:
          error_m = error_const[order - 1] * D[order]
          error_m_norm = norm(error_m / scale)
      else:
          error_m_norm = np.inf
  
      if order < MAX_ORDER:
          error_p = error_const[order + 1] * D[order + 2]
          error_p_norm = norm(error_p / scale)
      else:
          error_p_norm = np.inf
  
      error_norms = np.array([error_m_norm, error_norm, error_p_norm])
      factors = error_norms ** (-1 / np.arange(order, order + 3))
  
      delta_order = np.argmax(factors) - 1
      order += delta_order
      self.order = order
  
      factor = min(MAX_FACTOR, safety * np.max(factors))
      self.h_abs *= factor
      change_D(D, order, factor)
      self.n_equal_steps = 0
      self.LU = None
  
      return True, None
  
imports: from __future__ import division, print_function, absolute_import
  import numpy
  from scipy.linalg import lu_factor, lu_solve
  from scipy.sparse import issparse, csc_matrix, eye
  from scipy.sparse.linalg import splu
  from scipy.optimize._numdiff import group_columns
  from .common import validate_max_step, validate_tol, select_initial_step, norm, EPS, num_jac, validate_first_step, warn_extraneous
  from .base import OdeSolver, DenseOutput
label_desc: Use image node_id 3 for calling the BDF obj's underlying member method code with example usage: obj._step_impl() and returns: True,  None, True,  None, False,  self