Backward differentiation formula (BDF) solver for stiff ODEs.
Implements the solver described in [Shampine and Reichelt (1997)], a variable step size, variable order (VSVO) BDF integrator with order varying between 1 and 5.
Each step involves solving the following nonlinear equation by Newton's method:
0 = 1/1 * BDF(1, y)[n+1] + ... + 1/k * BDF(k, y)[n+1] - h ode_fn(t[n+1], y[n+1]) - bdf_coefficients[k-1] * (1/1 + ... + 1/k) * (y[n+1] - y[n] - BDF(1, y)[n] - ... - BDF(k, y)[n])
k >= 1 is the current order of the integrator,
h is the current step
bdf_coefficients is a list of numbers that parameterizes the method,
BDF(m, y) is the
m-th order backward difference of the vector
BDF(0, y)[n] = y[n] and
BDF(m + 1, y)[n] = BDF(m, y)[n] - BDF(m, y)[n - 1] for
m >= 0.
Newton's method can fail because * the method has exceeded the maximum number of iterations, * the method is converging too slowly, or * the method is not expected to converge (the last two conditions are determined by approximating the Lipschitz constant associated with the iteration).
True, the solver avoids evaluating the
Jacobian of the dynamics function as much as possible. In particular, Newton's
method will try to use the Jacobian from a previous integration step. If
Newton's method fails with an out-of-date Jacobian, the Jacobian is
re-evaluated and Newton's method is restarted. If Newton's method fails and
the Jacobian is already up-to-date, then the step size is decreased and
Newton's method is restarted.
Even if Newton's method converges, the solution it generates can still be rejected if it exceeds the specified error tolerance due to truncation error. In this case, the step size is decreased and Newton's method is restarted.
: Lawrence F. Shampine and Mark W. Reichelt. The MATLAB ODE Suite. SIAM Journal on Scientific Computing 18(1):1-22, 1997.
__init__( rtol=0.001, atol=1e-06, first_step_size=None, safety_factor=0.9, min_step_size_factor=0.1, max_step_size_factor=10.0, max_num_steps=None, max_order=bdf_util.MAX_ORDER, max_num_newton_iters=4, newton_tol_factor=0.1, newton_step_size_factor=0.5, bdf_coefficients=(-0.185, -1.0 / 9.0, -0.0823, -0.0415, 0.0), evaluate_jacobian_lazily=False, use_pfor_to_compute_jacobian=True, validate_args=False, name='bdf' )
Initializes the solver.
rtol: Optional float
Tensorspecifying an upper bound on relative error, per element of the dependent variable. The error tolerance for the next step is
tol = atol + rtol * abs(state)where
stateis the computed state at the current step (see also
atol). The next step is rejected if it incurs a local truncation error larger than
tol. Default value:
atol: Optional float
Tensorspecifying an upper bound on absolute error, per element of the dependent variable (see also
rtol). Default value:
first_step_size: Optional scalar float
Tensorspecifying the size of the first step. If unspecified, the size is chosen automatically. Default value:
safety_factor: Scalar positive float
Tensor. When Newton's method converges, the solver may choose to update the step size by applying a multiplicative factor to the current step size. This factor is
factor = clamp(factor_unclamped, min_step_size_factor, max_step_size_factor)where
factor_unclamped = error_ratio**(-1. / (order + 1)) * safety_factor(see also
max_step_size_factor). A small (respectively, large) value for the safety factor causes the solver to take smaller (respectively, larger) step sizes. A value larger than one, though not explicitly prohibited, is discouraged. Default value:
min_step_size_factor: Scalar float
safety_factor). Default value:
max_step_size_factor: Scalar float
safety_factor). Default value:
max_num_steps: Optional scalar integer
Tensorspecifying the maximum number of steps allowed (including rejected steps). If unspecified, there is no upper bound on the number of steps. Default value:
max_order: Scalar integer
Tensortaking values between 1 and 5 (inclusive) specifying the maximum BDF order. Default value:
max_num_newton_iters: Optional scalar integer
Tensorspecifying the maximum number of iterations per invocation of Newton's method. If unspecified, there is no upper bound on the number iterations. Default value:
newton_tol_factor: Scalar float
Tensorused to determine the stopping condition for Newton's method. In particular, Newton's method terminates when the distance to the root is estimated to be less than
newton_tol_factor * norm(atol + rtol * abs(state))where
stateis the computed state at the current step. Default value:
newton_step_size_factor: Scalar float
Tensorspecifying a multiplicative factor applied to the size of the integration step when Newton's method fails. Default value:
bdf_coefficients: 1-D float
Tensorwith 5 entries that parameterize the solver. The default values are those proposed in . Default value:
(-0.1850, -1. / 9., -0.0823, -0.0415, 0.).
evaluate_jacobian_lazily: Optional boolean specifying whether the Jacobian should be evaluated at each point in time or as needed (i.e., lazily). Default value:
use_pfor_to_compute_jacobian: Boolean specifying whether or not to use parallel for in computing the Jacobian when
jacobian_fnis not specified. Default value:
validate_args: Whether to validate input with asserts. If
Falseand the inputs are invalid, correct behavior is not guaranteed. Default value:
strname prefixed to Ops created by this function. Default value:
solve( ode_fn, initial_time, initial_state, solution_times, jacobian_fn=None, jacobian_sparsity=None, batch_ndims=None, previous_solver_internal_state=None )