Back

EulerImplicitSolver

The Euler Implicit Solver is a time integrator that uses the implicit backward Euler scheme to solve dynamic systems. This method is unconditionally stable, making it suitable for stiff systems where explicit methods might fail due to stability constraints.

abstract
The EulerImplicitSolver implements an implicit backward Euler scheme for solving first- and second-order ODEs in SOFA simulations, supporting Rayleigh damping and trapezoidal rule options.
sheet
# EulerImplicitSolver ## Overview EulerImplicitSolver is a time integrator that uses the implicit backward Euler scheme to solve ordinary differential equations (ODEs) of first or second order. It supports Rayleigh damping and can optionally use the trapezoidal rule for higher accuracy. ## Mathematical Model For a system described by an ODE, the implicit backward Euler method approximates the state at each time step as follows: - **Second Order** \[ x_{t+h} = x_t + h v_{t+h}, \] \[ v_{t+h} = v_t + h a_{t+h}. \] The unknown is \( dv = v_{t+h} - v_t \). Newton's law yields: \[ M dv = h f(t+h), \] where \( M \) is the mass matrix, and \( K = df/dx \) and \( B = df/dv \) are stiffness and damping matrices, respectively. Rayleigh damping adds terms \( -r_M M + r_K K \). The final equation system is: \[ P ( (1+h r_M) M - h B - h(h + r_K) K ) P dv = P h ( f(t) + (h+r_K) K v + B v - r_M M v ). \] - **First Order** For a first-order system, the integration scheme is: \[ x_{t+h} = x_t + h v_{t+h}, \] where \( M v_{t+h} = f_{ext}(t) - h K v_{t+h} \), leading to: \[ (M + h K) v_{t+h} = f_{ext}. \] - **Trapezoidal Rule** For the trapezoidal scheme, the equation for second-order systems becomes: \[ P ( (1+h/2 r_M) M - h/2 B - h/2 (h + r_K) K ) P dv = P h/2 ( 2 f(t) + (h + r_K) K v + B v - r_M M v ). \] ## Parameters and Data - **rayleighStiffness**: Rayleigh damping coefficient related to stiffness, > 0. - **rayleighMass**: Rayleigh damping coefficient related to mass, > 0. - **vdamping**: Velocity decay coefficient (no decay if null). - **firstOrder**: Use backward Euler scheme for first-order ODE system (default: false). - **trapezoidalScheme**: Boolean to use the trapezoidal scheme instead of implicit Euler (default: false). - **solveConstraint**: Apply ConstraintSolver (requires a ConstraintSolver in the same node, default: false). - **threadSafeVisitor**: If true, do not use realloc and free visitors in fwdInteractionForceField. - **computeResidual**: If true, compute residual at the end of solving. ## Practical Notes The EulerImplicitSolver is unconditionally stable but may require careful tuning of parameters such as Rayleigh damping coefficients. The trapezoidal scheme can improve accuracy but increases computational cost.
title
Euler Implicit Solver
description
The Euler Implicit Solver is a time integrator that uses the implicit backward Euler scheme to solve dynamic systems. This method is unconditionally stable, making it suitable for stiff systems where explicit methods might fail due to stability constraints.
parameters
  • {'name': 'linearSolver', 'type': 'core::solver::BaseLinearSolver*', 'description': 'The linear solver used to solve the system of equations arising from the implicit discretization.'}
  • {'name': 'f_rayleighStiffness', 'type': 'Real', 'defaultValue': '0.0', 'description': 'Rayleigh stiffness coefficient, which is part of the Rayleigh damping model that adds both mass-proportional and stiffness-proportional damping to the system.', 'deprecated': True}
  • {'name': 'f_rayleighMass', 'type': 'Real', 'defaultValue': '0.0', 'description': 'Rayleigh mass coefficient, another component of the Rayleigh damping model.', 'deprecated': True}
  • {'name': 'rayleighStiffness', 'type': 'Real', 'defaultValue': '0.0', 'description': 'Deprecation placeholder for the Rayleigh stiffness coefficient.'}
  • {'name': 'rayleighMass', 'type': 'Real', 'defaultValue': '0.0', 'description': 'Deprecation placeholder for the Rayleigh mass coefficient.'}
  • {'name': 'f_printLog', 'type': 'bool', 'defaultValue': 'false', 'description': 'Boolean flag to enable or disable logging of position, velocity, and force at each time step.'}
  • {'name': 'velocityDamping', 'type': 'Real', 'defaultValue': '0.0', 'description': 'Coefficient for exponential damping applied to the velocities after integration.', 'deprecated': True}
  • {'name': 'f_velocityDamping', 'type': 'Real', 'defaultValue': '0.0', 'description': 'Deprecation placeholder for the velocity damping coefficient.'}
  • {'name': 'rayleigh', 'type': 'core::behavior::Rayleigh*', 'description': 'Deprecated pointer to a Rayleigh object used in previous versions for defining Rayleigh damping.', 'deprecated': True}
methods
  • {'name': 'getPositionIntegrationFactor()', 'returnType': 'SReal', 'parameters': [], 'description': 'Returns the integration factor for position based on the current time step.'}
  • {'name': 'getIntegrationFactor(inputDerivative, outputDerivative)', 'returnType': 'SReal', 'parameters': [{'type': 'int', 'name': 'inputDerivative'}, {'type': 'int', 'name': 'outputDerivative'}], 'description': 'Returns the integration factor for converting derivatives from one order to another based on the current time step.'}
  • {'name': 'getSolutionIntegrationFactor(outputDerivative)', 'returnType': 'SReal', 'parameters': [{'type': 'int', 'name': 'outputDerivative'}], 'description': 'Returns the integration factor applied to the solution vector for a given derivative order based on the current time step.'}
maths
The Euler Implicit Solver is a numerical method used for solving ordinary differential equations (ODEs) that arise in the simulation of dynamic systems. It employs the implicit backward Euler scheme, which provides unconditional stability and hence is particularly useful for stiff problems where explicit methods can suffer from severe stability constraints. ### Mathematical Description: Consider a system described by an ODE of the form: \[ \frac{d extbf{x}}{dt} = extbf{f}(t, extbf{x}), \] where \( extbf{x}\) is the state vector representing positions and velocities, and \( extbf{f}\) is a function that describes how the system evolves over time. The Euler Implicit Solver approximates this ODE using the implicit backward Euler method. At each time step \(t_n = n imes dt\), where \(dt\) is the time step size, the state at the next time step \( extbf{x}_{n+1}\) is determined by solving: \[ \frac{ extbf{x}_{n+1} - extbf{x}_n}{dt} = extbf{f}(t_{n+1}, extbf{x}_{n+1}). \] This can be rewritten as a fixed point problem: \[ extbf{x}_{n+1} = extbf{x}_n + dt imes extbf{f}(t_{n+1}, extbf{x}_{n+1}). \] In the context of mechanics, this ODE often takes the form of Newton's second law: \[ M rac{d^2 extbf{x}}{dt^2} + B rac{d extbf{x}}{dt} + K extbf{x} = extbf{F}(t), \] where \(M\) is the mass matrix, \(B\) and \(K\) are damping and stiffness matrices respectively, and \( extbf{F}(t)\) represents external forces. Using the implicit backward Euler method, this becomes: \[ M rac{ extbf{x}_{n+1} - 2 extbf{x}_n + extbf{x}_{n-1}}{dt^2} + B rac{ extbf{x}_{n+1} - extbf{x}_n}{dt} + K extbf{x}_{n+1} = extbf{F}(t_{n+1}). \] This leads to the following linear system at each time step: \[ (M/dt^2 + B/dt + K) extbf{x}_{n+1} = M rac{ extbf{x}_n - 2 extbf{x}_{n-1}}{dt^2} - B rac{ extbf{x}_n}{dt} + extbf{F}(t_{n+1}). \] ### Physical Description: In physical terms, the implicit backward Euler method ensures that both velocity and position are updated in a way that accounts for forces at the new time step. This makes it suitable for simulating systems with strong damping or stiffness, as well as systems where external forces change rapidly. #### Integration Factors: - **Position**: The position update is given by \( extbf{x}_{n+1} = extbf{x}_n + dt imes extbf{v}_{n+1}\), where \( extbf{v}_{n+1}\) is the velocity at time step \(t_{n+1}\). - **Velocity**: The velocity update is given by solving the system \((M/dt^2 + B/dt + K) extbf{x}_{n+1} = extbf{RHS}(dt, t_n, extbf{x}_n)\), where RHS denotes the right-hand side of the equation. The solver ensures that all forces, including constraints and external forces like gravity, are accounted for at each step. Additionally, it supports velocity damping through an exponential decay factor to simulate energy dissipation over time. 
{
  "name": "EulerImplicitSolver",
  "main": {
    "name": "EulerImplicitSolver",
    "namespace": "sofa::component::odesolver::backward",
    "module": "Sofa.Component.ODESolver.Backward",
    "include": "sofa/component/odesolver/backward/EulerImplicitSolver.h",
    "doc": "Time integrator using implicit backward Euler scheme.\n\nSemi-implicit time integrator using backward Euler scheme for first and second degree ODEs. (default: second)\n** 2nd Order ***\nThis is based on [Baraff and Witkin, Large Steps in Cloth Simulation, SIGGRAPH 1998]\nThe integration scheme is based on the following equations:\n  \\f$x_{t+h} = x_t + h v_{t+h}\\f$\n  \\f$v_{t+h} = v_t + h a_{t+h}\\f$\n  The unknown is\n  \\f$v_{t+h} - v_t = dv\\f$\n  Newton's law is\n  \\f$ M dv = h f(t+h) \\f$\n  \\f$ M dv = h ( f(t) + K dx     + (B - r_M M + r_K K) (v+dv) )\\f$\n  \\f$ M dv = h ( f(t) + K h (v+dv) + (B - r_M M + r_K K) (v+dv) )\\f$\n  \\f$ M \\f$ is the mass matrix.\n  \\f$ K = df/dx \\f$ is the stiffness implemented (or not) by the force fields.\n  \\f$ B = df/dv \\f$ is the damping implemented (or not) by the force fields.\n  An additional, uniform Rayleigh damping  \\f$- r_M M + r_K K\\f$ is imposed by the solver.\nThis corresponds to the following equation system:\n  \\f$ ( (1+h r_M) M - h B - h(h + r_K) K ) dv = h ( f(t) + (h+r_K) K v + B v - r_M M v )\\f$\nMoreover, the projective constraints filter out the forbidden motions.\nThis is equivalent with multiplying vectors with a projection matrix \\f$P\\f$.\nFinally, the equation system set by this ode solver is:\n  \\f$ P ( (1+h r_M) M - h B - h(h + r_K) K ) P dv = P h ( f(t) + (h + r_K) K v + B v - r_M M v )\\f$\n** 1st Order ***\nThis integration scheme is based on the following equation:\n  \\f$x_{t+h} = x_t + h v_{t+h}\\f$\nApplied to this mechanical system:\n  \\f$ M v_t = f_{ext} \\f$\n  \\f$ M v_{t+h} = f_{ext_{t+h}} \\f$\n  \\f$           = f_{ext_{t}} + h (df_{ext}/dt)_{t+h} \\f$\n  \\f$           = f_{ext_{t}} + h (df_{ext}/dx)_{t+h} v_{t+h} \\f$\n  \\f$           = f_{ext_{t}} - h K v_{t+h} \\f$\n  \\f$ ( M + h K ) v_{t+h} = f_{ext} \\f$\n** Trapezoidal Rule ***\nThe trapezoidal scheme is based on\n  \\f$v_{t+h} = h/2 ( f(t+h) + f(t) )\\f$\nWith this and the same techniques as for the implicit Euler scheme we receive for *** 2nd Order *** equations\n  \\f$ P ( (1+h/2 r_M) M - h/2 B - h/2 (h + r_K) K ) P dv = P h/2 ( 2 f(t) + (h + r_K) K v + B v - r_M M v )\\f$\nand for *** 1st Order ***\n  \\f$ ( M + h/2 K ) v_{t+h} = f_{ext} \\f$",
    "inherits": [
      "LinearSolverAccessor",
      "OdeSolver"
    ],
    "templates": [],
    "data_fields": [
      {
        "name": "d_rayleighStiffness",
        "type": "SReal",
        "xmlname": "rayleighStiffness",
        "help": "Rayleigh damping coefficient related to stiffness, > 0"
      },
      {
        "name": "d_rayleighMass",
        "type": "SReal",
        "xmlname": "rayleighMass",
        "help": "Rayleigh damping coefficient related to mass, > 0"
      },
      {
        "name": "d_velocityDamping",
        "type": "SReal",
        "xmlname": "vdamping",
        "help": "Velocity decay coefficient (no decay if null)"
      },
      {
        "name": "d_firstOrder",
        "type": "bool",
        "xmlname": "firstOrder",
        "help": "Use backward Euler scheme for first order ODE system, which means that only the first derivative of the DOFs (state) appears in the equation. Higher derivatives are absent"
      },
      {
        "name": "d_trapezoidalScheme",
        "type": "bool",
        "xmlname": "trapezoidalScheme",
        "help": "Boolean to use the trapezoidal scheme instead of the implicit Euler scheme and get second order accuracy in time (false by default)"
      },
      {
        "name": "d_solveConstraint",
        "type": "bool",
        "xmlname": "solveConstraint",
        "help": "Apply ConstraintSolver (requires a ConstraintSolver in the same node as this solver, disabled by by default for now)"
      },
      {
        "name": "d_threadSafeVisitor",
        "type": "bool",
        "xmlname": "threadSafeVisitor",
        "help": "If true, do not use realloc and free visitors in fwdInteractionForceField."
      },
      {
        "name": "d_computeResidual",
        "type": "bool",
        "xmlname": "computeResidual",
        "help": "If true, the residual is computed at the end of the solving"
      },
      {
        "name": "d_residual",
        "type": "SReal",
        "xmlname": "residual",
        "help": "Residual norm at the end of the free-motion solving"
      }
    ],
    "links": [],
    "methods": [
      {
        "name": "init",
        "return_type": "void",
        "params": [],
        "is_virtual": true,
        "is_pure_virtual": false,
        "is_static": false,
        "access": "public"
      },
      {
        "name": "cleanup",
        "return_type": "void",
        "params": [],
        "is_virtual": true,
        "is_pure_virtual": false,
        "is_static": false,
        "access": "public"
      },
      {
        "name": "solve",
        "return_type": "void",
        "params": [
          {
            "name": "params",
            "type": "const core::ExecParams *"
          },
          {
            "name": "dt",
            "type": "SReal"
          },
          {
            "name": "xResult",
            "type": "sofa::core::MultiVecCoordId"
          },
          {
            "name": "vResult",
            "type": "sofa::core::MultiVecDerivId"
          }
        ],
        "is_virtual": true,
        "is_pure_virtual": false,
        "is_static": false,
        "access": "public"
      },
      {
        "name": "getVelocityIntegrationFactor",
        "return_type": "SReal",
        "params": [],
        "is_virtual": true,
        "is_pure_virtual": false,
        "is_static": false,
        "access": "public"
      },
      {
        "name": "getPositionIntegrationFactor",
        "return_type": "SReal",
        "params": [],
        "is_virtual": true,
        "is_pure_virtual": false,
        "is_static": false,
        "access": "public"
      },
      {
        "name": "getPositionIntegrationFactor",
        "return_type": "SReal",
        "params": [
          {
            "name": "dt",
            "type": "SReal"
          }
        ],
        "is_virtual": true,
        "is_pure_virtual": false,
        "is_static": false,
        "access": "public"
      },
      {
        "name": "getIntegrationFactor",
        "return_type": "SReal",
        "params": [
          {
            "name": "inputDerivative",
            "type": "int"
          },
          {
            "name": "outputDerivative",
            "type": "int"
          }
        ],
        "is_virtual": true,
        "is_pure_virtual": false,
        "is_static": false,
        "access": "public"
      },
      {
        "name": "getIntegrationFactor",
        "return_type": "SReal",
        "params": [
          {
            "name": "inputDerivative",
            "type": "int"
          },
          {
            "name": "outputDerivative",
            "type": "int"
          },
          {
            "name": "dt",
            "type": "SReal"
          }
        ],
        "is_virtual": false,
        "is_pure_virtual": false,
        "is_static": false,
        "access": "public"
      },
      {
        "name": "getSolutionIntegrationFactor",
        "return_type": "SReal",
        "params": [
          {
            "name": "outputDerivative",
            "type": "int"
          }
        ],
        "is_virtual": true,
        "is_pure_virtual": false,
        "is_static": false,
        "access": "public"
      },
      {
        "name": "getSolutionIntegrationFactor",
        "return_type": "SReal",
        "params": [
          {
            "name": "outputDerivative",
            "type": "int"
          },
          {
            "name": "dt",
            "type": "SReal"
          }
        ],
        "is_virtual": false,
        "is_pure_virtual": false,
        "is_static": false,
        "access": "public"
      },
      {
        "name": "reallocSolutionVector",
        "return_type": "void",
        "params": [
          {
            "name": "vop",
            "type": "sofa::simulation::common::VectorOperations *"
          }
        ],
        "is_virtual": false,
        "is_pure_virtual": false,
        "is_static": false,
        "access": "protected"
      },
      {
        "name": "reallocRightHandSideVector",
        "return_type": "void",
        "params": [
          {
            "name": "vop",
            "type": "sofa::simulation::common::VectorOperations *"
          }
        ],
        "is_virtual": false,
        "is_pure_virtual": false,
        "is_static": false,
        "access": "protected"
      },
      {
        "name": "reallocResidualVector",
        "return_type": "void",
        "params": [
          {
            "name": "vop",
            "type": "sofa::simulation::common::VectorOperations *"
          }
        ],
        "is_virtual": false,
        "is_pure_virtual": false,
        "is_static": false,
        "access": "protected"
      }
    ]
  },
  "desc": {
    "title": "Euler Implicit Solver",
    "description": "The Euler Implicit Solver is a time integrator that uses the implicit backward Euler scheme to solve dynamic systems. This method is unconditionally stable, making it suitable for stiff systems where explicit methods might fail due to stability constraints.",
    "parameters": [
      {
        "name": "linearSolver",
        "type": "core::solver::BaseLinearSolver*",
        "description": "The linear solver used to solve the system of equations arising from the implicit discretization."
      },
      {
        "name": "f_rayleighStiffness",
        "type": "Real",
        "defaultValue": "0.0",
        "description": "Rayleigh stiffness coefficient, which is part of the Rayleigh damping model that adds both mass-proportional and stiffness-proportional damping to the system.",
        "deprecated": true
      },
      {
        "name": "f_rayleighMass",
        "type": "Real",
        "defaultValue": "0.0",
        "description": "Rayleigh mass coefficient, another component of the Rayleigh damping model.",
        "deprecated": true
      },
      {
        "name": "rayleighStiffness",
        "type": "Real",
        "defaultValue": "0.0",
        "description": "Deprecation placeholder for the Rayleigh stiffness coefficient."
      },
      {
        "name": "rayleighMass",
        "type": "Real",
        "defaultValue": "0.0",
        "description": "Deprecation placeholder for the Rayleigh mass coefficient."
      },
      {
        "name": "f_printLog",
        "type": "bool",
        "defaultValue": "false",
        "description": "Boolean flag to enable or disable logging of position, velocity, and force at each time step."
      },
      {
        "name": "velocityDamping",
        "type": "Real",
        "defaultValue": "0.0",
        "description": "Coefficient for exponential damping applied to the velocities after integration.",
        "deprecated": true
      },
      {
        "name": "f_velocityDamping",
        "type": "Real",
        "defaultValue": "0.0",
        "description": "Deprecation placeholder for the velocity damping coefficient."
      },
      {
        "name": "rayleigh",
        "type": "core::behavior::Rayleigh*",
        "description": "Deprecated pointer to a Rayleigh object used in previous versions for defining Rayleigh damping.",
        "deprecated": true
      }
    ],
    "methods": [
      {
        "name": "getPositionIntegrationFactor()",
        "returnType": "SReal",
        "parameters": [],
        "description": "Returns the integration factor for position based on the current time step."
      },
      {
        "name": "getIntegrationFactor(inputDerivative, outputDerivative)",
        "returnType": "SReal",
        "parameters": [
          {
            "type": "int",
            "name": "inputDerivative"
          },
          {
            "type": "int",
            "name": "outputDerivative"
          }
        ],
        "description": "Returns the integration factor for converting derivatives from one order to another based on the current time step."
      },
      {
        "name": "getSolutionIntegrationFactor(outputDerivative)",
        "returnType": "SReal",
        "parameters": [
          {
            "type": "int",
            "name": "outputDerivative"
          }
        ],
        "description": "Returns the integration factor applied to the solution vector for a given derivative order based on the current time step."
      }
    ]
  },
  "maths": {
    "maths": "The Euler Implicit Solver is a numerical method used for solving ordinary differential equations (ODEs) that arise in the simulation of dynamic systems. It employs the implicit backward Euler scheme, which provides unconditional stability and hence is particularly useful for stiff problems where explicit methods can suffer from severe stability constraints.\n\n### Mathematical Description:\n\nConsider a system described by an ODE of the form:\n\\[ \\frac{d\textbf{x}}{dt} = \textbf{f}(t, \textbf{x}), \\]\nwhere \\(\textbf{x}\\) is the state vector representing positions and velocities, and \\(\textbf{f}\\) is a function that describes how the system evolves over time.\n\nThe Euler Implicit Solver approximates this ODE using the implicit backward Euler method. At each time step \\(t_n = n \times dt\\), where \\(dt\\) is the time step size, the state at the next time step \\(\textbf{x}_{n+1}\\) is determined by solving:\n\\[ \\frac{\textbf{x}_{n+1} - \textbf{x}_n}{dt} = \textbf{f}(t_{n+1}, \textbf{x}_{n+1}). \\]\n\nThis can be rewritten as a fixed point problem:\n\\[ \textbf{x}_{n+1} = \textbf{x}_n + dt \times \textbf{f}(t_{n+1}, \textbf{x}_{n+1}). \\]\n\nIn the context of mechanics, this ODE often takes the form of Newton's second law:\n\\[ M\frac{d^2\textbf{x}}{dt^2} + B\frac{d\textbf{x}}{dt} + K\textbf{x} = \textbf{F}(t), \\]\nwhere \\(M\\) is the mass matrix, \\(B\\) and \\(K\\) are damping and stiffness matrices respectively, and \\(\textbf{F}(t)\\) represents external forces.\n\nUsing the implicit backward Euler method, this becomes:\n\\[ M\frac{\textbf{x}_{n+1} - 2\textbf{x}_n + \textbf{x}_{n-1}}{dt^2} + B\frac{\textbf{x}_{n+1} - \textbf{x}_n}{dt} + K\textbf{x}_{n+1} = \textbf{F}(t_{n+1}). \\]\n\nThis leads to the following linear system at each time step:\n\\[ (M/dt^2 + B/dt + K) \textbf{x}_{n+1} = M\frac{\textbf{x}_n - 2\textbf{x}_{n-1}}{dt^2} - B\frac{\textbf{x}_n}{dt} + \textbf{F}(t_{n+1}). \\]\n\n### Physical Description:\n\nIn physical terms, the implicit backward Euler method ensures that both velocity and position are updated in a way that accounts for forces at the new time step. This makes it suitable for simulating systems with strong damping or stiffness, as well as systems where external forces change rapidly.\n\n#### Integration Factors:\n- **Position**: The position update is given by \\(\textbf{x}_{n+1} = \textbf{x}_n + dt \times \textbf{v}_{n+1}\\), where \\(\textbf{v}_{n+1}\\) is the velocity at time step \\(t_{n+1}\\).\n- **Velocity**: The velocity update is given by solving the system \\((M/dt^2 + B/dt + K)\textbf{x}_{n+1} = \textbf{RHS}(dt, t_n, \textbf{x}_n)\\), where RHS denotes the right-hand side of the equation.\n\nThe solver ensures that all forces, including constraints and external forces like gravity, are accounted for at each step. Additionally, it supports velocity damping through an exponential decay factor to simulate energy dissipation over time."
  },
  "summary": {
    "abstract": "The EulerImplicitSolver implements an implicit backward Euler scheme for solving first- and second-order ODEs in SOFA simulations, supporting Rayleigh damping and trapezoidal rule options.",
    "sheet": "# EulerImplicitSolver\n\n## Overview\n\nEulerImplicitSolver is a time integrator that uses the implicit backward Euler scheme to solve ordinary differential equations (ODEs) of first or second order. It supports Rayleigh damping and can optionally use the trapezoidal rule for higher accuracy.\n\n## Mathematical Model\n\nFor a system described by an ODE, the implicit backward Euler method approximates the state at each time step as follows:\n\n- **Second Order**\n  \\[ x_{t+h} = x_t + h v_{t+h}, \\]\n  \\[ v_{t+h} = v_t + h a_{t+h}. \\]\n  The unknown is \\( dv = v_{t+h} - v_t \\). Newton's law yields:\n  \\[ M dv = h f(t+h), \\]\n  where \\( M \\) is the mass matrix, and \\( K = df/dx \\) and \\( B = df/dv \\) are stiffness and damping matrices, respectively. Rayleigh damping adds terms \\( -r_M M + r_K K \\). The final equation system is:\n  \\[ P ( (1+h r_M) M - h B - h(h + r_K) K ) P dv = P h ( f(t) + (h+r_K) K v + B v - r_M M v ). \\]\n\n- **First Order**\n  For a first-order system, the integration scheme is:\n  \\[ x_{t+h} = x_t + h v_{t+h}, \\]\n  where \\( M v_{t+h} = f_{ext}(t) - h K v_{t+h} \\), leading to:\n  \\[ (M + h K) v_{t+h} = f_{ext}. \\]\n\n- **Trapezoidal Rule**\n  For the trapezoidal scheme, the equation for second-order systems becomes:\n  \\[ P ( (1+h/2 r_M) M - h/2 B - h/2 (h + r_K) K ) P dv = P h/2 ( 2 f(t) + (h + r_K) K v + B v - r_M M v ). \\]\n\n## Parameters and Data\n\n- **rayleighStiffness**: Rayleigh damping coefficient related to stiffness, > 0.\n- **rayleighMass**: Rayleigh damping coefficient related to mass, > 0.\n- **vdamping**: Velocity decay coefficient (no decay if null).\n- **firstOrder**: Use backward Euler scheme for first-order ODE system (default: false).\n- **trapezoidalScheme**: Boolean to use the trapezoidal scheme instead of implicit Euler (default: false).\n- **solveConstraint**: Apply ConstraintSolver (requires a ConstraintSolver in the same node, default: false).\n- **threadSafeVisitor**: If true, do not use realloc and free visitors in fwdInteractionForceField.\n- **computeResidual**: If true, compute residual at the end of solving.\n\n## Practical Notes\n\nThe EulerImplicitSolver is unconditionally stable but may require careful tuning of parameters such as Rayleigh damping coefficients. The trapezoidal scheme can improve accuracy but increases computational cost."
  }
}