Euler implicit (a.k.a backward Euler) ode solver (see OdeSolverEulerImplicit.dox for more details). More...
#include <OdeSolverEulerImplicit.h>
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 ¤tState, 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. | |
| 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) | |
| double | m_epsilonConvergence |
| Newton-Raphson convergence criteria (variation of the solution over time) | |
| Vector | m_previousSolution |
| Newton-Raphson previous solution (we solve a problem to find deltaV, the variation in velocity) | |
| 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) | |
| std::shared_ptr< LinearSparseSolveAndInverse > | m_linearSolver |
| The specialized linear solver to use when solving the ode equation. | |
| 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) | |
| Vector | m_rhs |
| Matrix | m_complianceMatrix |
| Compliance matrix which is the inverse of the system matrix, including boundary conditions. | |
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()
|
explicit |
Constructor.
- Parameters
-
equation The ode equation to be solved
Member Function Documentation
§ assembleLinearSystem()
|
overrideprotectedvirtual |
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()
|
virtual |
- Parameters
-
maximumIteration The Newton-Raphson algorithm maximum number of iterations
Reimplemented in SurgSim::Math::OdeSolverLinearEulerImplicit.
§ solve()
|
overridevirtual |
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.
The documentation for this class was generated from the following files:
- SurgSim/Math/OdeSolverEulerImplicit.h
- SurgSim/Math/OdeSolverEulerImplicit.cpp
