SurgSim::Math::OdeSolverEulerImplicit Class Reference

Euler implicit (a.k.a backward Euler) ode solver (see OdeSolverEulerImplicit.dox for more details). More...

`#include <SurgSim/Math/OdeSolverEulerImplicit.h>`

Inheritance diagram for SurgSim::Math::OdeSolverEulerImplicit:

## Public Member Functions

| OdeSolverEulerImplicit (OdeEquation *equation)
| Constructor. More...
virtual void | setNewtonRaphsonMaximumIteration (size_t maximumIteration)
size_t | getNewtonRaphsonMaximumIteration () const
void | setNewtonRaphsonEpsilonConvergence (double epsilonConvergence)
double | getNewtonRaphsonEpsilonConvergence () const
void | solve (double dt, const OdeState &currentState, OdeState *newState, bool computeCompliance=true) override
| Solves the equation. More...
Public Member Functions inherited from SurgSim::Math::OdeSolver
| OdeSolver (OdeEquation *equation)
| Constructor. More...
virtual | ~OdeSolver ()
| Virtual destructor. More...
const std::string | getName () const
| Gets the solver's name. More...
void | setLinearSolver (std::shared_ptr< LinearSparseSolveAndInverse > linearSolver)
| Sets the specialized linear solver to use with this Ode solver. More...
std::shared_ptr< LinearSparseSolveAndInverse > | getLinearSolver () const
| Gets the specialized linear solver used with this Ode solver. More...
void | computeMatrices (double dt, const OdeState &state, bool computeCompliance=true)
| Computes the system and compliance matrices for a given state. More...
const SparseMatrix & | getSystemMatrix () const
| Queries the current system matrix. More...
const Matrix & | getComplianceMatrix () const

## Protected Member Functions

void | assembleLinearSystem (double dt, const OdeState &state, const OdeState &newState, bool computeRHS=true) override
| Assemble the linear system (A.x=b) to be solved for the state and new states (useful for certain ode solver). More...
Protected Member Functions inherited from SurgSim::Math::OdeSolver
void | computeComplianceMatrixFromSystemMatrix (const OdeState &state)
| Helper method computing the compliance matrix from the system matrix and setting the boundary conditions. More...

## Protected Attributes

size_t | m_maximumIteration
| Newton-Raphson maximum number of iteration (1 => linearization) More...
double | m_epsilonConvergence
| Newton-Raphson convergence criteria (variation of the solution over time) More...
Vector | m_previousSolution
| Newton-Raphson previous solution (we solve a problem to find deltaV, the variation in velocity) More...
Protected Attributes inherited from SurgSim::Math::OdeSolver
std::string | m_name
| Name for this solver. More...
OdeEquation & | m_equation
| The ode equation (API providing the necessary evaluation methods and the initial state) More...
std::shared_ptr< LinearSparseSolveAndInverse > | m_linearSolver
| The specialized linear solver to use when solving the ode equation. More...
SparseMatrix | m_systemMatrix
| Linear system matrix (can be M, K, combination of MDK depending on the solver), including boundary conditions. More...
Vector | m_solution
| Linear system solution and rhs vectors (including boundary conditions) More...
Vector | m_rhs
Matrix | m_complianceMatrix
| Compliance matrix which is the inverse of the system matrix, including boundary conditions. More...

## Detailed Description

Euler implicit (a.k.a backward Euler) ode solver (see OdeSolverEulerImplicit.dox for more details).

The ODE to be solved is of the form:

\[ \displaystyle M.a(t) = f(t, x(t), v(t)) \]

.

Note
  \(f\) is supposed to be continuous with partial derivatives \(\displaystyle\frac{\partial f}{\partial x} = -K\) and \(\displaystyle\frac{\partial f}{\partial v} = -D\).
  Also, we define the system matrix \(\displaystyle S = \frac{M}{dt} + D + dt.K\).

ODE of order 1
This ODE is solved as an ODE of order 1 by defining the state vector \(y = \left(\begin{array}{c} x \\ v \end{array}\right)\):

\[ \displaystyle y' = \left(\begin{array}{c} x' \\ v' \end{array}\right) = \left(\begin{array}{c} v \\ M^{-1}.f(t, x, v) \end{array}\right) = f(t, y) \]

Integrating this equation gives:

\[ \displaystyle y(t+dt) - y(t) = \int_{t}^{t+dt} f(t, y) dt \]

Backward Euler to numerically integrate the ODE
Euler Implicit uses a rectangular numerical integration evaluated on the right: \(\int_{t}^{t+dt} f(t, y) dt \simeq dt. f(t+dt, y(t+dt))\). Unlike Explicit Euler (a.k.a. forward Euler) which uses a rectangular numerical integration evaluated on the left: \(\int_{t}^{t+dt} f(t, y) dt \simeq dt. f(t, y(t))\).
Euler implicit is also called backward Euler because the tangent evaluation is numerically evaluated using a backward evaluation: \(y'(t) = (y(t) - y(t-dt)) / dt\) which leads to the well known integration scheme (written using the state vector \(y\) or the variables \(x\) and \(v\):

\[ \begin{array}{ccc} y(t+dt) = y(t) + dt. f(t+dt, y(t+dt)) & or & \left\{ \begin{array}{ccccl} x(t+dt) &=& x(t) &+& dt.v(t+dt) \\ v(t+dt) &=& v(t) &+& dt.M^{-1}.f(t+dt, x(t+dt), v(t+dt)) \end{array} \right. \end{array} \]

This system of equations is non-linear w.r.t. to the unknown \(y(t+dt)\).

Newton-Raphson to solve the non-linear resulting equations
We use Newton-Raphson to solve this non-liner problem http://en.wikipedia.org/wiki/Newton%27s_method. Let's pose the non-linear problem to solve as:

\[ F(y) = y - y(t) - dt.f(t+dt, y) = 0 \]

The solution to this non-linear equation is \(y(t+dt)\).

The Newton-Raphson algorithm iteratively calculates a series of state \(y_n\) converging to the solution \(y(t+dt)\) (if the problem has a solution and the initial solution is close enough).
Each iteration computes:

\[ y_{n+1} = y_{n} - \frac{\partial F}{\partial y}(y_{n})^{-1} F(y_{n}) \]

where \(y_n = \left(\begin{array}{c}x_n\\v_n\end{array}\right)\) and \(y_{n+1} = \left(\begin{array}{c}x_{n+1}\\v_{n+1}\end{array}\right)\) are 2 successive approximations of the solution. They are neither \(y(t) = \left(\begin{array}{c}x(t)\\v(t)\end{array}\right)\) nor the solution \(y(t+dt) = \left(\begin{array}{c}x(t+dt)\\v(t+dt)\end{array}\right)\). Except on the very first iteration where \(y_0 = y(t)\).

Note
  The initial solution \(y_0\) is set to be \(y(t)\) as it is the latest solution and should by nature be close to the solution \(y(t+dt)\).
  Therefore, on the first iteration, \(F(y_0)\) will evaluate to exactly \(-dt.f(t+dt, y(t))\) which is the force evaluation from the previous time step.

\[ \frac{\partial F}{\partial y} = \left( I - dt.\frac{\partial f}{\partial y}\right) = \left( \begin{array}{cc} I & -dt.I \\ dt.M^{-1}K & I + dt.M^{-1}D \end{array} \right) \]

It can be easily shown using a Gauss pivot technique that the inverse is:

\[ \frac{\partial F}{\partial y}^{-1} = \left( \begin{array}{cc} I - dt.S^{-1}K & S^{-1}M \\ -S^{-1}K & \frac{S^{-1}M}{dt} \end{array} \right) \]

Therefore, we can formally develop the solution:

\[ y_{n+1} = \left\{ \begin{array}{ccccl} x_{n+1} &=& x_n &-& \left[\left(I - dt.S^{-1}K\right).\left(x_n - x(t) - dt.v_n\right) + S^{-1}M.\left(v_n - v(t) - dt.M^{-1}.f(x_n, v_n)\right)\right] \\ &=& x_n &-& x_n + x(t) + dt.v_n + dt.S^{-1}K\left(x_n - x(t) - dt.v_n\right) - S^{-1}M.\left(v_n - v(t) - dt.M^{-1}.f(x_n, v_n)\right) \\ &=& x(t) &+& dt.\left[v_n + S^{-1}\left(f(x_n, v_n) + K(x_n - x(t) - dt.v_n) - \frac{M}{dt}(v_n - v(t))\right)\right] \\ &=& x(t) &+& dt.v_{n+1} \\ v_{n+1} &=& v_n &-& \left[-S^{-1}K.\left(x_n - x(t) - dt.v_n\right) + \frac{S^{-1}M}{dt}.\left(v_n - v(t) - dt.M^{-1}.f(x_n, v_n)\right)\right] \\ &=& v_n &+& S^{-1}\left(f(x_n, v_n) + K(x_n - x(t) - dt.v_n) - \frac{M}{dt}(v_n - v(t)) \right) \end{array} \right. \]

We simply need to solve the system to find the velocity variation and we can deduct the new position from there.

## Constructor & Destructor Documentation

## OdeSolverEulerImplicit()

SurgSim::Math::OdeSolverEulerImplicit::OdeSolverEulerImplicit | ( | OdeEquation * | equation | ) |

explicit

Constructor.

Parameters

  equation | The ode equation to be solved

## Member Function Documentation

## assembleLinearSystem()

void SurgSim::Math::OdeSolverEulerImplicit::assembleLinearSystem | ( | double | dt,
|  | const OdeState & | state,
|  | const OdeState & | newState,
|  | bool | computeRHS = `true`
| ) |  |

override protected virtual

Assemble the linear system (A.x=b) to be solved for the state and new states (useful for certain ode solver).

Parameters

  dt | The time step used in the system
  state,newState | The state and newState to be used to evaluate the system
  computeRHS | True to compute the RHS vector, False otherwise

Note
  The method should fill up the LHS matrix in m_systemMatrix and the RHS vector in m_b (if requested)
  The method should take care of the boundary conditions properly on both the matrix and the vector.
  The method should prepare the linear solver m_linearSolver to be used with the m_systemMatrix

Implements SurgSim::Math::OdeSolver.

## getNewtonRaphsonEpsilonConvergence()

double SurgSim::Math::OdeSolverEulerImplicit::getNewtonRaphsonEpsilonConvergence | ( |  | ) | const

Returns
  The Newton-Raphson algorithm epsilon convergence

## getNewtonRaphsonMaximumIteration()

size_t SurgSim::Math::OdeSolverEulerImplicit::getNewtonRaphsonMaximumIteration | ( |  | ) | const

Returns
  The Newton-Raphson algorithm maximum number of iterations

## setNewtonRaphsonEpsilonConvergence()

void SurgSim::Math::OdeSolverEulerImplicit::setNewtonRaphsonEpsilonConvergence | ( | double | epsilonConvergence | ) |

Parameters

  epsilonConvergence | The Newton-Raphson algorithm epsilon convergence

## setNewtonRaphsonMaximumIteration()

void SurgSim::Math::OdeSolverEulerImplicit::setNewtonRaphsonMaximumIteration | ( | size_t | maximumIteration | ) |

virtual

Parameters

  maximumIteration | The Newton-Raphson algorithm maximum number of iterations

Reimplemented in SurgSim::Math::OdeSolverLinearEulerImplicit.

## solve()

void SurgSim::Math::OdeSolverEulerImplicit::solve | ( | double | dt,
|  | const OdeState & | currentState,
|  | OdeState * | newState,
|  | bool | computeCompliance = `true`
| ) |  |

override virtual

Solves the equation.

Parameters

  | dt | The time step
  | currentState | State at time t
  [out] | newState | State at time t+dt
  | computeCompliance | True to explicitly compute the compliance matrix, False otherwise

Implements SurgSim::Math::OdeSolver.

Reimplemented in SurgSim::Math::OdeSolverLinearEulerImplicit.

## Member Data Documentation

## m_epsilonConvergence

double SurgSim::Math::OdeSolverEulerImplicit::m_epsilonConvergence

protected

Newton-Raphson convergence criteria (variation of the solution over time)

## m_maximumIteration

size_t SurgSim::Math::OdeSolverEulerImplicit::m_maximumIteration

protected

Newton-Raphson maximum number of iteration (1 => linearization)

## m_previousSolution

Vector SurgSim::Math::OdeSolverEulerImplicit::m_previousSolution

protected

Newton-Raphson previous solution (we solve a problem to find deltaV, the variation in velocity)

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

- SurgSim/Math/OdeSolverEulerImplicit.h
- SurgSim/Math/OdeSolverEulerImplicit.cpp