Hide keyboard shortcuts

Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1import numpy as np 

2from scipy.integrate import ode 

3from .common import validate_tol, validate_first_step, warn_extraneous 

4from .base import OdeSolver, DenseOutput 

5 

6 

7class LSODA(OdeSolver): 

8 """Adams/BDF method with automatic stiffness detection and switching. 

9 

10 This is a wrapper to the Fortran solver from ODEPACK [1]_. It switches 

11 automatically between the nonstiff Adams method and the stiff BDF method. 

12 The method was originally detailed in [2]_. 

13 

14 Parameters 

15 ---------- 

16 fun : callable 

17 Right-hand side of the system. The calling signature is ``fun(t, y)``. 

18 Here ``t`` is a scalar, and there are two options for the ndarray ``y``: 

19 It can either have shape (n,); then ``fun`` must return array_like with 

20 shape (n,). Alternatively it can have shape (n, k); then ``fun`` 

21 must return an array_like with shape (n, k), i.e. each column 

22 corresponds to a single column in ``y``. The choice between the two 

23 options is determined by `vectorized` argument (see below). The 

24 vectorized implementation allows a faster approximation of the Jacobian 

25 by finite differences (required for this solver). 

26 t0 : float 

27 Initial time. 

28 y0 : array_like, shape (n,) 

29 Initial state. 

30 t_bound : float 

31 Boundary time - the integration won't continue beyond it. It also 

32 determines the direction of the integration. 

33 first_step : float or None, optional 

34 Initial step size. Default is ``None`` which means that the algorithm 

35 should choose. 

36 min_step : float, optional 

37 Minimum allowed step size. Default is 0.0, i.e., the step size is not 

38 bounded and determined solely by the solver. 

39 max_step : float, optional 

40 Maximum allowed step size. Default is np.inf, i.e., the step size is not 

41 bounded and determined solely by the solver. 

42 rtol, atol : float and array_like, optional 

43 Relative and absolute tolerances. The solver keeps the local error 

44 estimates less than ``atol + rtol * abs(y)``. Here `rtol` controls a 

45 relative accuracy (number of correct digits). But if a component of `y` 

46 is approximately below `atol`, the error only needs to fall within 

47 the same `atol` threshold, and the number of correct digits is not 

48 guaranteed. If components of y have different scales, it might be 

49 beneficial to set different `atol` values for different components by 

50 passing array_like with shape (n,) for `atol`. Default values are 

51 1e-3 for `rtol` and 1e-6 for `atol`. 

52 jac : None or callable, optional 

53 Jacobian matrix of the right-hand side of the system with respect to 

54 ``y``. The Jacobian matrix has shape (n, n) and its element (i, j) is 

55 equal to ``d f_i / d y_j``. The function will be called as 

56 ``jac(t, y)``. If None (default), the Jacobian will be 

57 approximated by finite differences. It is generally recommended to 

58 provide the Jacobian rather than relying on a finite-difference 

59 approximation. 

60 lband, uband : int or None 

61 Parameters defining the bandwidth of the Jacobian, 

62 i.e., ``jac[i, j] != 0 only for i - lband <= j <= i + uband``. Setting 

63 these requires your jac routine to return the Jacobian in the packed format: 

64 the returned array must have ``n`` columns and ``uband + lband + 1`` 

65 rows in which Jacobian diagonals are written. Specifically 

66 ``jac_packed[uband + i - j , j] = jac[i, j]``. The same format is used 

67 in `scipy.linalg.solve_banded` (check for an illustration). 

68 These parameters can be also used with ``jac=None`` to reduce the 

69 number of Jacobian elements estimated by finite differences. 

70 vectorized : bool, optional 

71 Whether `fun` is implemented in a vectorized fashion. A vectorized 

72 implementation offers no advantages for this solver. Default is False. 

73 

74 Attributes 

75 ---------- 

76 n : int 

77 Number of equations. 

78 status : string 

79 Current status of the solver: 'running', 'finished' or 'failed'. 

80 t_bound : float 

81 Boundary time. 

82 direction : float 

83 Integration direction: +1 or -1. 

84 t : float 

85 Current time. 

86 y : ndarray 

87 Current state. 

88 t_old : float 

89 Previous time. None if no steps were made yet. 

90 nfev : int 

91 Number of evaluations of the right-hand side. 

92 njev : int 

93 Number of evaluations of the Jacobian. 

94 

95 References 

96 ---------- 

97 .. [1] A. C. Hindmarsh, "ODEPACK, A Systematized Collection of ODE 

98 Solvers," IMACS Transactions on Scientific Computation, Vol 1., 

99 pp. 55-64, 1983. 

100 .. [2] L. Petzold, "Automatic selection of methods for solving stiff and 

101 nonstiff systems of ordinary differential equations", SIAM Journal 

102 on Scientific and Statistical Computing, Vol. 4, No. 1, pp. 136-148, 

103 1983. 

104 """ 

105 def __init__(self, fun, t0, y0, t_bound, first_step=None, min_step=0.0, 

106 max_step=np.inf, rtol=1e-3, atol=1e-6, jac=None, lband=None, 

107 uband=None, vectorized=False, **extraneous): 

108 warn_extraneous(extraneous) 

109 super(LSODA, self).__init__(fun, t0, y0, t_bound, vectorized) 

110 

111 if first_step is None: 

112 first_step = 0 # LSODA value for automatic selection. 

113 else: 

114 first_step = validate_first_step(first_step, t0, t_bound) 

115 

116 first_step *= self.direction 

117 

118 if max_step == np.inf: 

119 max_step = 0 # LSODA value for infinity. 

120 elif max_step <= 0: 

121 raise ValueError("`max_step` must be positive.") 

122 

123 if min_step < 0: 

124 raise ValueError("`min_step` must be nonnegative.") 

125 

126 rtol, atol = validate_tol(rtol, atol, self.n) 

127 

128 solver = ode(self.fun, jac) 

129 solver.set_integrator('lsoda', rtol=rtol, atol=atol, max_step=max_step, 

130 min_step=min_step, first_step=first_step, 

131 lband=lband, uband=uband) 

132 solver.set_initial_value(y0, t0) 

133 

134 # Inject t_bound into rwork array as needed for itask=5. 

135 solver._integrator.rwork[0] = self.t_bound 

136 solver._integrator.call_args[4] = solver._integrator.rwork 

137 

138 self._lsoda_solver = solver 

139 

140 def _step_impl(self): 

141 solver = self._lsoda_solver 

142 integrator = solver._integrator 

143 

144 # From lsoda.step and lsoda.integrate itask=5 means take a single 

145 # step and do not go past t_bound. 

146 itask = integrator.call_args[2] 

147 integrator.call_args[2] = 5 

148 solver._y, solver.t = integrator.run( 

149 solver.f, solver.jac or (lambda: None), solver._y, solver.t, 

150 self.t_bound, solver.f_params, solver.jac_params) 

151 integrator.call_args[2] = itask 

152 

153 if solver.successful(): 

154 self.t = solver.t 

155 self.y = solver._y 

156 # From LSODA Fortran source njev is equal to nlu. 

157 self.njev = integrator.iwork[12] 

158 self.nlu = integrator.iwork[12] 

159 return True, None 

160 else: 

161 return False, 'Unexpected istate in LSODA.' 

162 

163 def _dense_output_impl(self): 

164 iwork = self._lsoda_solver._integrator.iwork 

165 rwork = self._lsoda_solver._integrator.rwork 

166 

167 order = iwork[14] 

168 h = rwork[11] 

169 yh = np.reshape(rwork[20:20 + (order + 1) * self.n], 

170 (self.n, order + 1), order='F').copy() 

171 

172 return LsodaDenseOutput(self.t_old, self.t, h, order, yh) 

173 

174 

175class LsodaDenseOutput(DenseOutput): 

176 def __init__(self, t_old, t, h, order, yh): 

177 super(LsodaDenseOutput, self).__init__(t_old, t) 

178 self.h = h 

179 self.yh = yh 

180 self.p = np.arange(order + 1) 

181 

182 def _call_impl(self, t): 

183 if t.ndim == 0: 

184 x = ((t - self.t) / self.h) ** self.p 

185 else: 

186 x = ((t - self.t) / self.h) ** self.p[:, None] 

187 

188 return np.dot(self.yh, x)