BDF (Backward Differentiation Formula) time integration scheme for a system of ODEs for an initial value problem (IVP) More...

#include <Bdf.h>

Public Types
using	ScalarType = TScalar
	Floating point scalar type.

using	IndexType = TIndex
	Integer index type.

Public Member Functions
	Bdf (int step=1, int order=2)
	Construct a `step`-step BDF (backward differentiation formula) time integration scheme for a system of ODEs (ordinary differential equation) of order `order`.

auto	Order () const
	Order of the ODE.

auto	Step () const
	Step `s` of the `s`-step BDF scheme.

auto	Dimensions () const
	Number of ODEs.

auto	TimeStep () const
	Time step size.

auto	Inertia (int o) const
	Inertia of the BDF scheme for the $ o^\text{th} $ state derivative.

auto	State (int k, int o=0) const -> decltype(xt.col(0))
	$ o^\text{th} $ state derivative

auto	State (int k, int o=0) -> decltype(xt.col(0))
	$ o^\text{th} $ state derivative

auto	CurrentState (int o=0) const -> decltype(xt.col(0))
	Current state derivative.

auto	CurrentState (int o=0) -> decltype(xt.col(0))
	Current state derivative.

auto	Alpha () const
	Interpolation coefficients $ \alpha_k $ except $ \alpha_s $.

auto	Beta () const
	Forcing term coefficient $ \beta $ s.t. $ \tilde{\beta} = h \beta $.

auto	BetaTilde () const
	Time-step scaled forcing term coefficient $ \tilde{\beta} = h \beta $.

void	SetOrder (int order)
	Set the order of the ODE system.

void	SetStep (int step)
	Set the step of the BDF scheme.

void	SetTimeStep (ScalarType dt)
	Set the time step size.

void	ConstructEquations ()
	Construct the BDF equations, i.e. compute $ \tilde{x^{(o)}} = \sum_{k=0}^{s-1} \alpha_k x^{(o)}_{t_i - s + k} $ for all $ o = 0, ..., \text{order}-1 $.

template<class TDerivedX>
void	SetInitialConditions (Eigen::DenseBase< TDerivedX > const &x0)
	Set the initial conditions for the initial value problem.

template<class... TDerivedX>
void	SetInitialConditions (Eigen::DenseBase< TDerivedX > const &... x0)
	Set the initial conditions for the initial value problem.

template<class TDerivedX>
void	Step (Eigen::DenseBase< TDerivedX > const &x)
	Advance the BDF scheme by one time step.

template<class... TDerivedX>
void	Step (Eigen::DenseBase< TDerivedX > const &... xs)
	Advance the BDF scheme by one time step.

void	Tick ()
	Advance the BDF scheme by one time step.

void	Serialize (io::Archive &archive) const
	Serialize to HDF5 group.

void	Deserialize (io::Archive const &archive)
	Deserialize from HDF5 group.

Public Attributes
Eigen::Matrix< ScalarType, Eigen::Dynamic, Eigen::Dynamic >	xt

Eigen::Matrix< ScalarType, Eigen::Dynamic, Eigen::Dynamic >	xtilde

IndexType	ti
	Current time index s.t. $t = t_0 + h t_i$.

ScalarType	h
	Time step size $ h $.

Detailed Description

template<class TScalar = Scalar, class TIndex = Index>
class pbat::sim::integration::Bdf< TScalar, TIndex >

BDF (Backward Differentiation Formula) time integration scheme for a system of ODEs for an initial value problem (IVP)

Refer to background here.

An order $ p $ system of ODEs $ x^{(p)} = f(t, \frac{d^{p-1}}{dt^{p-1}} x, \dots, x) $ can be transformed into a system of $ p $ first-order ODEs using slack variables $ x^{(o)} = \frac{d^o}{dt^o} x $ for $ o = 0, \dots, p-1 $, such that the IVP can be rewritten as

\[ \begin{align*} \frac{d}{dt} x_{p-1} &= f(t, x^{(p-1)}, \dots, x^{(o)}) \\ \frac{d}{dt} x^{(p-j)} &= x^{(p-j+1)}, \quad j = 2, \dots, p \end{align*} . \]

Since BDF discretizes each equation as $ \sum_{k=0}^s \alpha_k x_{n-s+k} = h \beta f(t, x_n^{(p-1)}, \dots, x_n) $, we have that

\[ \begin{align*} x^{(p-1)}_n + \tilde{x}^{(p-1)}_\text{BDFS} - \tilde{\beta}_\text{BDFS} f(t, x_n^{(p-1)}, \dots, x_n) &= 0 \\ x^{(o)}_n + \tilde{x}^{(o)}_\text{BDFS} - \tilde{\beta}_\text{BDFS} x^{(o+1)}_n &= 0 \end{align*} \]

for $ o=0,\dots,p-2 $, where $ \tilde{\beta} = h \beta $ and $ n $ is the next time step index.

One performs time integration by solving for the states and their derivatives $ x_n^{(o)} $ for $ o=0,\dots,p-1 $ at each time step using root-finding of the above equations. Alternatively, the root-finding equations can be treated as stationarity conditions $ \nabla f = 0 $ for some objective function $ f $, and solved via numerical optimization of $ f $.

This class encapsulates the storage of the past states and their lower order derivatives $ x_{n-s+k}^{(o)} $, the construction of so-called inertias $ \tilde{x}^{(o)}_\text{BDFS} $, the interpolation coefficients $ \alpha_k $, forcing term coefficient $ \beta $ and generalize the BDF scheme to various ODE systems of different orders. Specific ODEs are then entirely defined by the forcing function $ f(t, \frac{d^{p-1}}{dt^{p-1}} x, \dots, x) $, and it is up to the user to derive their specific equations to solve.

bdf.SetInitialConditions(x0, v0)
for ti in steps:
  bdf.ConstructEquations()
  xn, vn = userSolve(bdf)
  bdf.Step(xn, vn)

Constructor & Destructor Documentation

◆ Bdf()

template<class TScalar, class TIndex>

pbat::sim::integration::Bdf< TScalar, TIndex >::Bdf	(	int	step = 1,
		int	order = 2 )

Construct a step-step BDF (backward differentiation formula) time integration scheme for a system of ODEs (ordinary differential equation) of order order.

Parameters

step	`0 < s < 7` backward differentiation scheme
order	`order > 0` order of the ODE

Precondition: 0 < step < 7; order > 0

Member Function Documentation

◆ Alpha()

template<class TScalar = Scalar, class TIndex = Index>

auto pbat::sim::integration::Bdf< TScalar, TIndex >::Alpha ( ) const

inline

Interpolation coefficients $ \alpha_k $ except $ \alpha_s $.

Returns: s x 1 vector of interpolation coefficients $ \alpha_k $

◆ Beta()

template<class TScalar = Scalar, class TIndex = Index>

auto pbat::sim::integration::Bdf< TScalar, TIndex >::Beta ( ) const

inline

Forcing term coefficient $ \beta $ s.t. $ \tilde{\beta} = h \beta $.

Returns: Forcing term coefficient $ \beta $

◆ BetaTilde()

template<class TScalar = Scalar, class TIndex = Index>

auto pbat::sim::integration::Bdf< TScalar, TIndex >::BetaTilde ( ) const

inline

Time-step scaled forcing term coefficient $ \tilde{\beta} = h \beta $.

Returns: Time-step scaled forcing term coefficient $ \tilde{\beta} = h \beta $

◆ CurrentState() [1/2]

template<class TScalar, class TIndex>

auto pbat::sim::integration::Bdf< TScalar, TIndex >::CurrentState ( int o = 0 ) -> decltype(xt.col(0))

Current state derivative.

Parameters

o	Order of the state derivative $ o = 0, ..., \text{order}-1 $

Returns: n x 1 current state derivative vector $ x^{(o)}_{t_i} $

◆ CurrentState() [2/2]

template<class TScalar, class TIndex>

auto pbat::sim::integration::Bdf< TScalar, TIndex >::CurrentState ( int o = 0 ) const -> decltype(xt.col(0))

Current state derivative.

Parameters

o	Order of the state derivative $ o = 0, ..., \text{order}-1 $

Returns: n x 1 current state derivative vector $ x^{(o)}_{t_i} $

◆ Deserialize()

template<class TScalar, class TIndex>

void pbat::sim::integration::Bdf< TScalar, TIndex >::Deserialize ( io::Archive const & archive )

Deserialize from HDF5 group.

Parameters

archive Archive to deserialize from

◆ Dimensions()

template<class TScalar = Scalar, class TIndex = Index>

auto pbat::sim::integration::Bdf< TScalar, TIndex >::Dimensions ( ) const

inline

Number of ODEs.

Returns: Number of ODEs

◆ Inertia()

template<class TScalar = Scalar, class TIndex = Index>

auto pbat::sim::integration::Bdf< TScalar, TIndex >::Inertia ( int o ) const

inline

Inertia of the BDF scheme for the $ o^\text{th} $ state derivative.

Parameters

o	Order of the state derivative $ o = 0, ..., \text{order}-1 $

Returns: n x 1 inertia vector for the $ o^\text{th} $ state derivative

◆ Order()

template<class TScalar = Scalar, class TIndex = Index>

auto pbat::sim::integration::Bdf< TScalar, TIndex >::Order ( ) const

inline

Order of the ODE.

Returns: Order of the ODE

◆ Serialize()

template<class TScalar, class TIndex>

void pbat::sim::integration::Bdf< TScalar, TIndex >::Serialize ( io::Archive & archive ) const

Serialize to HDF5 group.

Parameters

archive Archive to serialize to

◆ SetInitialConditions() [1/2]

template<class TScalar, class TIndex>

template<class... TDerivedX>

void pbat::sim::integration::Bdf< TScalar, TIndex >::SetInitialConditions ( Eigen::DenseBase< TDerivedX > const &... x0 )

inline

Set the initial conditions for the initial value problem.

Parameters

x0	`order` vectors of `n x 1` initial conditions $ x^{(o)}_{t_0} $

Precondition: sizeof...(TDerivedX) == order

Postcondition: ti == 0

◆ SetInitialConditions() [2/2]

template<class TScalar, class TIndex>

template<class TDerivedX>

void pbat::sim::integration::Bdf< TScalar, TIndex >::SetInitialConditions ( Eigen::DenseBase< TDerivedX > const & x0 )

inline

Set the initial conditions for the initial value problem.

Parameters

x0	`n x order` matrix of initial conditions s.t. `x0.col(o) = \f$ x^{(o)}_{t_0} \f$` for $ o = 0, ..., \text{order}-1 $

Precondition: x0.cols() == order

Postcondition: ti == 0

◆ SetOrder()

template<class TScalar, class TIndex>

void pbat::sim::integration::Bdf< TScalar, TIndex >::SetOrder ( int order )

Set the order of the ODE system.

Parameters

order Order of the ODE system $ \text{order} > 0 $

Precondition: order > 0

◆ SetStep()

template<class TScalar, class TIndex>

void pbat::sim::integration::Bdf< TScalar, TIndex >::SetStep ( int step )

Set the step of the BDF scheme.

Parameters

step	Step of the BDF scheme $ 0 < s < 7 $

Precondition: 0 < step < 7

◆ SetTimeStep()

template<class TScalar, class TIndex>

void pbat::sim::integration::Bdf< TScalar, TIndex >::SetTimeStep ( ScalarType dt )

Set the time step size.

Parameters

dt	Time step size $ dt > 0 $

Precondition: dt > 0

◆ State() [1/2]

template<class TScalar, class TIndex>

auto pbat::sim::integration::Bdf< TScalar, TIndex >::State	(	int	k,
		int	o = 0 ) -> decltype(xt.col(0))

$ o^\text{th} $ state derivative

Parameters

k	State index $ k = 0, ..., s $ for the vector $ x^{(o)}_{t_i - s + k} $
o	Order of the state derivative $ o = 0, ..., \text{order}-1 $

Returns: n x 1 state derivative vector $ x^{(o)}_{t_i - s + k} $

◆ State() [2/2]

template<class TScalar, class TIndex>

auto pbat::sim::integration::Bdf< TScalar, TIndex >::State	(	int	k,
		int	o = 0 ) const -> decltype(xt.col(0))

$ o^\text{th} $ state derivative

Parameters

k	State index $ k = 0, ..., s $ for the vector $ x^{(o)}_{t_i - s + k} $
o	Order of the state derivative $ o = 0, ..., \text{order}-1 $

Returns: n x 1 state derivative vector $ x^{(o)}_{t_i - s + k} $

◆ Step() [1/3]

template<class TScalar = Scalar, class TIndex = Index>

auto pbat::sim::integration::Bdf< TScalar, TIndex >::Step ( ) const

inline

Step s of the s-step BDF scheme.

Returns: Step s of the s-step BDF scheme

◆ Step() [2/3]

template<class TScalar, class TIndex>

template<class... TDerivedX>

void pbat::sim::integration::Bdf< TScalar, TIndex >::Step ( Eigen::DenseBase< TDerivedX > const &... xs )

inline

Advance the BDF scheme by one time step.

Template Parameters

TDerivedX Derived type of the input matrix

Parameters

xs	`order` vectors of `n x 1` current state derivatives $ x_{t_i}^{(o)} $

Precondition: sizeof...(TDerivedX) == order

◆ Step() [3/3]

template<class TScalar, class TIndex>

template<class TDerivedX>

void pbat::sim::integration::Bdf< TScalar, TIndex >::Step ( Eigen::DenseBase< TDerivedX > const & x )

inline

Advance the BDF scheme by one time step.

Template Parameters

TDerivedX Derived type of the input matrix

Parameters

x	`n x order` matrix of the current state derivatives $ x_{t_i}^{(o)} $

◆ TimeStep()

template<class TScalar = Scalar, class TIndex = Index>

auto pbat::sim::integration::Bdf< TScalar, TIndex >::TimeStep ( ) const

inline

Time step size.

Returns: Time step size

Member Data Documentation

◆ xt

template<class TScalar = Scalar, class TIndex = Index>

Eigen::Matrix<ScalarType, Eigen::Dynamic, Eigen::Dynamic> pbat::sim::integration::Bdf< TScalar, TIndex >::xt

n x |s*order| matrix of n-dimensional states and their time derivatives s.t. $ xt.col(o*s + k) = x^(k)(t - k*dt) $ for $ k = 0, ..., s $ and $ o = 0, ..., \text{order}-1 $

◆ xtilde

template<class TScalar = Scalar, class TIndex = Index>

Eigen::Matrix<ScalarType, Eigen::Dynamic, Eigen::Dynamic> pbat::sim::integration::Bdf< TScalar, TIndex >::xtilde

n x order matrix of n-dimensional aggregated past states and time derivatives s.t. xtilde.col(o) = $ \frac{1}{\alpha_s} \sum_{k=t_i-s}^{s-1} < \alpha_k x_k $ for $ o = 0, ..., \text{order}-1 $

The documentation for this class was generated from the following file:

C:/git/PhysicsBasedAnimationToolkit/source/pbat/sim/integration/Bdf.h

k	State index \( k = 0, ..., s \) for the vector \( x^{(o)}_{t_i - s + k} \)
o	Order of the state derivative \( o = 0, ..., \text{order}-1 \)

Public Types

Public Member Functions

Public Attributes

Detailed Description

Constructor & Destructor Documentation

◆ Bdf()

Member Function Documentation

◆ Alpha()

◆ Beta()

◆ BetaTilde()

◆ CurrentState() [1/2]

◆ CurrentState() [2/2]

◆ Deserialize()

◆ Dimensions()

◆ Inertia()

◆ Order()

◆ Serialize()

◆ SetInitialConditions() [1/2]

◆ SetInitialConditions() [2/2]

◆ SetOrder()

◆ SetStep()

◆ SetTimeStep()

◆ State() [1/2]

◆ State() [2/2]

◆ Step() [1/3]

◆ Step() [2/3]

◆ Step() [3/3]

◆ TimeStep()

Member Data Documentation

◆ xt

◆ xtilde