29.5. Energy conservation of integration schemes#

Energy conservation can be discussed with:

  • the Jacobian of the transformation induced by one-step of the method, \(\mathbf{z}_n \mapsto \mathbf{z}_{n+1}\)

  • the change of energy \(\Delta E = E_{n+1} - E_n\)

Energy conservation is discussed here for Hamiltonian systems, see simplectic integrators

\[\begin{split}\left\{\begin{aligned} \dot{\mathbf{q}} & = \partial_{\mathbf{p}} H \\ \dot{\mathbf{p}} & = - \partial_{\mathbf{q}} H \\ \end{aligned}\right.\end{split}\]

or for second-order mechanical systems

\[\mathbf{M} \ddot{\mathbf{x}} = \mathbf{f}(\mathbf{x},t)\]

with \(\mathbf{f} = - \partial_{\mathbf{x}} V(\mathbf{x})\).

The Hamiltonian function \(H(\mathbf{q}, \mathbf{p})\) can be thought as the energy as a function of the generalized coordinates \(\mathbf{q}\) and the generalized momenta \(\mathbf{p}\). For Hamiltoninan system, where the Hamiltonian function is independent from time, the Hamiltonian (the energy) is conserved along the trajectories of the system,

\[\begin{split}\begin{aligned} \dot{H} & = \dot{\mathbf{q}} \cdot \partial_{\mathbf{q}} H + \dot{\mathbf{p}} \cdot \partial_{\mathbf{p}} H = \\ & = \partial_{\mathbf{p}} H \cdot \partial_{\mathbf{q}} H - \partial_{\mathbf{q}} H \cdot \partial_{\mathbf{p}} H = 0 \ . \end{aligned}\end{split}\]

The second-order mechanical system is conservative if with \(\mathbf{f} = - \partial_{\mathbf{x}} V(\mathbf{x})\), as the total mechanical energy

\[E = \frac{1}{2} \dot{\mathbf{x}}^* \mathbf{M} \dot{\mathbf{x}} + V(\mathbf{x})\]

is conserved along the evolution of the system. This second-order system can be written - if \(\mathbf{M}\) is invertible - as a first-order system as

\[\begin{split}\dfrac{d}{dt} \begin{bmatrix} \mathbf{x} \\ \mathbf{v} \end{bmatrix} = \begin{bmatrix} \mathbf{v} \\ \mathbf{M}^{-1} \mathbf{f}(\mathbf{x},t) \end{bmatrix} \ ,\end{split}\]

or as a first-order system with Hamiltonian form, with \(H(\mathbf{p}, \mathbf{q}) = \frac{1}{2} \mathbf{p}^* \mathbf{M}^{-1} \mathbf{p} + V(\mathbf{q})\), having defined \(\mathbf{p} := \mathbf{M} \dot{\mathbf{x}}\), \(\mathbf{q}:= \mathbf{x}\),

\[\begin{split}\dfrac{d}{dt} \begin{bmatrix} \mathbf{q} \\ \mathbf{p} \end{bmatrix} = \begin{bmatrix} \mathbf{M}^{-1} \mathbf{p} \\ - \partial_{\mathbf{q}}V(\mathbf{q}) \end{bmatrix} \ .\end{split}\]

.

29.5.1. Examples#

29.5.1.1. Verlet#

For the mechanical system governed by the differential equation \(\ddot{\mathbf{x}} = \mathbf{f}(\mathbf{x},t)\), the numerical scheme reads

\[\begin{split}\begin{aligned} \widetilde{\mathbf{x}}_{n+1} & = \mathbf{x}_n + \frac{h}{2} \mathbf{v}_n \\ \mathbf{v}_{n+1} & = \mathbf{v}_n + h \mathbf{f}(\widetilde{\mathbf{x}}_{n+1})\\ \mathbf{x}_{n+1} & = \widetilde{\mathbf{x}}_{n+1} + \frac{h}{2} \mathbf{v}_{n+1} \end{aligned}\end{split}\]

so that the transformation \(\mathbf{z}_{n} \mapsto \mathbf{z}_{n+1}\) reads

\[\begin{split}\begin{aligned} \mathbf{v}_{n+1} & = \mathbf{v}_n + h \mathbf{f}(\widetilde{\mathbf{x}}_{n+1}(\mathbf{x}_n, \mathbf{v}_n) )\\ \mathbf{x}_{n+1} & = \widetilde{\mathbf{x}}_{n+1}(\mathbf{x}_n, \mathbf{v}_n) + \frac{h}{2} \mathbf{v}_{n+1} = \widetilde{\mathbf{x}}_{n+1}(\mathbf{x}_n, \mathbf{v}_n) + \frac{h}{2} \left[ \mathbf{v}_n + h \mathbf{f}(\widetilde{\mathbf{x}}_{n+1}(\mathbf{x}_n, \mathbf{v}_n) ) \right] \end{aligned}\end{split}\]

or

\[\begin{split}\begin{aligned} \mathbf{v}_{n+1} & = \mathbf{v}_n + h \mathbf{f}\left( \mathbf{x}_n + \frac{h}{2} \mathbf{v}_n \right) \\ \mathbf{x}_{n+1} & = \mathbf{x}_{n} + h \mathbf{v}_n + \frac{h^2}{2} \mathbf{f}\left( \mathbf{x}_n + \frac{h}{2} \mathbf{v}_n \right) \end{aligned}\end{split}\]

29.5.1.1.1. Jacobian of the one-step transformation#

Thus, the Jacobian reads

\[\begin{split}\mathbf{J} = \partial_{\mathbf{z}_{n}} \mathbf{z}_{n+1} = \begin{bmatrix} \mathbf{I} + \frac{h^2}{2} \partial_{\widetilde{\mathbf{x}}} \mathbf{f} & h + \frac{h^3}{4} \partial_{\widetilde{\mathbf{x}}} \mathbf{f} \\ h \partial_{\widetilde{\mathbf{x}}} \mathbf{f} & \mathbf{I} + \frac{h^2}{2} \partial_{\widetilde{\mathbf{x}}} \mathbf{f} \end{bmatrix}\end{split}\]

and the value of its determinant can be evaluated - for \(h\) “small enough” the upper-left block is approximately the identity and thus invertible; all the blocks commute, as they’re linear combinations of zero and first power of \(\partial_{\widetilde{\mathbf{x}}} \mathbf{f}\) - as

\[ | \mathbf{J} | = \left| \left( \mathbf{I} + \frac{h^2}{2} \partial_{\widetilde{\mathbf{x}}} \mathbf{f} \right)^2 -h^2 \partial_{\widetilde{\mathbf{x}}} \mathbf{f} \left( \mathbf{I} + \frac{h^2}{4} \partial_{\widetilde{\mathbf{x}}} \mathbf{f} \right) \right| = | \mathbf{I} | \ .\]
Determinant of block matrices - with Schur complement

The determinant of the block matrices - with Schur complement, if it exists -

\[\begin{split}\begin{bmatrix} A & B \\ C & D \end{bmatrix} = \begin{bmatrix} I & 0 \\ C A^{-1} & I \end{bmatrix} \begin{bmatrix} A & 0 \\ 0 & D - C A^{-1} B \end{bmatrix}\begin{bmatrix} I & A^{-1} B \\ 0 & I \end{bmatrix}\end{split}\]

reads - just take the determinants and exploit the property of determinant of the product is equal to the product of the determinants, and the determinant of a block-diagonal matrix -

\[\begin{split}\left|\begin{bmatrix} A & B \\ C & D \end{bmatrix} \right| = |A||D - C A^{-1} B| \ .\end{split}\]

If \(A\), \(B\), \(C\), \(D\) have proper dimensions and \(AC = CA\), i.e. commute, it follows

\[\begin{split}\begin{aligned} \left|\begin{bmatrix} A & B \\ C & D \end{bmatrix} \right| & = |A||D - C A^{-1} B| = \\ & = |A D - A C A^{-1} B | = \\ & = |A D - C A A^{-1} B | = \\ & = |A D - C B | \ . \end{aligned}\end{split}\]
Proof of inverse of block matrix by direct computation
\[\begin{split}\begin{aligned} \begin{bmatrix} I & 0 \\ C A^{-1} & I \end{bmatrix} \begin{bmatrix} A & 0 \\ 0 & D - C A^{-1} B \end{bmatrix}\begin{bmatrix} I & A^{-1} B \\ 0 & I \end{bmatrix} & = \begin{bmatrix} A & 0 \\ C & D - C A^{-1} B \end{bmatrix}\begin{bmatrix} I & A^{-1} B \\ 0 & I \end{bmatrix} = \\ & = \begin{bmatrix} A & B \\ C & D \end{bmatrix} \end{aligned}\end{split}\]

29.5.1.1.2. One-step energy difference#

29.5.1.2. Explicit Euler#

For the mechanical system governed by the differential equation \(\ddot{\mathbf{x}} = \mathbf{f}(\mathbf{x},t)\), or by the first-order system

\[\begin{split}\dfrac{d}{dt}\begin{bmatrix} \mathbf{x} \\ \mathbf{v} \end{bmatrix} = \begin{bmatrix} \mathbf{v} \\ \mathbf{f}(\mathbf{x}) \end{bmatrix}\end{split}\]

the numerical scheme, and the transformation \(\mathbf{z}_{n} \mapsto \mathbf{z}_{n+1}\), reads

\[\begin{split}\begin{aligned} \mathbf{x}_{n+1} & = \mathbf{x}_{n} + h \mathbf{v}_{n} \\ \mathbf{v}_{n+1} & = \mathbf{v}_n + h \mathbf{f}( \mathbf{x}_{n})\\ \end{aligned}\end{split}\]

29.5.1.2.1. Jacobian of the one-step transformation#

Thus, the Jacobian reads

\[\begin{split}\mathbf{J} = \partial_{\mathbf{z}_{n}} \mathbf{z}_{n+1} = \begin{bmatrix} \mathbf{I} & h \mathbf{I} \\ h \partial_{\widetilde{\mathbf{x}}} \mathbf{f} & \mathbf{I} \end{bmatrix}\end{split}\]

and the value of its determinant can be evaluated - for \(h\) “small enough” the upper-left block is approximately the identity and thus invertible; all the blocks commute, as they’re linear combinations of zero and first power of \(\partial_{\widetilde{\mathbf{x}}} \mathbf{f}\) - as

\[ | \mathbf{J} | = \left| \mathbf{I} - h^2 \partial_{\widetilde{\mathbf{x}}} \mathbf{f} \right| \ .\]

First-derivative of the force w.r.t. \(\mathbf{x}\) can be related to the stiffness matrix around equilibria,

\[\mathbf{f}(\overline{\mathbf{x}} + \delta \mathbf{x} ) \simeq \mathbf{f}(\overline{\mathbf{x}}) + \delta \mathbf{x} \cdot \partial_{\mathbf{x}} \mathbf{f}(\overline{\mathbf{x}}) = \mathbf{K} \delta \mathbf{x} \ .\]

For marginally-stable systems, the stiffness matrix has imaginary eigenvalues, \(\lambda_k(\mathbf{K}) = \mp i \Omega_k\). For the numerical integration of system around an equilibrium, it follows

\[ | \mathbf{J} | = \left| \mathbf{I} - h^2 \mathbf{K} \right| \ , \]

and thus - recalling that the determinant is equal to the product of the eigenvalues -

\[|\mathbf{J}| = \prod_k \lambda_k( \mathbf{I} - h^2 \mathbf{K} ) = \prod_k \left( 1 - h^2 \lambda_k( \mathbf{K} ) \right) = \prod_k ( 1 + h^2 \Omega_k^2 ) \ge 1 \ .\]
Eigenproblem of \(\ \mathbf{A} + b \mathbf{I}\)

If \(\mathbf{A} \mathbf{v}_k = \alpha_k \mathbf{v}_k\), then the matrix \(\mathbf{A} + b \mathbf{I}\) has the same eigenvectors \(\mathbf{v}_k\) with eigenvalues \(\lambda_k = \alpha_k + b\)

29.5.1.2.2. One-step energy difference#

29.5.1.3. Implicit Euler#

For the mechanical system governed by the differential equation \(\ddot{\mathbf{x}} = \mathbf{f}(\mathbf{x},t)\), or by the first-order system

\[\begin{split}\dfrac{d}{dt}\begin{bmatrix} \mathbf{x} \\ \mathbf{v} \end{bmatrix} = \begin{bmatrix} \mathbf{v} \\ \mathbf{f}(\mathbf{x}) \end{bmatrix}\end{split}\]

the numerical scheme, and the transformation \(\mathbf{z}_{n} \mapsto \mathbf{z}_{n+1}\), reads

\[\begin{split}\begin{aligned} \mathbf{x}_{n+1} & = \mathbf{x}_{n} + h \mathbf{v}_{n+1} \\ \mathbf{v}_{n+1} & = \mathbf{v}_n + h \mathbf{f}( \mathbf{x}_{n+1})\\ \end{aligned}\end{split}\]

29.5.1.3.1. Jacobian of the one-step transformation#

The Jacobian of the inverse transformation \(\mathbf{z}_{n+1} \mapsto \mathbf{z}_n\) reads

\[\begin{split}\mathbf{J}^{-1} = \partial_{\mathbf{z}_{n+1}} \mathbf{z}_{n} = \begin{bmatrix} \mathbf{I} & -h \mathbf{I} \\ -h \partial_{\widetilde{\mathbf{x}}} \mathbf{f} & \mathbf{I} \end{bmatrix}\end{split}\]

and the value of its determinant can be evaluated - for \(h\) “small enough” the upper-left block is approximately the identity and thus invertible; all the blocks commute, as they’re linear combinations of zero and first power of \(\partial_{\widetilde{\mathbf{x}}} \mathbf{f}\) - as

\[ | \mathbf{J}^{-1} | = \left| \mathbf{I} - h^2 \partial_{\widetilde{\mathbf{x}}} \mathbf{f} \right| \ .\]

First-derivative of the force w.r.t. \(\mathbf{x}\) can be related to the stiffness matrix around equilibria,

\[\mathbf{f}(\overline{\mathbf{x}} + \delta \mathbf{x} ) \simeq \mathbf{f}(\overline{\mathbf{x}}) + \delta \mathbf{x} \cdot \partial_{\mathbf{x}} \mathbf{f}(\overline{\mathbf{x}}) = \mathbf{K} \delta \mathbf{x} \ .\]

For marginally-stable systems, the stiffness matrix has imaginary eigenvalues, \(\lambda_k(\mathbf{K}) = \mp i \Omega_k\). For the numerical integration of system around an equilibrium, it follows

\[ | \mathbf{J}^{-1} | = \left| \mathbf{I} - h^2 \mathbf{K} \right| \ , \]

and thus - recalling that the determinant is equal to the product of the eigenvalues -

\[|\mathbf{J}^{-1}| = \prod_k \lambda_k( \mathbf{I} - h^2 \mathbf{K} ) = \prod_k \left( 1 - h^2 \lambda_k( \mathbf{K} ) \right) = \prod_k ( 1 + h^2 \Omega_k^2 ) \ge 1 \ ,\]

and thus

\[\left| \mathbf{J} \right| = \left| \mathbf{J}^{-1} \right|^{-1} = \left[ \prod_k ( 1 + h^2 \Omega_k^2 ) \right]^{-1} \le 1 \ .\]

29.5.1.3.2. One-step energy difference#

29.5.1.4. Crank-Nicolson#