Incremental Potential and Dissipation

Dissipative Systems

Hamilton’s least action principle always gives us conservative forces.

$${\left({\int }_0^T\left(K(\mathbf{q},\dot{\mathbf{q}})-V(\mathbf{q})\right)dt\right)}^{\circ }=0⟹\underset{\begin{array}{c}\text{change of}\\ \text{momentum}\end{array}}{\underbrace{\frac{d}{dt}\frac{\partial K}{\partial \dot{\mathbf{q}}}}}\underset{\begin{array}{c}\text{fictitious}\\ \text{force}\end{array}}{\underbrace{-\frac{\partial K}{\partial \mathbf{q}}}}=\underset{\begin{array}{c}\text{conservative}\\ \text{force}\end{array}}{\underbrace{-\frac{\partial V}{\partial \mathbf{q}}}}$$

This works for conservative systems, but what happens when we want to simulate a block sliding to a halt due to friction? The standard action functional doesn’t work, because the friction force is non-conservative.

Traditionally, D’Alembert’s Principle is a generalization of the least action that allows any force, including non-conservative forces like friction. The idea is to add “virtual work” done from any additional force. So,

$${\left({\int }_0^T\left(K(\mathbf{q},\dot{\mathbf{q}})-V(\mathbf{q})\right)dt\right)}^{\circ }=\underset{\begin{array}{c}\text{virtual work from}\\ \text{any additional}\\ \text{forces}\end{array}}{\underbrace{{\int }_0^T⟨{\mathbf{f}}_{nc}(\mathbf{q},\dot{\mathbf{q}})\mid \mathring{\mathbf{q}}⟩dt}}⟹\frac{d}{dt}\frac{\partial K}{\partial \dot{\mathbf{q}}}-\frac{\partial K}{\partial \mathbf{q}}=-\frac{\partial V}{\partial \mathbf{q}}+{\mathbf{f}}_{nc}$$

Definition (Rayleigh Dissipation Function)

Can’t dissipative forces be written as a derivative of some function? In 1873 Rayleigh observed the following: Suppose we have some “linear viscous”¹ force.

$$m\ddot{\mathbf{q}}=\underset{\begin{array}{c}\text{conservative}\\ \text{force}\end{array}}{\underbrace{-\frac{\partial V}{\partial \mathbf{q}}}}\underset{\begin{array}{c}\text{damping force}\\ \text{linear in}\\ \text{velocity}\end{array}}{\underbrace{-\mathbf{B}\dot{\mathbf{q}}}}$$

This linear damping is actually the variation of some (quadratic) function:

$$\mathbf{B}\dot{\mathbf{q}}=\frac{\partial }{\partial \dot{\mathbf{q}}}R(\dot{\mathbf{q}})$$

where

$$R(\dot{\mathbf{q}})=\frac{1}{2}{\dot{\mathbf{q}}}^{\top }\mathbf{B}\dot{\mathbf{q}}$$

is called the Rayleigh Dissipation Function. $\mathbf{B}$ here is a constant number (or matrix). The thicker the fluid, the “larger” $\mathbf{B}$. Indeed, the force depends on the velocity $\dot{\mathbf{q}}$. Here we defined $R$ via the power rule.

Rayleigh’s Euler-Lagrange Equation

We can derive the equations of motion for a system with Rayleigh dissipation by modifying the least action principle. Given kinetic energy $K(\mathbf{q},\dot{\mathbf{q}})$, potential energy $V(\mathbf{q})$, and a dissipation function $R(\mathbf{q},\dot{\mathbf{q}})$, we have

$$\underset{\begin{array}{c}\text{change of}\\ \text{momentum}\end{array}}{\underbrace{\frac{d}{dt}\frac{\partial K}{\partial \dot{\mathbf{q}}}}}\underset{\begin{array}{c}\text{fictitious}\\ \text{force}\end{array}}{\underbrace{-\frac{\partial K}{\partial \mathbf{q}}}}=\underset{\begin{array}{c}\text{conservative}\\ \text{force}\end{array}}{\underbrace{-\frac{\partial V}{\partial \mathbf{q}}}}\underset{\begin{array}{c}\text{dissipative}\\ \text{force}\end{array}}{\underbrace{-\frac{\partial R}{\partial \dot{\mathbf{q}}}}}$$

with respect to velocity. This construction accurately reflects reality if the mapping $\dot{\mathbf{q}}\mapsto R(\mathbf{q},\dot{\mathbf{q}})$ satsifies the following properties:

1. The mapping must be convex (second derivative is nonnegative). As velocity increases, the dissipative force must reliably increase or remain constant.
2. It must be zero at $\dot{\mathbf{q}}=0$. The dissipative force only acts when there is some velocity.
3. It must be nonnegative otherwise. This ensures adherence to the second law of thermodynamics. The dissipative force should always be doing negative work, so the system should be losing energy.

Related Works

• Helmholtz’s minimal dissipation principle
• Originally, Rayleigh only considered quadratic $R$.
• After 1970’s (Moreau) discovered that other non-quadratic $R$ can recover nonlinear forces like Coulomb friction and plasticity.
• There is a natural time discretization (Kane, Marsden, Ortiz, West 1999).

Remarks

We can apply Noether’s theorem for time independence

$$\frac{d}{dt}\left(\frac{\partial K}{\partial \dot{\mathbf{q}}}\dot{\mathbf{q}}-K+V\right)=-\frac{\partial R}{\partial \dot{\mathbf{q}}}\dot{\mathbf{q}}$$

If $K$ and $R$ are quadratic (highest order of the variable in exactly $2$) in $\dot{\mathbf{q}}$, then

$$\frac{\partial K}{\partial \dot{\mathbf{q}}}\dot{\mathbf{q}}=2K\text{and}\frac{\partial R}{\partial \dot{\mathbf{q}}}\dot{\mathbf{q}}=2R$$

such that

$$\frac{d}{dt}\left(K+V\right)=-2R$$

To design dissipation, just model the rate of energy dissipation.

Interlude (Equations of Motion for Discrete Time Dissipation)

Let’s try dropping the continuous time derivative and model the discrete time energy loss. For simplicity, assume $K(\mathbf{q},\dot{\mathbf{q}})=\frac{1}{2}{\dot{\mathbf{q}}}^{\top }\mathbf{M}\dot{\mathbf{q}}$. The equation of motion is

$$\mathbf{M}\ddot{\mathbf{q}}=-\frac{\partial V}{\partial \mathbf{q}}-\frac{\partial R}{\partial \dot{\mathbf{q}}}$$

Suppose the space of positions is given by $Q={\mathbb{R}}^m$ where each point has coordinate $\mathbf{q}=(q_1,\dots ,q_m{)}^{\top }$. Suppose the inertia is independent of position $\mathbf{q}$ and

$$\text{KineticEnergy}(\dot{\mathbf{q}})=\frac{1}{2}{\dot{\mathbf{q}}}^{\top }\mathbf{M}\dot{\mathbf{q}}$$

and the potential energy is $U=U(q_1,\dots ,q_m)$. Then the equation of motion is

$$(\mathbf{M}\ddot{\mathbf{q}}{)}_i=-(dU{)}_i=-\frac{\partial U}{\partial q_i}$$

Symplectic & Backward Euler with Dissipation

We can discretize time $t^{(n)}=n\Delta t$. Denote the state at the $n-$th step as ${\mathbf{q}}^{(n)}$. We can approximate the second time derivative as

$$(\ddot{\mathbf{q}}{)}^{(n)}\approx \frac{1}{\Delta t^2}\left({\mathbf{q}}^{(n-1)}-2{\mathbf{q}}^{(n)}+{\mathbf{q}}^{(n+1)}\right)$$

Given ${\mathbf{q}}^{(n-1)},{\mathbf{q}}^{(n)}$, solve for ${\mathbf{q}}^{(n+1)}$. We use the Euler method.

Symplectic (explicit):

$$\frac{1}{\Delta t^2}\left({\mathbf{q}}^{(n-1)}-2{\mathbf{q}}^{(n)}+{\mathbf{q}}^{(n+1)}\right)=-{\mathbf{M}}^{-1}(dU){|}_{{\mathbf{q}}^{(n)}}$$

Backward (implicit):

$$\frac{1}{\Delta t^2}\left({\mathbf{q}}^{(n-1)}-2{\mathbf{q}}^{(n)}+{\mathbf{q}}^{(n+1)}\right)=-{\mathbf{M}}^{-1}(dU){|}_{{\mathbf{q}}^{(n+1)}}$$

Example 1: Damped Harmonic Oscillator

We can try to simulate a damped spring. Recall the equations of motion for a spring and its solved form:

$$\{\begin{array}{l}\ddot{q}+{\omega }^{2}q=0\\ q(t)=a\cos (\omega t)+b\sin (\omega t)\end{array}$$

From interlude, ${\mathbf{M}}^{-1}dU{|}_{\mathbf{q}}={\omega }^2\mathbf{q}$. The symplectic form is

$$\begin{array}{rl}\frac{1}{\Delta t^2}\left({\mathbf{q}}^{(n-1)}-2{\mathbf{q}}^{(n)}+{\mathbf{q}}^{(n+1)}\right)& =-{\mathbf{M}}^{-1}(dU){|}_{{\mathbf{q}}^{(n)}}\\ {\mathbf{q}}^{(n-1)}-2{\mathbf{q}}^{(n)}+{\mathbf{q}}^{(n+1)}& =-\Delta t^2\left({\omega }^2{\mathbf{q}}^{(n)}\right)\\ {\mathbf{q}}^{(n+1)}& =-\Delta t^2\left({\omega }^2{\mathbf{q}}^{(n)}\right)+2{\mathbf{q}}^{(n)}-{\mathbf{q}}^{(n-1)}\\ {\mathbf{q}}^{(n+1)}& =-{\mathbf{q}}^{(n-1)}+\left(2-\Delta t^2{\omega }^2\right){\mathbf{q}}^{(n)}\end{array}$$

We can express this as a matrix multiplication:

$$\left[\begin{array}{c}{\mathbf{q}}^{(n)}\\ {\mathbf{q}}^{(n+1)}\end{array}\right]=\left[\begin{array}{cc}0& 1\\ -1& 2-\Delta t^2{\omega }^2\end{array}\right]\left[\begin{array}{c}{\mathbf{q}}^{(n-1)}\\ {\mathbf{q}}^{(n)}\end{array}\right]$$

The determinant of the square matrix is $1$, so the symplectic form is volume-preserving (area). Both eigenvalue norms are $1$ when $\Delta t^2{\omega }^2<4$ (conditional stability).

The backward form is

$$\begin{array}{rl}\frac{1}{\Delta t^2}\left({\mathbf{q}}^{(n-1)}-2{\mathbf{q}}^{(n)}+{\mathbf{q}}^{(n+1)}\right)& =-{\mathbf{M}}^{-1}(dU){|}_{{\mathbf{q}}^{(n+1)}}\\ {\mathbf{q}}^{(n-1)}-2{\mathbf{q}}^{(n)}+{\mathbf{q}}^{(n+1)}& =-\Delta t^2\left({\omega }^2{\mathbf{q}}^{(n+1)}\right)\\ {\mathbf{q}}^{(n+1)}& =-{\mathbf{q}}^{(n-1)}+2{\mathbf{q}}^{(n)}-\Delta t^2\left({\omega }^2{\mathbf{q}}^{(n+1)}\right)\\ {\mathbf{q}}^{(n+1)}& =\frac{-{\mathbf{q}}^{(n-1)}+2{\mathbf{q}}^{(n)}}{1+\Delta t^2{\omega }^2}\end{array}$$

and likewise,

$$\left[\begin{array}{c}{\mathbf{q}}^{(n)}\\ {\mathbf{q}}^{(n+1)}\end{array}\right]=\left[\begin{array}{cc}0& 1\\ -\frac{1}{1+\Delta t^2{\omega }^2}& \frac{2}{1+\Delta t^2{\omega }^2}\end{array}\right]\left[\begin{array}{c}{\mathbf{q}}^{(n-1)}\\ {\mathbf{q}}^{(n)}\end{array}\right]$$

The determinant of the square matrix is $<1$ and so the area shrinks. All eigenvalue norms are $<1$ for all $\Delta t^2{\omega }^2$ (unconditional stability).

Minimal Incremental Potential Principle

We can exploit the unconditional stability of the backward Euler method to design a variational principle for discrete time dissipation. Recall the backward Euler method:

$$\frac{1}{\Delta t^2}\left({\mathbf{q}}^{(n-1)}-2{\mathbf{q}}^{(n)}+{\mathbf{q}}^{(n+1)}\right)=-(dU){|}_{{\mathbf{q}}^{(n+1)}}$$

Let’s rearrange with ${\mathbf{q}}_{\text{pred}}:=2{\mathbf{q}}^{(n)}-{\mathbf{q}}^{(n-1)}$.

$$\frac{1}{\Delta t^2}\left({\mathbf{q}}^{(n+1)}-{\mathbf{q}}_{\text{pred}}\right)+(dU){|}_{{\mathbf{q}}^{(n+1)}}=0$$

This is actually an optimality condition for the following optimization problem:

$${\mathbf{q}}^{(n+1)}=\underset{\mathbf{q}\in {\mathbb{R}}^m}{\mathrm{argmin}}\left(\frac{1}{2\Delta t^2}\Vert \mathbf{q}-{\mathbf{q}}_{\text{pred}}{\Vert }_{\mathbf{M}}^2+U(\mathbf{q})\right)$$

Sneakily, $\text{pred}$ stands for “predicted position”.

$"\\usepackage{amsmath}\n\\usepackage{amssymb}\n\\usepackage{amsfonts}\n\\usepackage{xcolor}\n\n\\begin{document}\n\\begin{tikzpicture}[>=stealth, scale=1]\n\n$

source code

In terms of velocity, we can write

$$\begin{array}{rl}{\mathbf{v}}_{\text{new}}& =\frac{1}{\Delta t}\left({\mathbf{q}}^{(n+1)}-{\mathbf{q}}^{(n)}\right)\\ {\mathbf{v}}_{\text{old}}& =\frac{1}{\Delta t}\left({\mathbf{q}}^{(n)}-{\mathbf{q}}^{(n-1)}\right)\end{array}$$

with the following visualization:

$"\\usepackage{amsmath}\n\\usepackage{amssymb}\n\\usepackage{amsfonts}\n\\usepackage{xcolor}\n\n\\begin{document}\n\\begin{tikzpicture}[>=stealth, scale=1]\n\n$

source code

So, the physical system will decide its new velocity by

$${\mathbf{v}}_{\text{new}}=\underset{\mathbf{v}\in {\mathbb{R}}^m}{\mathrm{argmin}}\left(\frac{1}{2}\Vert \mathbf{v}-{\mathbf{v}}_{\text{old}}{\Vert }_{\mathbf{M}}^2+U\left({\mathbf{q}}^{(n)}+\Delta t\mathbf{v}\right)\right)$$

This equation is called the minimal incremental potential principle. In particular, the terms compete with one another:

1. The intertia term. It encourages the new velocity to be close to the old velocity. So any deviation from ${\mathbf{v}}_{\text{old}}$ is penalized, scaled by the mass matrix $\mathbf{M}$.
2. The potential term. The object is drawn toward lower energy states like falling under gravity or a spring returning to its rest position.

The true physical path is simply the exact velocity vector that minimizes the sum of these two penalties. Indeed, this means calculating the velocity for every time step requires a good numerical optimizer.

For collision and contact, we can add “smooth barrier” functions in potential and perform optimization properly.

$$b(d,\hat{d})=\{\begin{array}{ll}-{\left(d-\hat{d}\right)}^2\ln \left(\frac{d}{\hat{d}}\right)& \text{if }d<\hat{d}\\ 0& \text{if }d\ge \hat{d}\end{array}$$

$"\\usepackage{amsmath}\n\\usepackage{amssymb}\n\\usepackage{xcolor}\n\n\\begin{document}\n\\begin{tikzpicture}[x=4.5cm, y=1.6cm, >=stealth]\n\n$

source code

Adding Viscous Dissipation

Recall the backward Euler update on conservative systems:

$${\mathbf{q}}^{(n+1)}=\underset{\mathbf{q}\in {\mathbb{R}}^m}{\mathrm{argmin}}\left(\frac{1}{2\Delta t^2}\Vert \mathbf{q}-{\mathbf{q}}_{\text{pred}}{\Vert }_{\mathbf{M}}^2+U(\mathbf{q})\right)$$

where ${\mathbf{q}}_{\text{pred}}=2{\mathbf{q}}^{(n)}-{\mathbf{q}}^{(n-1)}$. We can add dissipation by adding a Rayleigh dissipation funtion:

$${\mathbf{q}}^{(n+1)}=\underset{\mathbf{q}\in {\mathbb{R}}^m}{\mathrm{argmin}}\left(\frac{1}{2\Delta t^2}\Vert \mathbf{q}-{\mathbf{q}}_{\text{pred}}{\Vert }_{\mathbf{M}}^2+\underset{\text{Effective Incremental Potential}}{\underbrace{U(\mathbf{q})+(\Delta t)R\left({\mathbf{q}}^{(n)},\frac{\mathbf{q}-{\mathbf{q}}^{(n)}}{\Delta t}\right)}}\right)$$

Equivalently,

$${\mathbf{v}}_{\text{new}}=\underset{\mathbf{v}\in {\mathbb{R}}^m}{\mathrm{argmin}}\left(\frac{1}{2}\Vert \mathbf{v}-{\mathbf{v}}_{\text{old}}{\Vert }_{\mathbf{M}}^2+\underset{\text{Effective Incremental Potential}}{\underbrace{U\left({\mathbf{q}}^{(n)}+\Delta t\mathbf{v}\right)+(\Delta t)R\left({\mathbf{q}}^{(n)},\mathbf{v}\right)}}\right)$$

The third term also competes with the first two terms. It is the dissipation penalty. The object wants to minimize the energy it bleeds out into the environment (i.e. heat from friction). Because $R$ is zero when velocity is zero, this term effectively pulls the optimal velocity closer to zero.

Writing the system in $F=ma$,

$$\mathbf{M}\frac{{\mathbf{v}}_{\text{new}}-{\mathbf{v}}_{\text{old}}}{\Delta t}=-\frac{\partial U}{\partial \mathbf{q}}-\frac{\partial R}{\partial \mathbf{v}}$$

Places to use a quadratic dissipation function

1. Car’s shock absorber. It is a piston pushing through oil (or a viscous fluid). It restructs motion proportional to velocity.
2. Lubricated sliding. Two metal plates separated by a thin layer of oil. It creates a smooth, linear drag force.

Quasi-Static System

To study a dissipative system, we often consdier a quasi-static regime, where inertia is negligible. In this case, the system is in balance with potential force and external force at all times. So,

$$\begin{array}{rl}\cancel{\bcancel{\mathbf{M}\frac{{\mathbf{v}}_{\text{new}}-{\mathbf{v}}_{\text{old}}}{\Delta t}}}& =-\frac{\partial U}{\partial \mathbf{q}}-\frac{\partial R}{\partial \mathbf{v}}\\ \frac{\partial R}{\partial \mathbf{v}}& =-\frac{\partial U}{\partial \mathbf{q}}+{\mathbf{f}}_{\text{ext}}\end{array}$$

For quadratic $R$, this determines a “terminal velocity” that the system will eventually reach.

Intuition

We can think of dragging a heavy block through thick mud (or an extremely viscous fluid). The block never really accelerates nor does it carry momentum. It’s speed is dictated by how hard it is pulled against the mud’s resistance. When the block is no longer being pulled, it stops.

This is as if the mass $\mathbf{M}$ approaches zero (because we basically have zero momentum). This is why it is “quasi-static”, because it’s kind of still (static) but not quite (quasi).

Importantly, the crossed out penalty term vanishes.

Traditionally, for solving a general force, we solve some kind of relation between $\mathbf{f},\mathbf{q},\mathbf{v}$. Here, ${\mathbf{f}}_{\text{ext}}$ is an ireversible force that is not a function of $\mathbf{q}$.

Dry Friction

$"\\usepackage{amsmath}\n\\usepackage{amssymb}\n\\usepackage{xcolor}\n\n\\begin{document}\n\\begin{tikzpicture}[>=stealth, scale=1.3]\n\n$

source code

We have some laws of friction:

• Amonton’s 1st law: friction force is proportional to the normal force
• Amonton’s 2nd law: friction force is independent of the contact area
• Coulomb’s law: Once the motion starts, the friction force is independent of velocity.

$"\\usepackage{amsmath}\n\\usepackage{amssymb}\n\\usepackage{amsfonts}\n\\usepackage{xcolor}\n\n\\begin{document}\n\\begin{tikzpicture}[x=1cm, y=1cm, scale=0.8, >=stealth]\n\n$

source code

The force ${\mathbf{f}}_c$ at contact lies in a friction cone (in the dual space at contact, or the tangent cone).

• At each point of contact, we have an outward normal (covector) $\mathbf{n}$.
• The relative velocity between contact should satisfy $⟨\mathbf{n}|\mathbf{v}⟩\ge 0$.
• The normal (${\mathbf{f}}^{\perp }$) and tangent part (${\mathbf{f}}^{||}$) of the contact force, and tangent velocity (${\mathbf{v}}^{||}$) should satisfy the following relations: $$|{\mathbf{f}}^{||}|\le \mu {\mathbf{f}}^{\perp }|{\mathbf{f}}^{||}|<\mu {\mathbf{f}}^{\perp }⟺{\mathbf{v}}^{||}=0$$ where $\mu$ is the coefficient of friction. In the image above, $\mu$ dictates the width of the friction cone.
• When the tangent velocity is nonzero, the tangent force is in the same direction with it $\alpha {\mathbf{f}}^{||}={♭}_{{\mathbb{R}}^3}\mathbf{v},\alpha \ge 0$

From physics, these two relations are merely:

1. Static Friction (stuck): If the required tangential force (the force pushing the object to the right) is strictly less than the boundary of the cone, the object is perfectly stuck and thus has zero tangential velocity.
2. Kinetic Friction (sliding): If the external push exceeds the cone’s edge, the friction force maxes out at the cone boundary $|{\mathbf{f}}^{||}|=\mu |{\mathbf{f}}^{\perp }|$. The object begins to slide and the friction force points in the exact opposite direction of the motion.

Classical Approach for Contact

Suppose we want to model how a rigid body comes into contact with a wall or some other object. We first establish a set of points ${\mathbf{p}}_i$ that will be in contact. Then

$"\\usepackage{amsmath}\n\\usepackage{amssymb}\n\\usepackage{amsfonts}\n\\usepackage{xcolor}\n\n\\begin{document}\n\\begin{tikzpicture}[>=stealth, scale=1.0]\n\n$

source code

This is computationally miserable. When a 3D object hits the floor, we have to keep track of the object’s state $(\mathbf{c},\mathbf{R})\in Q$ and the contact space where the objects touch $({\mathbf{p}}_1,\dots ,{\mathbf{p}}_k)\in {\mathbb{R}}^3\times \cdots \times {\mathbb{R}}^3$. To know if an object is sliding or stuck, we need to view the object’s linear and angular velocity $(\dot{\mathbf{c}},\dot{\mathbf{R}})$ and map that to 3D velocities at every single contact point $({\dot{\mathbf{p}}}_1,\dots ,{\dot{\mathbf{p}}}_k)$.

This translation is done via the Jacobian denoted by $d\varphi$. Conversely, the floor pushes back on those chair legs by friction and normal forces $({\mathbf{f}}_1,\dots ,{\mathbf{f}}_k)$. To know how those forces affect the chair’s motion, we need to map the net force and net torque back to the cneter of pass. This is done via the dual (tranpose, since ${\mathbb{R}}^3$) of the Jacobian, denoted by $d{\varphi }^{\ast }$.

We then need to solve for velocity and contact forces so that all conditions are satisfied. Because the body is rigid, all the particles are connected.

1. The friction force at point ${\mathbf{p}}_1$ depends on its velocity.
2. But its velocity depends on the total rotation of the chair.
3. The total rotation of the chair depends on the forces pushing up on points ${\mathbf{p}}_2,{\mathbf{p}}_3,{\mathbf{p}}_4$.

Obviously, we cannot calculate them all in closed form. This is called the Linear Complementarity Problem (LCP). The solver would have to guess if each point is stuck, sliding, or lifting off the ground.

Dissipation Function for Dry Friction

We can take a variational approach to dry friction.

$$\begin{array}{rl}\mathbf{M}\frac{{\mathbf{v}}_{\text{new}}-{\mathbf{v}}_{\text{old}}}{\Delta t}& =-\frac{\partial U}{\partial \mathbf{q}}-\frac{\partial R}{\partial \mathbf{v}}\\ R(\mathbf{v})& =\mu |\mathbf{v}|\end{array}$$

Newmark Algorithm

We initally chose the Backward Euler method for its simplicity and unconditional stability. However, it is only “first order” accurate, meaning it introduces a lot of numerical error over time.

Recall the backward Euler update:

$${\mathbf{q}}^{(n+1)}=\underset{\mathbf{q}\in {\mathbb{R}}^m}{\mathrm{argmin}}\left(\frac{1}{2\Delta t^2}\Vert \mathbf{q}-{\mathbf{q}}_{\text{pred}}{\Vert }_{\mathbf{M}}^2+U(\mathbf{q})+\Delta tR\left({\mathbf{q}}^{(n)},\frac{\mathbf{q}-{\mathbf{q}}^{(n)}}{\Delta t}\right)\right)$$

where ${\mathbf{q}}_{\text{pred}}:=2{\mathbf{q}}^{(n)}-{\mathbf{q}}^{(n-1)}$. The ${\mathbf{q}}_{\text{pred}}$ assumes the object will continue moving at the same velocity as before, ignoring forces and acceleration.

Now, let the incremental potential be

$${\tilde{U}}^{(n)}\left({\mathbf{q}}^{(n+1)}\right)=U\left({\mathbf{q}}^{(n+1)}\right)+\Delta tR\left((1-s){\mathbf{q}}^{(n)}+s{\mathbf{q}}^{(n+1)},\frac{{\mathbf{q}}^{(n+1)}-{\mathbf{q}}^{(n)}}{\Delta t}\right)$$

The acceleration at current time can be read off from

$${\mathbf{a}}^{(n)}={\mathbf{M}}^{-1}\frac{\partial {\tilde{U}}^{(n)}\left({\mathbf{q}}^{(n)}\right)}{\partial {\mathbf{q}}^{(n)}}$$

Velocity can be kept track of by

$${\mathbf{v}}^{(n)}={\mathbf{v}}^{(n-1)}+\Delta t\left((1-\gamma ){\mathbf{a}}^{(n-1)}+\gamma {\mathbf{a}}^{(n)}\right)$$

$\gamma$ is a parameter that controls how much we trust the current acceleration ${\mathbf{a}}^{(n)}$ vs the previous acceleration ${\mathbf{a}}^{(n-1)}$. The prediction of the next time step

$${\mathbf{q}}_{\text{pred}}^{(n+1)}:={\mathbf{q}}^{(n)}+\Delta t{\mathbf{v}}^{(n)}+\frac{\Delta t^2}{2}(1-2\beta ){\mathbf{a}}^{(n)}$$

If we ignore the $\beta$ parameter for a moment, this prediction is actually the classic kinematic equation from basic physiscs, $x=x_0+v+\frac{1}{2}a{t}^2$. We are basically predicting the next position using the current velocity and current acceleration.

We then solve an optimization problem like before,

$${\mathbf{q}}^{(n+1)}=\underset{\mathbf{q}\in {\mathbb{R}}^m}{\mathrm{argmin}}\left(\frac{1}{2\Delta t^2}\Vert \mathbf{q}-{\mathbf{q}}_{\text{pred}}{\Vert }_{\mathbf{M}}^2+\beta {\tilde{U}}^{(n)}(\mathbf{q})\right)$$

This is the Newmark algorithm.

Numerical Optimization for Contact

Since we have a variational formulation for contact and have shown that is simply an optimization problem, we can solve for contact forces and velocities using numerical optimization. This is a more direct approach than the classical LCP formulation.

Using the smooth barrier and dissipation function, every time step is actually one unconstrained optimization problem

$$\underset{\mathbf{x}\in {\mathbb{R}}^m}{\mathrm{minimize}}\mathcal{L}(\mathbf{x})$$

Here, $\mathbf{x}$ may be velocity or position

$$\begin{array}{rl}{\mathbf{q}}^{(n+1)}& =\underset{\mathbf{q}\in {\mathbb{R}}^m}{\mathrm{argmin}}\left(\frac{1}{2\Delta t^2}\Vert \mathbf{q}-{\mathbf{q}}_{\text{pred}}{\Vert }_{\mathbf{M}}^2+U(\mathbf{q})\right)\\ {\mathbf{v}}_{\text{new}}& =\underset{\mathbf{v}\in {\mathbb{R}}^m}{\mathrm{argmin}}\left(\frac{1}{2}\Vert \mathbf{v}-{\mathbf{v}}_{\text{old}}{\Vert }_{\mathbf{M}}^2+U\left({\mathbf{q}}^{(n)}+\Delta t\mathbf{v}\right)\right)\end{array}$$

Note that the initial guess (from the state of the previous time frame) for optimization is usually very close to the optimizer. We can use gradient descent using some metric $♭=\mathbf{H}$.

$${\mathbf{x}}^{(n+1)}\leftarrow {\mathbf{x}}^{(n)}-\alpha ♯(d\mathcal{L}{)}_{{\mathbf{x}}^{(n)}}={\mathbf{x}}^{(n)}-\alpha {\mathbf{H}}^{-1}\left[\begin{array}{c}\partial \mathcal{L}/\partial x_1\\ ⋮\\ \partial \mathcal{L}/\partial x_m\end{array}\right]$$

Now, we just need to choose a good $\mathbf{H}$ and a step size $\alpha >0$.

• Classic gradient descent uses $\mathbf{H}=\mathbf{I}$ and a small step size $\alpha$. This is very stable but can be slow to converge.
• Newton’s method uses $\mathbf{H}=\text{Hessian}(\mathcal{L})$.
• Quasi-Newton methods use an approximation of the Hessian.
• Use line search for choosing $\alpha$.

Footnotes

1. Imagine running your hand through water. As you move faster, the force pushing back on your hand gets stronger. The same principle is true for air resistance and friction between surfaces. ↩