Coverage for /usr/lib/python3/dist-packages/scipy/integrate/

1import numpy as np

2from scipy.linalg import lu_factor, lu_solve

3from scipy.sparse import issparse, csc_matrix, eye

4from scipy.sparse.linalg import splu

5from scipy.optimize._numdiff import group_columns

6from .common import (validate_max_step, validate_tol, select_initial_step,

7 norm, EPS, num_jac, validate_first_step,

8 warn_extraneous)

9from .base import OdeSolver, DenseOutput

12MAX_ORDER = 5

13NEWTON_MAXITER = 4

14MIN_FACTOR = 0.2

15MAX_FACTOR = 10

18def compute_R(order, factor):

19 """Compute the matrix for changing the differences array."""

20 I = np.arange(1, order + 1)[:, None]

21 J = np.arange(1, order + 1)

22 M = np.zeros((order + 1, order + 1))

23 M[1:, 1:] = (I - 1 - factor * J) / I

24 M[0] = 1

25 return np.cumprod(M, axis=0)

28def change_D(D, order, factor):

29 """Change differences array in-place when step size is changed."""

30 R = compute_R(order, factor)

31 U = compute_R(order, 1)

32 RU = R.dot(U)

33 D[:order + 1] = np.dot(RU.T, D[:order + 1])

36def solve_bdf_system(fun, t_new, y_predict, c, psi, LU, solve_lu, scale, tol):

37 """Solve the algebraic system resulting from BDF method."""

38 d = 0

39 y = y_predict.copy()

40 dy_norm_old = None

41 converged = False

42 for k in range(NEWTON_MAXITER):

43 f = fun(t_new, y)

44 if not np.all(np.isfinite(f)):

45 break

47 dy = solve_lu(LU, c * f - psi - d)

48 dy_norm = norm(dy / scale)

50 if dy_norm_old is None:

51 rate = None

52 else:

53 rate = dy_norm / dy_norm_old

55 if (rate is not None and (rate >= 1 or

56 rate ** (NEWTON_MAXITER - k) / (1 - rate) * dy_norm > tol)):

57 break

59 y += dy

60 d += dy

62 if (dy_norm == 0 or

63 rate is not None and rate / (1 - rate) * dy_norm < tol):

64 converged = True

65 break

67 dy_norm_old = dy_norm

69 return converged, k + 1, y, d

72class BDF(OdeSolver):

73 """Implicit method based on backward-differentiation formulas.

75 This is a variable order method with the order varying automatically from

76 1 to 5. The general framework of the BDF algorithm is described in [1]_.

77 This class implements a quasi-constant step size as explained in [2]_.

78 The error estimation strategy for the constant-step BDF is derived in [3]_.

79 An accuracy enhancement using modified formulas (NDF) [2]_ is also implemented.

81 Can be applied in the complex domain.

83 Parameters

84 ----------

85 fun : callable

86 Right-hand side of the system: the time derivative of the state ``y``

87 at time ``t``. The calling signature is ``fun(t, y)``, where ``t`` is a

88 scalar and ``y`` is an ndarray with ``len(y) = len(y0)``. ``fun`` must

89 return an array of the same shape as ``y``. See `vectorized` for more

90 information.

91 t0 : float

92 Initial time.

93 y0 : array_like, shape (n,)

94 Initial state.

95 t_bound : float

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

97 determines the direction of the integration.

98 first_step : float or None, optional

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

100 should choose.

101 max_step : float, optional

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

103 bounded and determined solely by the solver.

104 rtol, atol : float and array_like, optional

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

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

107 relative accuracy (number of correct digits), while `atol` controls

108 absolute accuracy (number of correct decimal places). To achieve the

109 desired `rtol`, set `atol` to be smaller than the smallest value that

110 can be expected from ``rtol * abs(y)`` so that `rtol` dominates the

111 allowable error. If `atol` is larger than ``rtol * abs(y)`` the

112 number of correct digits is not guaranteed. Conversely, to achieve the

113 desired `atol` set `rtol` such that ``rtol * abs(y)`` is always smaller

114 than `atol`. If components of y have different scales, it might be

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

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

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

118 jac : {None, array_like, sparse_matrix, callable}, optional

119 Jacobian matrix of the right-hand side of the system with respect to y,

120 required by this method. The Jacobian matrix has shape (n, n) and its

121 element (i, j) is equal to ``d f_i / d y_j``.

122 There are three ways to define the Jacobian:

123

124 * If array_like or sparse_matrix, the Jacobian is assumed to

125 be constant.

126 * If callable, the Jacobian is assumed to depend on both

127 t and y; it will be called as ``jac(t, y)`` as necessary.

128 For the 'Radau' and 'BDF' methods, the return value might be a

129 sparse matrix.

130 * If None (default), the Jacobian will be approximated by

131 finite differences.

132

133 It is generally recommended to provide the Jacobian rather than

134 relying on a finite-difference approximation.

135 jac_sparsity : {None, array_like, sparse matrix}, optional

136 Defines a sparsity structure of the Jacobian matrix for a

137 finite-difference approximation. Its shape must be (n, n). This argument

138 is ignored if `jac` is not `None`. If the Jacobian has only few non-zero

139 elements in *each* row, providing the sparsity structure will greatly

140 speed up the computations [4]_. A zero entry means that a corresponding

141 element in the Jacobian is always zero. If None (default), the Jacobian

142 is assumed to be dense.

143 vectorized : bool, optional

144 Whether `fun` can be called in a vectorized fashion. Default is False.

145

146 If ``vectorized`` is False, `fun` will always be called with ``y`` of

147 shape ``(n,)``, where ``n = len(y0)``.

148

149 If ``vectorized`` is True, `fun` may be called with ``y`` of shape

150 ``(n, k)``, where ``k`` is an integer. In this case, `fun` must behave

151 such that ``fun(t, y)[:, i] == fun(t, y[:, i])`` (i.e. each column of

152 the returned array is the time derivative of the state corresponding

153 with a column of ``y``).

154

155 Setting ``vectorized=True`` allows for faster finite difference

156 approximation of the Jacobian by this method, but may result in slower

157 execution overall in some circumstances (e.g. small ``len(y0)``).

158

159 Attributes

160 ----------

161 n : int

162 Number of equations.

163 status : string

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

165 t_bound : float

166 Boundary time.

167 direction : float

168 Integration direction: +1 or -1.

169 t : float

170 Current time.

171 y : ndarray

172 Current state.

173 t_old : float

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

175 step_size : float

176 Size of the last successful step. None if no steps were made yet.

177 nfev : int

178 Number of evaluations of the right-hand side.

179 njev : int

180 Number of evaluations of the Jacobian.

181 nlu : int

182 Number of LU decompositions.

183

184 References

185 ----------

186 .. [1] G. D. Byrne, A. C. Hindmarsh, "A Polyalgorithm for the Numerical

187 Solution of Ordinary Differential Equations", ACM Transactions on

188 Mathematical Software, Vol. 1, No. 1, pp. 71-96, March 1975.

189 .. [2] L. F. Shampine, M. W. Reichelt, "THE MATLAB ODE SUITE", SIAM J. SCI.

190 COMPUTE., Vol. 18, No. 1, pp. 1-22, January 1997.

191 .. [3] E. Hairer, G. Wanner, "Solving Ordinary Differential Equations I:

192 Nonstiff Problems", Sec. III.2.

193 .. [4] A. Curtis, M. J. D. Powell, and J. Reid, "On the estimation of

194 sparse Jacobian matrices", Journal of the Institute of Mathematics

195 and its Applications, 13, pp. 117-120, 1974.

196 """

197 def __init__(self, fun, t0, y0, t_bound, max_step=np.inf,

198 rtol=1e-3, atol=1e-6, jac=None, jac_sparsity=None,

199 vectorized=False, first_step=None, **extraneous):

200 warn_extraneous(extraneous)

201 super().__init__(fun, t0, y0, t_bound, vectorized,

202 support_complex=True)

203 self.max_step = validate_max_step(max_step)

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

205 f = self.fun(self.t, self.y)

206 if first_step is None:

207 self.h_abs = select_initial_step(self.fun, self.t, self.y, f,

208 self.direction, 1,

209 self.rtol, self.atol)

210 else:

211 self.h_abs = validate_first_step(first_step, t0, t_bound)

212 self.h_abs_old = None

213 self.error_norm_old = None

214

215 self.newton_tol = max(10 * EPS / rtol, min(0.03, rtol ** 0.5))

216

217 self.jac_factor = None

218 self.jac, self.J = self._validate_jac(jac, jac_sparsity)

219 if issparse(self.J):

220 def lu(A):

221 self.nlu += 1

222 return splu(A)

223

224 def solve_lu(LU, b):

225 return LU.solve(b)

226

227 I = eye(self.n, format='csc', dtype=self.y.dtype)

228 else:

229 def lu(A):

230 self.nlu += 1

231 return lu_factor(A, overwrite_a=True)

232

233 def solve_lu(LU, b):

234 return lu_solve(LU, b, overwrite_b=True)

235

236 I = np.identity(self.n, dtype=self.y.dtype)

237

238 self.lu = lu

239 self.solve_lu = solve_lu

240 self.I = I

241

242 kappa = np.array([0, -0.1850, -1/9, -0.0823, -0.0415, 0])

243 self.gamma = np.hstack((0, np.cumsum(1 / np.arange(1, MAX_ORDER + 1))))

244 self.alpha = (1 - kappa) * self.gamma

245 self.error_const = kappa * self.gamma + 1 / np.arange(1, MAX_ORDER + 2)

246

247 D = np.empty((MAX_ORDER + 3, self.n), dtype=self.y.dtype)

248 D[0] = self.y

249 D[1] = f * self.h_abs * self.direction

250 self.D = D

251

252 self.order = 1

253 self.n_equal_steps = 0

254 self.LU = None

255

256 def _validate_jac(self, jac, sparsity):

257 t0 = self.t

258 y0 = self.y

259

260 if jac is None:

261 if sparsity is not None:

262 if issparse(sparsity):

263 sparsity = csc_matrix(sparsity)

264 groups = group_columns(sparsity)

265 sparsity = (sparsity, groups)

266

267 def jac_wrapped(t, y):

268 self.njev += 1

269 f = self.fun_single(t, y)

270 J, self.jac_factor = num_jac(self.fun_vectorized, t, y, f,

271 self.atol, self.jac_factor,

272 sparsity)

273 return J

274 J = jac_wrapped(t0, y0)

275 elif callable(jac):

276 J = jac(t0, y0)

277 self.njev += 1

278 if issparse(J):

279 J = csc_matrix(J, dtype=y0.dtype)

280

281 def jac_wrapped(t, y):

282 self.njev += 1

283 return csc_matrix(jac(t, y), dtype=y0.dtype)

284 else:

285 J = np.asarray(J, dtype=y0.dtype)

286

287 def jac_wrapped(t, y):

288 self.njev += 1

289 return np.asarray(jac(t, y), dtype=y0.dtype)

290

291 if J.shape != (self.n, self.n):

292 raise ValueError("`jac` is expected to have shape {}, but "

293 "actually has {}."

294 .format((self.n, self.n), J.shape))

295 else:

296 if issparse(jac):

297 J = csc_matrix(jac, dtype=y0.dtype)

298 else:

299 J = np.asarray(jac, dtype=y0.dtype)

300

301 if J.shape != (self.n, self.n):

302 raise ValueError("`jac` is expected to have shape {}, but "

303 "actually has {}."

304 .format((self.n, self.n), J.shape))

305 jac_wrapped = None

306

307 return jac_wrapped, J

308

309 def _step_impl(self):

310 t = self.t

311 D = self.D

312

313 max_step = self.max_step

314 min_step = 10 * np.abs(np.nextafter(t, self.direction * np.inf) - t)

315 if self.h_abs > max_step:

316 h_abs = max_step

317 change_D(D, self.order, max_step / self.h_abs)

318 self.n_equal_steps = 0

319 elif self.h_abs < min_step:

320 h_abs = min_step

321 change_D(D, self.order, min_step / self.h_abs)

322 self.n_equal_steps = 0

323 else:

324 h_abs = self.h_abs

325

326 atol = self.atol

327 rtol = self.rtol

328 order = self.order

329

330 alpha = self.alpha

331 gamma = self.gamma

332 error_const = self.error_const

333

334 J = self.J

335 LU = self.LU

336 current_jac = self.jac is None

337

338 step_accepted = False

339 while not step_accepted:

340 if h_abs < min_step:

341 return False, self.TOO_SMALL_STEP

342

343 h = h_abs * self.direction

344 t_new = t + h

345

346 if self.direction * (t_new - self.t_bound) > 0:

347 t_new = self.t_bound

348 change_D(D, order, np.abs(t_new - t) / h_abs)

349 self.n_equal_steps = 0

350 LU = None

351

352 h = t_new - t

353 h_abs = np.abs(h)

354

355 y_predict = np.sum(D[:order + 1], axis=0)

356

357 scale = atol + rtol * np.abs(y_predict)

358 psi = np.dot(D[1: order + 1].T, gamma[1: order + 1]) / alpha[order]

359

360 converged = False

361 c = h / alpha[order]

362 while not converged:

363 if LU is None:

364 LU = self.lu(self.I - c * J)

365

366 converged, n_iter, y_new, d = solve_bdf_system(

367 self.fun, t_new, y_predict, c, psi, LU, self.solve_lu,

368 scale, self.newton_tol)

369

370 if not converged:

371 if current_jac:

372 break

373 J = self.jac(t_new, y_predict)

374 LU = None

375 current_jac = True

376

377 if not converged:

378 factor = 0.5

379 h_abs *= factor

380 change_D(D, order, factor)

381 self.n_equal_steps = 0

382 LU = None

383 continue

384

385 safety = 0.9 * (2 * NEWTON_MAXITER + 1) / (2 * NEWTON_MAXITER

386 + n_iter)

387

388 scale = atol + rtol * np.abs(y_new)

389 error = error_const[order] * d

390 error_norm = norm(error / scale)

391

392 if error_norm > 1:

393 factor = max(MIN_FACTOR,

394 safety * error_norm ** (-1 / (order + 1)))

395 h_abs *= factor

396 change_D(D, order, factor)

397 self.n_equal_steps = 0

398 # As we didn't have problems with convergence, we don't

399 # reset LU here.

400 else:

401 step_accepted = True

402

403 self.n_equal_steps += 1

404

405 self.t = t_new

406 self.y = y_new

407

408 self.h_abs = h_abs

409 self.J = J

410 self.LU = LU

411

412 # Update differences. The principal relation here is

413 # D^{j + 1} y_n = D^{j} y_n - D^{j} y_{n - 1}. Keep in mind that D

414 # contained difference for previous interpolating polynomial and

415 # d = D^{k + 1} y_n. Thus this elegant code follows.

416 D[order + 2] = d - D[order + 1]

417 D[order + 1] = d

418 for i in reversed(range(order + 1)):

419 D[i] += D[i + 1]

420

421 if self.n_equal_steps < order + 1:

422 return True, None

423

424 if order > 1:

425 error_m = error_const[order - 1] * D[order]

426 error_m_norm = norm(error_m / scale)

427 else:

428 error_m_norm = np.inf

429

430 if order < MAX_ORDER:

431 error_p = error_const[order + 1] * D[order + 2]

432 error_p_norm = norm(error_p / scale)

433 else:

434 error_p_norm = np.inf

435

436 error_norms = np.array([error_m_norm, error_norm, error_p_norm])

437 with np.errstate(divide='ignore'):

438 factors = error_norms ** (-1 / np.arange(order, order + 3))

439

440 delta_order = np.argmax(factors) - 1

441 order += delta_order

442 self.order = order

443

444 factor = min(MAX_FACTOR, safety * np.max(factors))

445 self.h_abs *= factor

446 change_D(D, order, factor)

447 self.n_equal_steps = 0

448 self.LU = None

449

450 return True, None

451

452 def _dense_output_impl(self):

453 return BdfDenseOutput(self.t_old, self.t, self.h_abs * self.direction,

454 self.order, self.D[:self.order + 1].copy())

455

456

457class BdfDenseOutput(DenseOutput):

458 def __init__(self, t_old, t, h, order, D):

459 super().__init__(t_old, t)

460 self.order = order

461 self.t_shift = self.t - h * np.arange(self.order)

462 self.denom = h * (1 + np.arange(self.order))

463 self.D = D

464

465 def _call_impl(self, t):

466 if t.ndim == 0:

467 x = (t - self.t_shift) / self.denom

468 p = np.cumprod(x)

469 else:

470 x = (t - self.t_shift[:, None]) / self.denom[:, None]

471 p = np.cumprod(x, axis=0)

472

473 y = np.dot(self.D[1:].T, p)

474 if y.ndim == 1:

475 y += self.D[0]

476 else:

477 y += self.D[0, :, None]

478

479 return y

Coverage for /usr/lib/python3/dist-packages/scipy/integrate/_ivp/bdf.py: 9%

244 statements