A firstorder, fully implicit integrator optimized for secondorder systems, with a secondorder error estimate.
The velocityimplicit Euler integrator is a variant of the firstorder implicit Euler that takes advantage of the simple mapping q̇ = N(q) v of second order systems to formulate a smaller problem in velocities (and miscellaneous states if any) only. For systems with secondorder dynamics, VelocityImplicitEulerIntegrator formulates a problem that is half as large as that formulated by Drake's ImplicitEulerIntegrator, resulting in improved runtime performance. Upon convergence of the resulting system of equations, this method provides the same discretization as ImplicitEulerIntegrator, but at a fraction of the computational cost.
This integrator requires a system of ordinary differential equations (ODEs) in state x = (q,v,z)
to be expressible in the following form:
q̇ = N(q) v; (1) ẏ = f_y(t,q,y), (2)
where q̇
and v
are linearly related via the kinematic mapping N(q)
, y = (v,z)
, and f_y
is a function that can depend on the time and state.
Implicit Euler uses the following update rule at time step n:
qⁿ⁺¹ = qⁿ + h N(qⁿ⁺¹) vⁿ⁺¹; (3) yⁿ⁺¹ = yⁿ + h f_y(tⁿ⁺¹,qⁿ⁺¹,yⁿ⁺¹). (4)
To solve the nonlinear system for (qⁿ⁺¹,yⁿ⁺¹)
, the velocityimplicit Euler integrator iterates with a modified Newton's method: At iteration k
, it finds a (qₖ₊₁,yₖ₊₁)
that attempts to satisfy
qₖ₊₁ = qⁿ + h N(qₖ) vₖ₊₁. (5) yₖ₊₁ = yⁿ + h f_y(tⁿ⁺¹,qₖ₊₁,yₖ₊₁); (6)
In this notation, the n
's index timesteps, while the k
's index the specific NewtonRaphson iterations within each time step.
Notice that we've intentionally lagged N(qₖ) one iteration behind in Eq (5). This allows it to substitute (5) into (6) to obtain a nonlinear system in y
only. Contrast this strategy with the one implemented by ImplicitEulerIntegrator, which solves a larger nonlinear system in the full state x.
To find a (qₖ₊₁,yₖ₊₁)
that approximately satisfies (56), we linearize the system (56) to compute a Newton step. Define
ℓ(y) = f_y(tⁿ⁺¹,qⁿ + h N(qₖ) v,y), (7) Jₗ(y) = ∂ℓ(y) / ∂y. (8)
To advance the Newton step, the velocityimplicit Euler integrator solves the following linear equation for Δy
:
(I  h Jₗ) Δy =  R(yₖ), (9)
where R(y) = y  yⁿ  h ℓ(y)
and Δy = yₖ₊₁  yₖ
. The Δy
solution directly gives us yₖ₊₁
. It then substitutes the vₖ₊₁
component of yₖ₊₁
in (5) to get qₖ₊₁
.
This implementation uses a Newton method and relies upon the convergence to a solution for y
in R(y) = 0
where R(y) = y  yⁿ  h ℓ(y)
as h
becomes sufficiently small. General implementational details for the Newton method were gleaned from Section IV.8 in [Hairer, 1996].
In this integrator, we simultaneously take a large step at the requested step size of h as well as two halfsized steps each with step size h/2
. The result from two halfsized steps is propagated as the solution, while the difference between the two results is used as the error estimate for the propagated solution. This error estimate is accurate to the second order.
To be precise, let x̅ⁿ⁺¹
be the computed solution from a large step, x̃ⁿ⁺¹
be the computed solution from two small steps, and xⁿ⁺¹
be the true solution. Since the integrator propagates x̃ⁿ⁺¹
as its solution, we denote the true error vector as ε = x̃ⁿ⁺¹  xⁿ⁺¹
. VelocityImplicitEulerIntegrator uses ε* = x̅ⁿ⁺¹  x̃ⁿ⁺¹
, the difference between the two solutions, as the secondorder error estimate, because for a smooth system, ‖ε*‖ = O(h²)
, and ‖ε  ε*‖ = O(h³)
. See the notes in VelocityImplicitEulerIntegrator<T>::get_error_estimate_order() for a detailed derivation of the error estimate's truncation error.
In this implementation, VelocityImplicitEulerIntegrator<T> attempts the large fullsized step before attempting the two small halfsized steps, because the large step is more likely to fail to converge, and if it is performed first, convergence failures are detected early, avoiding the unnecessary effort of computing potentiallysuccessful small steps.
h
is the smallest irrevocable timeincrement advanced by this integrator: if, for example, the second small halfsized step fails, this integrator revokes to the state before the first small step. This behavior is similar to other integrators with multistage evaluation: the stepcounting statistics treat a "step" as the combination of all the stages. T  The scalar type, which must be one of the default nonsymbolic scalars. 
#include <drake/systems/analysis/velocity_implicit_euler_integrator.h>
Public Member Functions  
~VelocityImplicitEulerIntegrator () override=default  
VelocityImplicitEulerIntegrator (const System< T > &system, Context< T > *context=nullptr)  
bool  supports_error_estimation () const final 
Returns true, because this integrator supports error estimation. More...  
int  get_error_estimate_order () const final 
Returns the asymptotic order of the difference between the large and small steps (from which the error estimate is computed), which is 2. More...  
Does not allow copy, move, or assignment  
VelocityImplicitEulerIntegrator (const VelocityImplicitEulerIntegrator &)=delete  
VelocityImplicitEulerIntegrator &  operator= (const VelocityImplicitEulerIntegrator &)=delete 
VelocityImplicitEulerIntegrator (VelocityImplicitEulerIntegrator &&)=delete  
VelocityImplicitEulerIntegrator &  operator= (VelocityImplicitEulerIntegrator &&)=delete 
Public Member Functions inherited from ImplicitIntegrator< T >  
virtual  ~ImplicitIntegrator () 
ImplicitIntegrator (const System< T > &system, Context< T > *context=nullptr)  
int  max_newton_raphson_iterations () const 
The maximum number of NewtonRaphson iterations to take before the NewtonRaphson process decides that convergence will not be attained. More...  
void  set_reuse (bool reuse) 
Sets whether the integrator attempts to reuse Jacobian matrices and iteration matrix factorizations (default is true ). More...  
bool  get_reuse () const 
Gets whether the integrator attempts to reuse Jacobian matrices and iteration matrix factorizations. More...  
void  set_use_full_newton (bool flag) 
Sets whether the method operates in "full Newton" mode, in which case Jacobian and iteration matrices are freshly computed on every NewtonRaphson iteration. More...  
bool  get_use_full_newton () const 
Gets whether this method is operating in "full Newton" mode. More...  
void  set_jacobian_computation_scheme (JacobianComputationScheme scheme) 
Sets the Jacobian computation scheme. More...  
JacobianComputationScheme  get_jacobian_computation_scheme () const 
int64_t  get_num_derivative_evaluations_for_jacobian () const 
Gets the number of ODE function evaluations (calls to EvalTimeDerivatives()) used only for computing the Jacobian matrices since the last call to ResetStatistics(). More...  
int64_t  get_num_jacobian_evaluations () const 
Gets the number of Jacobian computations (i.e., the number of times that the Jacobian matrix was reformed) since the last call to ResetStatistics(). More...  
int64_t  get_num_newton_raphson_iterations () const 
Gets the number of iterations used in the NewtonRaphson nonlinear systems of equation solving process since the last call to ResetStatistics(). More...  
int64_t  get_num_iteration_matrix_factorizations () const 
Gets the number of factorizations of the iteration matrix since the last call to ResetStatistics(). More...  
int64_t  get_num_error_estimator_derivative_evaluations () const 
Gets the number of ODE function evaluations (calls to EvalTimeDerivatives()) used only for the error estimation process since the last call to ResetStatistics(). More...  
int64_t  get_num_error_estimator_derivative_evaluations_for_jacobian () const 
int64_t  get_num_error_estimator_newton_raphson_iterations () const 
Gets the number of iterations used in the NewtonRaphson nonlinear systems of equation solving process for the error estimation process since the last call to ResetStatistics(). More...  
int64_t  get_num_error_estimator_jacobian_evaluations () const 
Gets the number of Jacobian matrix computations used only during the error estimation process since the last call to ResetStatistics(). More...  
int64_t  get_num_error_estimator_iteration_matrix_factorizations () const 
Gets the number of factorizations of the iteration matrix used only during the error estimation process since the last call to ResetStatistics(). More...  
Public Member Functions inherited from IntegratorBase< T >  
IntegratorBase (const System< T > &system, Context< T > *context=nullptr)  
Maintains references to the system being integrated and the context used to specify the initial conditions for that system (if any). More...  
virtual  ~IntegratorBase ()=default 
void  Reset () 
Resets the integrator to initial values, i.e., default construction values. More...  
void  Initialize () 
An integrator must be initialized before being used. More...  
StepResult  IntegrateNoFurtherThanTime (const T &publish_time, const T &update_time, const T &boundary_time) 
(Internal use only) Integrates the system forward in time by a single step with step size subject to integration error tolerances (assuming that the integrator supports error estimation). More...  
void  IntegrateWithMultipleStepsToTime (const T &t_final) 
Stepping function for integrators operating outside of Simulator that advances the continuous state exactly to t_final . More...  
bool  IntegrateWithSingleFixedStepToTime (const T &t_target) 
Stepping function for integrators operating outside of Simulator that advances the continuous state using a single step to t_target . More...  
const Context< T > &  get_context () const 
Returns a const reference to the internallymaintained Context holding the most recent state in the trajectory. More...  
Context< T > *  get_mutable_context () 
Returns a mutable pointer to the internallymaintained Context holding the most recent state in the trajectory. More...  
void  reset_context (Context< T > *context) 
Replace the pointer to the internallymaintained Context with a different one. More...  
const System< T > &  get_system () const 
Gets a constant reference to the system that is being integrated (and was provided to the constructor of the integrator). More...  
bool  is_initialized () const 
Indicates whether the integrator has been initialized. More...  
const T &  get_previous_integration_step_size () const 
Gets the size of the last (previous) integration step. More...  
IntegratorBase (const IntegratorBase &)=delete  
IntegratorBase &  operator= (const IntegratorBase &)=delete 
IntegratorBase (IntegratorBase &&)=delete  
IntegratorBase &  operator= (IntegratorBase &&)=delete 
void  set_target_accuracy (double accuracy) 
Request that the integrator attempt to achieve a particular accuracy for the continuous portions of the simulation. More...  
double  get_target_accuracy () const 
Gets the target accuracy. More...  
double  get_accuracy_in_use () const 
Gets the accuracy in use by the integrator. More...  
const ContinuousState< T > *  get_error_estimate () const 
Gets the error estimate (used only for integrators that support error estimation). More...  
const T &  get_ideal_next_step_size () const 
Return the step size the integrator would like to take next, based primarily on the integrator's accuracy prediction. More...  
void  set_fixed_step_mode (bool flag) 
Sets an integrator with error control to fixed step mode. More...  
bool  get_fixed_step_mode () const 
Gets whether an integrator is running in fixed step mode. More...  
const Eigen::VectorXd &  get_generalized_state_weight_vector () const 
Gets the weighting vector (equivalent to a diagonal matrix) applied to weighting both generalized coordinate and velocity state variable errors, as described in the group documentation. More...  
Eigen::VectorBlock< Eigen::VectorXd >  get_mutable_generalized_state_weight_vector () 
Gets a mutable weighting vector (equivalent to a diagonal matrix) applied to weighting both generalized coordinate and velocity state variable errors, as described in the group documentation. More...  
const Eigen::VectorXd &  get_misc_state_weight_vector () const 
Gets the weighting vector (equivalent to a diagonal matrix) for weighting errors in miscellaneous continuous state variables z . More...  
Eigen::VectorBlock< Eigen::VectorXd >  get_mutable_misc_state_weight_vector () 
Gets a mutable weighting vector (equivalent to a diagonal matrix) for weighting errors in miscellaneous continuous state variables z . More...  
void  request_initial_step_size_target (const T &step_size) 
Request that the first attempted integration step have a particular size. More...  
const T &  get_initial_step_size_target () const 
Gets the target size of the first integration step. More...  
void  set_maximum_step_size (const T &max_step_size) 
Sets the maximum step size that may be taken by this integrator. More...  
const T &  get_maximum_step_size () const 
Gets the maximum step size that may be taken by this integrator. More...  
double  get_stretch_factor () const 
Gets the stretch factor (> 1), which is multiplied by the maximum (typically userdesignated) integration step size to obtain the amount that the integrator is able to stretch the maximum time step toward hitting an upcoming publish or update event in IntegrateNoFurtherThanTime(). More...  
void  set_requested_minimum_step_size (const T &min_step_size) 
Sets the requested minimum step size h_min that may be taken by this integrator. More...  
const T &  get_requested_minimum_step_size () const 
Gets the requested minimum step size h_min for this integrator. More...  
void  set_throw_on_minimum_step_size_violation (bool throws) 
Sets whether the integrator should throw a std::exception when the integrator's step size selection algorithm determines that it must take a step smaller than the minimum step size (for, e.g., purposes of error control). More...  
bool  get_throw_on_minimum_step_size_violation () const 
Reports the current setting of the throw_on_minimum_step_size_violation flag. More...  
T  get_working_minimum_step_size () const 
Gets the current value of the working minimum step size h_work(t) for this integrator, which may vary with the current time t as stored in the integrator's context. More...  
void  ResetStatistics () 
Forget accumulated statistics. More...  
int64_t  get_num_substep_failures () const 
Gets the number of failed substeps (implying one or more step size reductions was required to permit solving the necessary nonlinear system of equations). More...  
int64_t  get_num_step_shrinkages_from_substep_failures () const 
Gets the number of step size shrinkages due to substep failures (e.g., integrator convergence failures) since the last call to ResetStatistics() or Initialize(). More...  
int64_t  get_num_step_shrinkages_from_error_control () const 
Gets the number of step size shrinkages due to failure to meet targeted error tolerances, since the last call to ResetStatistics or Initialize(). More...  
int64_t  get_num_derivative_evaluations () const 
Returns the number of ODE function evaluations (calls to CalcTimeDerivatives()) since the last call to ResetStatistics() or Initialize(). More...  
const T &  get_actual_initial_step_size_taken () const 
The actual size of the successful first step. More...  
const T &  get_smallest_adapted_step_size_taken () const 
The size of the smallest step taken as the result of a controlled integration step adjustment since the last Initialize() or ResetStatistics() call. More...  
const T &  get_largest_step_size_taken () const 
The size of the largest step taken since the last Initialize() or ResetStatistics() call. More...  
int64_t  get_num_steps_taken () const 
The number of integration steps taken since the last Initialize() or ResetStatistics() call. More...  
void  add_derivative_evaluations (double evals) 
Manually increments the statistic for the number of ODE evaluations. More...  
void  StartDenseIntegration () 
Starts dense integration, allocating a new dense output for this integrator to use. More...  
const trajectories::PiecewisePolynomial< T > *  get_dense_output () const 
Returns a const pointer to the integrator's current PiecewisePolynomial instance, holding a representation of the continuous state trajectory since the last StartDenseIntegration() call. More...  
std::unique_ptr< trajectories::PiecewisePolynomial< T > >  StopDenseIntegration () 
Stops dense integration, yielding ownership of the current dense output to the caller. More...  
Additional Inherited Members  
Public Types inherited from ImplicitIntegrator< T >  
enum  JacobianComputationScheme { kForwardDifference, kCentralDifference, kAutomatic } 
Public Types inherited from IntegratorBase< T >  
enum  StepResult { kReachedPublishTime = 1, kReachedZeroCrossing = 2, kReachedUpdateTime = 3, kTimeHasAdvanced = 4, kReachedBoundaryTime = 5, kReachedStepLimit = 6 } 
Status returned by IntegrateNoFurtherThanTime(). More...  
Protected Types inherited from ImplicitIntegrator< T >  
enum  ConvergenceStatus { kDiverged, kConverged, kNotConverged } 
Protected Member Functions inherited from ImplicitIntegrator< T >  
virtual int  do_max_newton_raphson_iterations () const 
Derived classes can override this method to change the number of NewtonRaphson iterations (10 by default) to take before the NewtonRaphson process decides that convergence will not be attained. More...  
bool  MaybeFreshenMatrices (const T &t, const VectorX< T > &xt, const T &h, int trial, const std::function< void(const MatrixX< T > &J, const T &h, typename ImplicitIntegrator< T >::IterationMatrix *)> &compute_and_factor_iteration_matrix, typename ImplicitIntegrator< T >::IterationMatrix *iteration_matrix) 
Computes necessary matrices (Jacobian and iteration matrix) for NewtonRaphson (NR) iterations, as necessary. More...  
void  FreshenMatricesIfFullNewton (const T &t, const VectorX< T > &xt, const T &h, const std::function< void(const MatrixX< T > &J, const T &h, typename ImplicitIntegrator< T >::IterationMatrix *)> &compute_and_factor_iteration_matrix, typename ImplicitIntegrator< T >::IterationMatrix *iteration_matrix) 
Computes necessary matrices (Jacobian and iteration matrix) for full NewtonRaphson (NR) iterations, if full NewtonRaphson method is activated (if it's not activated, this method is a noop). More...  
bool  IsUpdateZero (const VectorX< T > &xc, const VectorX< T > &dxc, double eps=1.0) const 
Checks whether a proposed update is effectively zero, indicating that the NewtonRaphson process converged. More...  
ConvergenceStatus  CheckNewtonConvergence (int iteration, const VectorX< T > &xtplus, const VectorX< T > &dx, const T &dx_norm, const T &last_dx_norm) const 
Checks a NewtonRaphson iteration process for convergence. More...  
virtual void  DoImplicitIntegratorReset () 
Derived classes can override this method to perform routines when Reset() is called. More...  
bool  IsBadJacobian (const MatrixX< T > &J) const 
Checks to see whether a Jacobian matrix is "bad" (has any NaN or Inf values) and needs to be recomputed. More...  
MatrixX< T > &  get_mutable_jacobian () 
void  DoResetStatistics () override 
Resets any statistics particular to a specific integrator. More...  
void  DoReset () final 
Derived classes can override this method to perform routines when Reset() is called. More...  
const MatrixX< T > &  CalcJacobian (const T &t, const VectorX< T > &x) 
void  ComputeForwardDiffJacobian (const System< T > &system, const T &t, const VectorX< T > &xt, Context< T > *context, MatrixX< T > *J) 
void  ComputeCentralDiffJacobian (const System< T > &system, const T &t, const VectorX< T > &xt, Context< T > *context, MatrixX< T > *J) 
void  ComputeAutoDiffJacobian (const System< T > &system, const T &t, const VectorX< T > &xt, const Context< T > &context, MatrixX< T > *J) 
void  increment_num_iter_factorizations () 
void  increment_jacobian_computation_derivative_evaluations (int count) 
void  increment_jacobian_evaluations () 
void  set_jacobian_is_fresh (bool flag) 
template<>  
void  ComputeAutoDiffJacobian (const System< AutoDiffXd > &, const AutoDiffXd &, const VectorX< AutoDiffXd > &, const Context< AutoDiffXd > &, MatrixX< AutoDiffXd > *) 
Protected Member Functions inherited from IntegratorBase< T >  
const ContinuousState< T > &  EvalTimeDerivatives (const Context< T > &context) 
Evaluates the derivative function and updates call statistics. More...  
template<typename U >  
const ContinuousState< U > &  EvalTimeDerivatives (const System< U > &system, const Context< U > &context) 
Evaluates the derivative function (and updates call statistics). More...  
void  set_accuracy_in_use (double accuracy) 
Sets the working ("in use") accuracy for this integrator. More...  
bool  StepOnceErrorControlledAtMost (const T &h_max) 
Default code for advancing the continuous state of the system by a single step of h_max (or smaller, depending on error control). More...  
T  CalcStateChangeNorm (const ContinuousState< T > &dx_state) const 
Computes the infinity norm of a change in continuous state. More...  
std::pair< bool, T >  CalcAdjustedStepSize (const T &err, const T &attempted_step_size, bool *at_minimum_step_size) const 
Calculates adjusted integrator step sizes toward keeping state variables within error bounds on the next integration step. More...  
trajectories::PiecewisePolynomial< T > *  get_mutable_dense_output () 
Returns a mutable pointer to the internallymaintained PiecewisePolynomial instance, holding a representation of the continuous state trajectory since the last time StartDenseIntegration() was called. More...  
bool  DoDenseStep (const T &h) 
Calls DoStep(h) while recording the resulting step in the dense output. More...  
ContinuousState< T > *  get_mutable_error_estimate () 
Gets an error estimate of the state variables recorded by the last call to StepOnceFixedSize(). More...  
void  set_actual_initial_step_size_taken (const T &h) 
void  set_smallest_adapted_step_size_taken (const T &h) 
Sets the size of the smalleststeptaken statistic as the result of a controlled integration step adjustment. More...  
void  set_largest_step_size_taken (const T &h) 
void  set_ideal_next_step_size (const T &h) 

delete 

delete 

overridedefault 

explicit 

finalvirtual 
Returns the asymptotic order of the difference between the large and small steps (from which the error estimate is computed), which is 2.
That is, the error estimate, ε* = x̅ⁿ⁺¹  x̃ⁿ⁺¹
has the property that ‖ε*‖ = O(h²)
, and it deviates from the true error, ε
, by ‖ε  ε*‖ = O(h³)
.
To derive the secondorder error estimate, let us first define the vector valued function e(tⁿ, h, xⁿ) = x̅ⁿ⁺¹  xⁿ⁺¹
, the local truncation error for a single, fullsized velocityimplicit Euler integration step, with initial conditions (tⁿ, xⁿ)
, and a step size of h
. Furthermore, use ẍ
to denote df/dt
, and ∇f
and ∇ẍ
to denote the Jacobians df/dx
and dẍ/dx
of the ODE system ẋ = f(t, x)
. Note that ẍ
uses a total time derivative, i.e., ẍ = ∂f/∂t + ∇f f
.
Let us use x*
to denote the true solution after a halfstep, x(tⁿ+½h)
, and x̃*
to denote the velocityimplicit Euler solution after a single halfsized step. Furthermore, let us use xⁿ*¹
to denote the true solution of the system at time t = tⁿ+h
if the system were at x̃*
when t = tⁿ+½h
. See the following diagram for an illustration.
Legend: ───── propagation along the true system :···· propagation using implicit Euler with a half step : propagation using implicit Euler with a full step Time tⁿ tⁿ+½h tⁿ+h State : x̅ⁿ⁺¹ <─── used for error estimation : : : : :·········· x̃ⁿ⁺¹ <─── propagated result : : :········· x̃* ─────── xⁿ*¹ : xⁿ ─────── x* ─────── xⁿ⁺¹ <─── true solution
We will use superscripts to denote evaluating an expression with x
at that subscript and t
at the corresponding time, e.g. ẍⁿ
denotes ẍ(tⁿ, xⁿ)
, and f*
denotes f(tⁿ+½h, x*)
. We first present a shortened derivation, followed by the longer, detailed version.
We know the local truncation error for the implicit Euler method is:
e(tⁿ, h, xⁿ) = x̅ⁿ⁺¹  xⁿ⁺¹ = ½ h²ẍⁿ + O(h³). (10)
The local truncation error ε from taking two half steps is composed of these two terms:
e₁ = xⁿ*¹  xⁿ⁺¹ = (1/8) h²ẍⁿ + O(h³), (15) e₂ = x̃ⁿ⁺¹  xⁿ*¹ = (1/8) h²ẍⁿ + O(h³). (20)
Taking the sum,
ε = x̃ⁿ⁺¹  xⁿ⁺¹ = e₁ + e₂ = (1/4) h²ẍⁿ + O(h³). (21)
These two estimations allow us to obtain an estimation of the local error from the difference between the available quantities x̅ⁿ⁺¹ and x̃ⁿ⁺¹:
ε* = x̅ⁿ⁺¹  x̃ⁿ⁺¹ = e(tⁿ, h, xⁿ)  ε, = (1/4) h²ẍⁿ + O(h³), (22)
and therefore our error estimate is second order.
Below we will show this derivation in detail along with the proof that ‖ε  ε*‖ = O(h³)
:
Let us look at a single velocityimplicit Euler step. Upon NewtonRaphson convergence, the truncation error for velocityimplicit Euler, which is the same as the truncation error for implicit Euler (because both methods solve Eqs. (34)), is
e(tⁿ, h, xⁿ) = ½ h²ẍⁿ⁺¹ + O(h³) = ½ h²ẍⁿ + O(h³). (10)
To see why the two are equivalent, we can Taylor expand about (tⁿ, xⁿ)
,
ẍⁿ⁺¹ = ẍⁿ + h dẍ/dtⁿ + O(h²) = ẍⁿ + O(h). e(tⁿ, h, xⁿ) = ½ h²ẍⁿ⁺¹ + O(h³) = ½ h²(ẍⁿ + O(h)) + O(h³) = ½ h²ẍⁿ + O(h³).
Moving on with our derivation, after one small halfsized implicit Euler step, the solution x̃*
is
x̃* = x* + e(tⁿ, ½h, xⁿ) = x* + (1/8) h²ẍⁿ + O(h³), x̃*  x* = (1/8) h²ẍⁿ + O(h³). (11)
Taylor expanding about t = tⁿ+½h
in this x = x̃*
alternate reality,
xⁿ*¹ = x̃* + ½h f(tⁿ+½h, x̃*) + O(h²). (12)
Similarly, Taylor expansions about t = tⁿ+½h
and the true solution x = x*
also give us
xⁿ⁺¹ = x* + ½h f* + O(h²), (13) f(tⁿ+½h, x̃*) = f* + (∇f*) (x̃*  x*) + O(‖x̃*  x*‖²) = f* + O(h²), (14)
where in the last line we substituted Eq. (11).
Eq. (12) minus Eq. (13) gives us,
xⁿ*¹  xⁿ⁺¹ = x̃*  x* + ½h(f(tⁿ+½h, x̃*)  f*) + O(h³), = x̃*  x* + O(h³),
where we just substituted in Eq. (14). Finally, substituting in Eq. (11),
e₁ = xⁿ*¹  xⁿ⁺¹ = (1/8) h²ẍⁿ + O(h³). (15)
After the second small step, the solution x̃ⁿ⁺¹
is
x̃ⁿ⁺¹ = xⁿ*¹ + e(tⁿ+½h, ½h, x̃*), = xⁿ*¹ + (1/8)h² ẍ(tⁿ+½h, x̃*) + O(h³). (16)
Taking Taylor expansions about (tⁿ, xⁿ)
,
x* = xⁿ + ½h fⁿ + O(h²) = xⁿ + O(h). (17) x̃*  xⁿ = (x̃*  x*) + (x*  xⁿ) = O(h), (18)
where we substituted in Eqs. (11) and (17), and
ẍ(tⁿ+½h, x̃*) = ẍⁿ + ½h ∂ẍ/∂tⁿ + ∇ẍⁿ (x̃*  xⁿ) + O(h ‖x̃*  xⁿ‖) = ẍⁿ + O(h), (19)
where we substituted in Eq. (18).
Substituting Eqs. (19) and (15) into Eq. (16),
x̃ⁿ⁺¹ = xⁿ*¹ + (1/8) h²ẍⁿ + O(h³) (20) = xⁿ⁺¹ + (1/4) h²ẍⁿ + O(h³),
therefore
ε = x̃ⁿ⁺¹  xⁿ⁺¹ = (1/4) h² ẍⁿ + O(h³). (21)
Subtracting Eq. (21) from Eq. (10),
e(tⁿ, h, xⁿ)  ε = (½  1/4) h²ẍⁿ + O(h³); ⇒ ε* = x̅ⁿ⁺¹  x̃ⁿ⁺¹ = (1/4) h²ẍⁿ + O(h³). (22)
Eq. (22) shows that our error estimate is secondorder. Since the first term on the RHS matches ε
(Eq. (21)),
ε* = ε + O(h³). (23)
Implements IntegratorBase< T >.

delete 

delete 

finalvirtual 
Returns true, because this integrator supports error estimation.
Implements IntegratorBase< T >.