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 csc_matrix, issparse, 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, num_jac, EPS, warn_extraneous,

8 validate_first_step)

9from .base import OdeSolver, DenseOutput

11S6 = 6 ** 0.5

13# Butcher tableau. A is not used directly, see below.

14C = np.array([(4 - S6) / 10, (4 + S6) / 10, 1])

15E = np.array([-13 - 7 * S6, -13 + 7 * S6, -1]) / 3

17# Eigendecomposition of A is done: A = T L T**-1. There is 1 real eigenvalue

18# and a complex conjugate pair. They are written below.

19MU_REAL = 3 + 3 ** (2 / 3) - 3 ** (1 / 3)

20MU_COMPLEX = (3 + 0.5 * (3 ** (1 / 3) - 3 ** (2 / 3))

21 - 0.5j * (3 ** (5 / 6) + 3 ** (7 / 6)))

23# These are transformation matrices.

24T = np.array([

25 [0.09443876248897524, -0.14125529502095421, 0.03002919410514742],

26 [0.25021312296533332, 0.20412935229379994, -0.38294211275726192],

27 [1, 1, 0]])

28TI = np.array([

29 [4.17871859155190428, 0.32768282076106237, 0.52337644549944951],

30 [-4.17871859155190428, -0.32768282076106237, 0.47662355450055044],

31 [0.50287263494578682, -2.57192694985560522, 0.59603920482822492]])

32# These linear combinations are used in the algorithm.

33TI_REAL = TI[0]

34TI_COMPLEX = TI[1] + 1j * TI[2]

36# Interpolator coefficients.

37P = np.array([

38 [13/3 + 7*S6/3, -23/3 - 22*S6/3, 10/3 + 5 * S6],

39 [13/3 - 7*S6/3, -23/3 + 22*S6/3, 10/3 - 5 * S6],

40 [1/3, -8/3, 10/3]])

43NEWTON_MAXITER = 6 # Maximum number of Newton iterations.

44MIN_FACTOR = 0.2 # Minimum allowed decrease in a step size.

45MAX_FACTOR = 10 # Maximum allowed increase in a step size.

48def solve_collocation_system(fun, t, y, h, Z0, scale, tol,

49 LU_real, LU_complex, solve_lu):

50 """Solve the collocation system.

52 Parameters

53 ----------

54 fun : callable

55 Right-hand side of the system.

56 t : float

57 Current time.

58 y : ndarray, shape (n,)

59 Current state.

60 h : float

61 Step to try.

62 Z0 : ndarray, shape (3, n)

63 Initial guess for the solution. It determines new values of `y` at

64 ``t + h * C`` as ``y + Z0``, where ``C`` is the Radau method constants.

65 scale : ndarray, shape (n)

66 Problem tolerance scale, i.e. ``rtol * abs(y) + atol``.

67 tol : float

68 Tolerance to which solve the system. This value is compared with

69 the normalized by `scale` error.

70 LU_real, LU_complex

71 LU decompositions of the system Jacobians.

72 solve_lu : callable

73 Callable which solves a linear system given a LU decomposition. The

74 signature is ``solve_lu(LU, b)``.

76 Returns

77 -------

78 converged : bool

79 Whether iterations converged.

80 n_iter : int

81 Number of completed iterations.

82 Z : ndarray, shape (3, n)

83 Found solution.

84 rate : float

85 The rate of convergence.

86 """

87 n = y.shape[0]

88 M_real = MU_REAL / h

89 M_complex = MU_COMPLEX / h

91 W = TI.dot(Z0)

92 Z = Z0

94 F = np.empty((3, n))

95 ch = h * C

97 dW_norm_old = None

98 dW = np.empty_like(W)

99 converged = False

100 rate = None

101 for k in range(NEWTON_MAXITER):

102 for i in range(3):

103 F[i] = fun(t + ch[i], y + Z[i])

104

105 if not np.all(np.isfinite(F)):

106 break

107

108 f_real = F.T.dot(TI_REAL) - M_real * W[0]

109 f_complex = F.T.dot(TI_COMPLEX) - M_complex * (W[1] + 1j * W[2])

110

111 dW_real = solve_lu(LU_real, f_real)

112 dW_complex = solve_lu(LU_complex, f_complex)

113

114 dW[0] = dW_real

115 dW[1] = dW_complex.real

116 dW[2] = dW_complex.imag

117

118 dW_norm = norm(dW / scale)

119 if dW_norm_old is not None:

120 rate = dW_norm / dW_norm_old

121

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

123 rate ** (NEWTON_MAXITER - k) / (1 - rate) * dW_norm > tol)):

124 break

125

126 W += dW

127 Z = T.dot(W)

128

129 if (dW_norm == 0 or

130 rate is not None and rate / (1 - rate) * dW_norm < tol):

131 converged = True

132 break

133

134 dW_norm_old = dW_norm

135

136 return converged, k + 1, Z, rate

137

138

139def predict_factor(h_abs, h_abs_old, error_norm, error_norm_old):

140 """Predict by which factor to increase/decrease the step size.

141

142 The algorithm is described in [1]_.

143

144 Parameters

145 ----------

146 h_abs, h_abs_old : float

147 Current and previous values of the step size, `h_abs_old` can be None

148 (see Notes).

149 error_norm, error_norm_old : float

150 Current and previous values of the error norm, `error_norm_old` can

151 be None (see Notes).

152

153 Returns

154 -------

155 factor : float

156 Predicted factor.

157

158 Notes

159 -----

160 If `h_abs_old` and `error_norm_old` are both not None then a two-step

161 algorithm is used, otherwise a one-step algorithm is used.

162

163 References

164 ----------

165 .. [1] E. Hairer, S. P. Norsett G. Wanner, "Solving Ordinary Differential

166 Equations II: Stiff and Differential-Algebraic Problems", Sec. IV.8.

167 """

168 if error_norm_old is None or h_abs_old is None or error_norm == 0:

169 multiplier = 1

170 else:

171 multiplier = h_abs / h_abs_old * (error_norm_old / error_norm) ** 0.25

172

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

174 factor = min(1, multiplier) * error_norm ** -0.25

175

176 return factor

177

178

179class Radau(OdeSolver):

180 """Implicit Runge-Kutta method of Radau IIA family of order 5.

181

182 The implementation follows [1]_. The error is controlled with a

183 third-order accurate embedded formula. A cubic polynomial which satisfies

184 the collocation conditions is used for the dense output.

185

186 Parameters

187 ----------

188 fun : callable

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

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

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

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

193 information.

194 t0 : float

195 Initial time.

196 y0 : array_like, shape (n,)

197 Initial state.

198 t_bound : float

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

200 determines the direction of the integration.

201 first_step : float or None, optional

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

203 should choose.

204 max_step : float, optional

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

206 bounded and determined solely by the solver.

207 rtol, atol : float and array_like, optional

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

209 estimates less than ``atol + rtol * abs(y)``. HHere `rtol` controls a

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

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

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

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

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

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

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

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

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

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

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

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

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

223 y, required by this method. The Jacobian matrix has shape (n, n) and

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

225 There are three ways to define the Jacobian:

226

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

228 be constant.

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

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

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

232 sparse matrix.

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

234 finite differences.

235

236 It is generally recommended to provide the Jacobian rather than

237 relying on a finite-difference approximation.

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

239 Defines a sparsity structure of the Jacobian matrix for a

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

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

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

243 speed up the computations [2]_. A zero entry means that a corresponding

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

245 is assumed to be dense.

246 vectorized : bool, optional

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

248

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

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

251

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

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

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

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

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

257

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

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

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

261

262 Attributes

263 ----------

264 n : int

265 Number of equations.

266 status : string

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

268 t_bound : float

269 Boundary time.

270 direction : float

271 Integration direction: +1 or -1.

272 t : float

273 Current time.

274 y : ndarray

275 Current state.

276 t_old : float

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

278 step_size : float

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

280 nfev : int

281 Number of evaluations of the right-hand side.

282 njev : int

283 Number of evaluations of the Jacobian.

284 nlu : int

285 Number of LU decompositions.

286

287 References

288 ----------

289 .. [1] E. Hairer, G. Wanner, "Solving Ordinary Differential Equations II:

290 Stiff and Differential-Algebraic Problems", Sec. IV.8.

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

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

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

294 """

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

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

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

298 warn_extraneous(extraneous)

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

300 self.y_old = None

301 self.max_step = validate_max_step(max_step)

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

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

304 # Select initial step assuming the same order which is used to control

305 # the error.

306 if first_step is None:

307 self.h_abs = select_initial_step(

308 self.fun, self.t, self.y, self.f, self.direction,

309 3, self.rtol, self.atol)

310 else:

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

312 self.h_abs_old = None

313 self.error_norm_old = None

314

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

316 self.sol = None

317

318 self.jac_factor = None

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

320 if issparse(self.J):

321 def lu(A):

322 self.nlu += 1

323 return splu(A)

324

325 def solve_lu(LU, b):

326 return LU.solve(b)

327

328 I = eye(self.n, format='csc')

329 else:

330 def lu(A):

331 self.nlu += 1

332 return lu_factor(A, overwrite_a=True)

333

334 def solve_lu(LU, b):

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

336

337 I = np.identity(self.n)

338

339 self.lu = lu

340 self.solve_lu = solve_lu

341 self.I = I

342

343 self.current_jac = True

344 self.LU_real = None

345 self.LU_complex = None

346 self.Z = None

347

348 def _validate_jac(self, jac, sparsity):

349 t0 = self.t

350 y0 = self.y

351

352 if jac is None:

353 if sparsity is not None:

354 if issparse(sparsity):

355 sparsity = csc_matrix(sparsity)

356 groups = group_columns(sparsity)

357 sparsity = (sparsity, groups)

358

359 def jac_wrapped(t, y, f):

360 self.njev += 1

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

362 self.atol, self.jac_factor,

363 sparsity)

364 return J

365 J = jac_wrapped(t0, y0, self.f)

366 elif callable(jac):

367 J = jac(t0, y0)

368 self.njev = 1

369 if issparse(J):

370 J = csc_matrix(J)

371

372 def jac_wrapped(t, y, _=None):

373 self.njev += 1

374 return csc_matrix(jac(t, y), dtype=float)

375

376 else:

377 J = np.asarray(J, dtype=float)

378

379 def jac_wrapped(t, y, _=None):

380 self.njev += 1

381 return np.asarray(jac(t, y), dtype=float)

382

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

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

385 "actually has {}."

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

387 else:

388 if issparse(jac):

389 J = csc_matrix(jac)

390 else:

391 J = np.asarray(jac, dtype=float)

392

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

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

395 "actually has {}."

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

397 jac_wrapped = None

398

399 return jac_wrapped, J

400

401 def _step_impl(self):

402 t = self.t

403 y = self.y

404 f = self.f

405

406 max_step = self.max_step

407 atol = self.atol

408 rtol = self.rtol

409

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

411 if self.h_abs > max_step:

412 h_abs = max_step

413 h_abs_old = None

414 error_norm_old = None

415 elif self.h_abs < min_step:

416 h_abs = min_step

417 h_abs_old = None

418 error_norm_old = None

419 else:

420 h_abs = self.h_abs

421 h_abs_old = self.h_abs_old

422 error_norm_old = self.error_norm_old

423

424 J = self.J

425 LU_real = self.LU_real

426 LU_complex = self.LU_complex

427

428 current_jac = self.current_jac

429 jac = self.jac

430

431 rejected = False

432 step_accepted = False

433 message = None

434 while not step_accepted:

435 if h_abs < min_step:

436 return False, self.TOO_SMALL_STEP

437

438 h = h_abs * self.direction

439 t_new = t + h

440

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

442 t_new = self.t_bound

443

444 h = t_new - t

445 h_abs = np.abs(h)

446

447 if self.sol is None:

448 Z0 = np.zeros((3, y.shape[0]))

449 else:

450 Z0 = self.sol(t + h * C).T - y

451

452 scale = atol + np.abs(y) * rtol

453

454 converged = False

455 while not converged:

456 if LU_real is None or LU_complex is None:

457 LU_real = self.lu(MU_REAL / h * self.I - J)

458 LU_complex = self.lu(MU_COMPLEX / h * self.I - J)

459

460 converged, n_iter, Z, rate = solve_collocation_system(

461 self.fun, t, y, h, Z0, scale, self.newton_tol,

462 LU_real, LU_complex, self.solve_lu)

463

464 if not converged:

465 if current_jac:

466 break

467

468 J = self.jac(t, y, f)

469 current_jac = True

470 LU_real = None

471 LU_complex = None

472

473 if not converged:

474 h_abs *= 0.5

475 LU_real = None

476 LU_complex = None

477 continue

478

479 y_new = y + Z[-1]

480 ZE = Z.T.dot(E) / h

481 error = self.solve_lu(LU_real, f + ZE)

482 scale = atol + np.maximum(np.abs(y), np.abs(y_new)) * rtol

483 error_norm = norm(error / scale)

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

485 + n_iter)

486

487 if rejected and error_norm > 1:

488 error = self.solve_lu(LU_real, self.fun(t, y + error) + ZE)

489 error_norm = norm(error / scale)

490

491 if error_norm > 1:

492 factor = predict_factor(h_abs, h_abs_old,

493 error_norm, error_norm_old)

494 h_abs *= max(MIN_FACTOR, safety * factor)

495

496 LU_real = None

497 LU_complex = None

498 rejected = True

499 else:

500 step_accepted = True

501

502 recompute_jac = jac is not None and n_iter > 2 and rate > 1e-3

503

504 factor = predict_factor(h_abs, h_abs_old, error_norm, error_norm_old)

505 factor = min(MAX_FACTOR, safety * factor)

506

507 if not recompute_jac and factor < 1.2:

508 factor = 1

509 else:

510 LU_real = None

511 LU_complex = None

512

513 f_new = self.fun(t_new, y_new)

514 if recompute_jac:

515 J = jac(t_new, y_new, f_new)

516 current_jac = True

517 elif jac is not None:

518 current_jac = False

519

520 self.h_abs_old = self.h_abs

521 self.error_norm_old = error_norm

522

523 self.h_abs = h_abs * factor

524

525 self.y_old = y

526

527 self.t = t_new

528 self.y = y_new

529 self.f = f_new

530

531 self.Z = Z

532

533 self.LU_real = LU_real

534 self.LU_complex = LU_complex

535 self.current_jac = current_jac

536 self.J = J

537

538 self.t_old = t

539 self.sol = self._compute_dense_output()

540

541 return step_accepted, message

542

543 def _compute_dense_output(self):

544 Q = np.dot(self.Z.T, P)

545 return RadauDenseOutput(self.t_old, self.t, self.y_old, Q)

546

547 def _dense_output_impl(self):

548 return self.sol

549

550

551class RadauDenseOutput(DenseOutput):

552 def __init__(self, t_old, t, y_old, Q):

553 super().__init__(t_old, t)

554 self.h = t - t_old

555 self.Q = Q

556 self.order = Q.shape[1] - 1

557 self.y_old = y_old

558

559 def _call_impl(self, t):

560 x = (t - self.t_old) / self.h

561 if t.ndim == 0:

562 p = np.tile(x, self.order + 1)

563 p = np.cumprod(p)

564 else:

565 p = np.tile(x, (self.order + 1, 1))

566 p = np.cumprod(p, axis=0)

567 # Here we don't multiply by h, not a mistake.

568 y = np.dot(self.Q, p)

569 if y.ndim == 2:

570 y += self.y_old[:, None]

571 else:

572 y += self.y_old

573

574 return y

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

261 statements