|
norm, num_jac, EPS, warn_extraneous)
# Butcher tableau. A is not used directly, see below.
# Eigendecomposition of A is done: A = T L T**-1. There is 1 real eigenvalue # and a complex conjugate pair. They are written below. - 0.5j * (3 ** (5 / 6) + 3 ** (7 / 6)))
# These are transformation matrices. [0.09443876248897524, -0.14125529502095421, 0.03002919410514742], [0.25021312296533332, 0.20412935229379994, -0.38294211275726192], [1, 1, 0]]) [4.17871859155190428, 0.32768282076106237, 0.52337644549944951], [-4.17871859155190428, -0.32768282076106237, 0.47662355450055044], [0.50287263494578682, -2.57192694985560522, 0.59603920482822492]]) # These linear combinations are used in the algorithm.
# Interpolator coefficients. [13/3 + 7*S6/3, -23/3 - 22*S6/3, 10/3 + 5 * S6], [13/3 - 7*S6/3, -23/3 + 22*S6/3, 10/3 - 5 * S6], [1/3, -8/3, 10/3]])
LU_real, LU_complex, solve_lu): """Solve the collocation system.
Parameters ---------- fun : callable Right-hand side of the system. t : float Current time. y : ndarray, shape (n,) Current state. h : float Step to try. Z0 : ndarray, shape (3, n) Initial guess for the solution. It determines new values of `y` at ``t + h * C`` as ``y + Z0``, where ``C`` is the Radau method constants. scale : float Problem tolerance scale, i.e. ``rtol * abs(y) + atol``. tol : float Tolerance to which solve the system. This value is compared with the normalized by `scale` error. LU_real, LU_complex LU decompositions of the system Jacobians. solve_lu : callable Callable which solves a linear system given a LU decomposition. The signature is ``solve_lu(LU, b)``.
Returns ------- converged : bool Whether iterations converged. n_iter : int Number of completed iterations. Z : ndarray, shape (3, n) Found solution. rate : float The rate of convergence. """ n = y.shape[0] M_real = MU_REAL / h M_complex = MU_COMPLEX / h
W = TI.dot(Z0) Z = Z0
F = np.empty((3, n)) ch = h * C
dW_norm_old = None dW = np.empty_like(W) converged = False for k in range(NEWTON_MAXITER): for i in range(3): F[i] = fun(t + ch[i], y + Z[i])
if not np.all(np.isfinite(F)): break
f_real = F.T.dot(TI_REAL) - M_real * W[0] f_complex = F.T.dot(TI_COMPLEX) - M_complex * (W[1] + 1j * W[2])
dW_real = solve_lu(LU_real, f_real) dW_complex = solve_lu(LU_complex, f_complex)
dW[0] = dW_real dW[1] = dW_complex.real dW[2] = dW_complex.imag
dW_norm = norm(dW / scale) if dW_norm_old is not None: rate = dW_norm / dW_norm_old else: rate = None
if (rate is not None and (rate >= 1 or rate ** (NEWTON_MAXITER - k) / (1 - rate) * dW_norm > tol)): break
W += dW Z = T.dot(W)
if (dW_norm == 0 or rate is not None and rate / (1 - rate) * dW_norm < tol): converged = True break
dW_norm_old = dW_norm
return converged, k + 1, Z, rate
"""Predict by which factor to increase/decrease the step size.
The algorithm is described in [1]_.
Parameters ---------- h_abs, h_abs_old : float Current and previous values of the step size, `h_abs_old` can be None (see Notes). error_norm, error_norm_old : float Current and previous values of the error norm, `error_norm_old` can be None (see Notes).
Returns ------- factor : float Predicted factor.
Notes ----- If `h_abs_old` and `error_norm_old` are both not None then a two-step algorithm is used, otherwise a one-step algorithm is used.
References ---------- .. [1] E. Hairer, S. P. Norsett G. Wanner, "Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems", Sec. IV.8. """ if error_norm_old is None or h_abs_old is None or error_norm == 0: multiplier = 1 else: multiplier = h_abs / h_abs_old * (error_norm_old / error_norm) ** 0.25
with np.errstate(divide='ignore'): factor = min(1, multiplier) * error_norm ** -0.25
return factor
"""Implicit Runge-Kutta method of Radau IIA family of order 5.
The implementation follows [1]_. The error is controlled with a third-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 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. 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 [2]_. 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] E. Hairer, G. Wanner, "Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems", Sec. IV.8. .. [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. """ rtol=1e-3, atol=1e-6, jac=None, jac_sparsity=None, vectorized=False, **extraneous): warn_extraneous(extraneous) super(Radau, self).__init__(fun, t0, y0, t_bound, vectorized) self.y_old = None self.max_step = validate_max_step(max_step) self.rtol, self.atol = validate_tol(rtol, atol, self.n) self.f = self.fun(self.t, self.y) # Select initial step assuming the same order which is used to control # the error. self.h_abs = select_initial_step( self.fun, self.t, self.y, self.f, self.direction, 3, self.rtol, self.atol) self.h_abs_old = None self.error_norm_old = None
self.newton_tol = max(10 * EPS / rtol, min(0.03, rtol ** 0.5)) self.sol = None
self.jac_factor = None self.jac, self.J = self._validate_jac(jac, jac_sparsity) if issparse(self.J): def lu(A): self.nlu += 1 return splu(A)
def solve_lu(LU, b): return LU.solve(b)
I = eye(self.n, format='csc') else: def lu(A): self.nlu += 1 return lu_factor(A, overwrite_a=True)
def solve_lu(LU, b): return lu_solve(LU, b, overwrite_b=True)
I = np.identity(self.n)
self.lu = lu self.solve_lu = solve_lu self.I = I
self.current_jac = True self.LU_real = None self.LU_complex = None self.Z = None
t0 = self.t y0 = self.y
if jac is None: if sparsity is not None: if issparse(sparsity): sparsity = csc_matrix(sparsity) groups = group_columns(sparsity) sparsity = (sparsity, groups)
def jac_wrapped(t, y, f): self.njev += 1 J, self.jac_factor = num_jac(self.fun_vectorized, t, y, f, self.atol, self.jac_factor, sparsity) return J J = jac_wrapped(t0, y0, self.f) elif callable(jac): J = jac(t0, y0) self.njev = 1 if issparse(J): J = csc_matrix(J)
def jac_wrapped(t, y, _=None): self.njev += 1 return csc_matrix(jac(t, y), dtype=float)
else: J = np.asarray(J, dtype=float)
def jac_wrapped(t, y, _=None): self.njev += 1 return np.asarray(jac(t, y), dtype=float)
if J.shape != (self.n, self.n): raise ValueError("`jac` is expected to have shape {}, but " "actually has {}." .format((self.n, self.n), J.shape)) else: if issparse(jac): J = csc_matrix(jac) else: J = np.asarray(jac, dtype=float)
if J.shape != (self.n, self.n): raise ValueError("`jac` is expected to have shape {}, but " "actually has {}." .format((self.n, self.n), J.shape)) jac_wrapped = None
return jac_wrapped, J
t = self.t y = self.y f = self.f
max_step = self.max_step atol = self.atol rtol = self.rtol
min_step = 10 * np.abs(np.nextafter(t, self.direction * np.inf) - t) if self.h_abs > max_step: h_abs = max_step h_abs_old = None error_norm_old = None elif self.h_abs < min_step: h_abs = min_step h_abs_old = None error_norm_old = None else: h_abs = self.h_abs h_abs_old = self.h_abs_old error_norm_old = self.error_norm_old
J = self.J LU_real = self.LU_real LU_complex = self.LU_complex
current_jac = self.current_jac jac = self.jac
rejected = False step_accepted = False message = None 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
h = t_new - t h_abs = np.abs(h)
if self.sol is None: Z0 = np.zeros((3, y.shape[0])) else: Z0 = self.sol(t + h * C).T - y
scale = atol + np.abs(y) * rtol
converged = False while not converged: if LU_real is None or LU_complex is None: LU_real = self.lu(MU_REAL / h * self.I - J) LU_complex = self.lu(MU_COMPLEX / h * self.I - J)
converged, n_iter, Z, rate = solve_collocation_system( self.fun, t, y, h, Z0, scale, self.newton_tol, LU_real, LU_complex, self.solve_lu)
if not converged: if current_jac: break
J = self.jac(t, y, f) current_jac = True LU_real = None LU_complex = None
if not converged: h_abs *= 0.5 LU_real = None LU_complex = None continue
y_new = y + Z[-1] ZE = Z.T.dot(E) / h error = self.solve_lu(LU_real, f + ZE) scale = atol + np.maximum(np.abs(y), np.abs(y_new)) * rtol error_norm = norm(error / scale) safety = 0.9 * (2 * NEWTON_MAXITER + 1) / (2 * NEWTON_MAXITER + n_iter)
if rejected and error_norm > 1: error = self.solve_lu(LU_real, self.fun(t, y + error) + ZE) error_norm = norm(error / scale)
if error_norm > 1: factor = predict_factor(h_abs, h_abs_old, error_norm, error_norm_old) h_abs *= max(MIN_FACTOR, safety * factor)
LU_real = None LU_complex = None rejected = True else: step_accepted = True
recompute_jac = jac is not None and n_iter > 2 and rate > 1e-3
factor = predict_factor(h_abs, h_abs_old, error_norm, error_norm_old) factor = min(MAX_FACTOR, safety * factor)
if not recompute_jac and factor < 1.2: factor = 1 else: LU_real = None LU_complex = None
f_new = self.fun(t_new, y_new) if recompute_jac: J = jac(t_new, y_new, f_new) current_jac = True elif jac is not None: current_jac = False
self.h_abs_old = self.h_abs self.error_norm_old = error_norm
self.h_abs = h_abs * factor
self.y_old = y
self.t = t_new self.y = y_new self.f = f_new
self.Z = Z
self.LU_real = LU_real self.LU_complex = LU_complex self.current_jac = current_jac self.J = J
self.t_old = t self.sol = self._compute_dense_output()
return step_accepted, message
Q = np.dot(self.Z.T, P) return RadauDenseOutput(self.t_old, self.t, self.y_old, Q)
return self.sol
super(RadauDenseOutput, self).__init__(t_old, t) self.h = t - t_old self.Q = Q self.order = Q.shape[1] - 1 self.y_old = y_old
x = (t - self.t_old) / self.h if t.ndim == 0: p = np.tile(x, self.order + 1) p = np.cumprod(p) else: p = np.tile(x, (self.order + 1, 1)) p = np.cumprod(p, axis=0) # Here we don't multiply by h, not a mistake. y = np.dot(self.Q, p) if y.ndim == 2: y += self.y_old[:, None] else: y += self.y_old
return y |