<EMAIL_ADDRESS><EMAIL_ADDRESS>
# Application of Explicit Symplectic Algorithms to Integration of Damping
Oscillators

Tianshu Luo  Yimu Guo  Institute of Applied Mechanics, Department of
Mechanics, Zhejiang University, Hangzhou, Zhejiang, 310027, P.R.China
Institute of Applied Mechanics, Department of Mechanics, Zhejiang University,
Hangzhou, Zhejiang, 310027, P.R.China  [ [

(Received: date / Accepted: date)

###### Abstract

In this paper an approach is outlined. With this approach some explicit
algorithms can be applied to solve the initial value problem of
$n-$dimensional damped oscillators. This approach is based upon following
structure: for any non-conservative classical mechanical system and arbitrary
initial conditions, there exists a conservative system; both systems share one
and only one common phase curve; and, the value of the Hamiltonian of the
conservative system is, up to an additive constant, equal to the total energy
of the non-conservative system on the aforementioned phase curve, the constant
depending on the initial conditions. A key way applying explicit symplectic
algorithms to damping oscillators is that by the Newton-Laplace principle the
nonconservative force can be reasonably assumed to be equal to a function of a
component of generalized coordinates $q_{i}$ along a phase curve, such that
the damping force can be represented as a function analogous to an elastic
restoring force numerically in advance. Two numerical examples are given to
demonstrate the good characteristics of the algorithms.

###### keywords:

Hamiltonian, dissipation, non-conservative system, damping, explicit
symplectic algorithm

††journal: Communications in Nonlinear Science and Numerical Simulation

##  1 Introduction

Feng[1, 2, 3, 4],Marsden[5],Neri[6] and Yoshida[7]had developed a series of
symplectic algorithms for Hamiltonian systems. These algorithms possess some
advantages. But it is difficult to apply these algorithms to damping dynamical
systems, because it has been stated in most classical textbooks that the
Hamiltonian formalism focuses on solving conservative problems. Damping
phenomena is very important in the modeling of dynamical systems, and can not
be avoided. Our aim is to apply some explicit canonical algorithms to
nonlinear damping dynamical systems, which is generated generally by FE-
method. These canonical algorithms reported in this paper can be readily
utilized for computing large-scale nonlinear damped dynamical systems.

Betch[8][9][10] attempted to apply directly some implicit algorithms to
damping systems. The implicit symplectic algorithms utilized by Betch[8]
possess a few good characteristic, e.g. energy-conservation, momentum-
consistence, etc… In terms of energy-conservation, implicit symplectic
algorithms might be better than explicit symplectic ones. But explicit
symplectic schemes might be more suitable for nonlinear problems.

If one needs to apply symplectic algorithms to a dissipative system, one must
convert the dissipative system into a Hamiltonian system or find some
relationship between the dissipative system and a conservative one.

In the literature[11], we have stated a proposition describing a relation
among a damping dynamical system and conservative ones:

###### Proposition 1.1

For any non-conservative classical mechanical system and arbitrary initial
condition, there exists a conservative system; both systems sharing one and
only one common phase curve; and the value of the Hamiltonian of the
conservative system is equal to the sum of the total energy of the non-
conservative system on the aforementioned phase curve and a constant depending
on the initial condition.

In other words, a dissipative ordinary equation and a conservative equation
may possess a common particular solution. In the next section, an analytical
examples are given to explain this proposition. Readers can find the detailed
proof of Proposition 1.1 in the reference[11]

In the Literature [12] a basic explicit canonical integrator is proposed.
Based on this basic scheme, Neri[6] constructed 4-order explicit canonical
integrator, and then Yoshida [7] proposed a general method to construct higher
order explicit symplectic integrator. Utilizing the Proposition 1.1, we apply
this class of explicit canonical integrators to damping dynamical systems.
This point will be in detail stated in sec. 3.

##  2 One-dimensional Analytical Example

Consider a special one-dimensional simple mechanical system:

$\ddot{x}+c\dot{x}=0,$ (1)

where $c$ is a constant. The exact solution of the equation above is

$x=A_{1}+A_{2}e^{-ct},$ (2)

where $A_{1},A_{2}$ are constants. Differentiation gives the velocity:

$\dot{x}=-cA_{2}e^{-ct}.$ (3)

From the initial condition $x_{0},\dot{x}_{0}$, we find
$A_{1}=x_{0}+\dot{x}_{0}/c,A_{2}=-\dot{x}_{0}/c$. Inverting Eq. (2) yields

$t=-\frac{1}{c}\ln\frac{x-A_{1}}{A_{2}}$ (4)

and by substituting into Eq. (3), such we have

$\dot{x}=-c(x-A_{1})$ (5)

The dissipative force $F$ in the dissipative system (1) is

$F=c\dot{x}.$ (6)

Substituting Eq. (5) into Eq. (6), the conservative force $\mathcal{F}$ is
expressed as

$\mathcal{F}=-c^{2}(x-A_{1});$ (7)

Clearly, the conservative force $\mathcal{F}$ depends on the initial condition
of the dissipative system (1), in other words, an initial condition determines
a conservative force. Consequently, a new conservative system yields

$\ddot{x}+\mathcal{F}=0\rightarrow\ddot{x}-c^{2}(x-A_{1})=0.$ (8)

The stiffness coefficient in this equation must be negative. One can readily
verify that the particular solution (2) of the dissipative system can satisfy
the conservative one (8). This point agrees with Proposition (1.1).

The potential of the conservative system (8) is

$V=\int_{0}^{x}\left[-c^{2}(x-A_{1})\right]\mathrm{d}x=-\frac{c^{2}}{2}x^{2}+c^{2}A_{1}x$

Therefore the Hamiltonian is

$\hat{H}=T+V=\frac{1}{2}p^{2}-\frac{c^{2}}{2}x^{2}+c^{2}A_{1}x,$

where $p=\dot{x}$.

Furthermore, Proposition (1.1) can be depicted by Fig. 1. The phase flow of
conservative system (2) transforms the red area in phase space to the purple
area; the phase flow of conservative system (8) transforms the red area to the
green area. The blue curve in Fig. 1 illustrates the common phase curve. If
one draws more common phase curves and phase flows, the picture will like a
flower, the phase flow of the nonconservative system likes a pistil and phase
flows conservative systems like petals.

Figure 1: Relationship between nonconservative system (1) and conservative one
(8)

##  3 Modification Symplectic Numerical Schemes

###  3.1 Basic Explicit Symplectic Numerical Schemes

In the paper[12][6][7] a symplectic algorithm based second kind generation
function was stated：

$\begin{array}[]{l}{{\bm{p}}^{i+1}}={{\bm{p}}^{i}}-\tau{H_{q}}({{\bm{p}}^{i+1}},{{\bm{q}}^{i}})\\\
{{\bm{q}}^{i+1}}={{\bm{q}}^{i}}+\tau{H_{p}}({{\bm{p}}^{i+1}},{{\bm{q}}^{i}}),\\\
\end{array}$ (9)

where the superscript $i$ denotes the $i$-th time node, $\bm{q}$ denotes
coordinates and $\bm{p}$ denotes canonical momenta, and $H$ denotes
Hamiltonian quantity, $H_{q}=\partial H/\partial\bm{q},\ \ H_{p}=\partial
H/\partial\bm{p}$. If the Hamiltonian is seperable, i.e.
$H=U(\bm{p})+V(\bm{q}),V_{q}=H_{q},U_{p}=H_{p}$, then the symplectic scheme(9)
above becomes an explicit symplectic scheme:

$\begin{array}[]{l}{{\bm{p}}^{i+1}}={{\bm{p}}^{i}}-\tau{V_{q}}({{\bm{q}}^{i}})\\\
{{\bm{q}}^{i+1}}={{\bm{q}}^{i}}+\tau{U_{p}}({{\bm{p}}^{i+1}}).\\\ \end{array}$
(10)

For some nonlinear vibration mechanical system,
$V_{q}=\mathsfsl{K}(\bm{q})\bm{q}$.

Let us consider an $n-$dimensional nonlinear oscillator:

$\ddot{\bm{q}}+\mathsfsl{C}\dot{\bm{q}}+\mathsfsl{K}\bm{q}=0,$ (11)

where $\mathsfsl{C}$ denotes a non-linear damping coefficient matrix which
depends on $\bm{q}$, and $\mathsfsl{K}$ denotes a non-linear stiffness matrix
which depends on $\bm{q}$ and consists of two parts
$\mathsfsl{K}=\mathsfsl{\check{K}}+\mathsfsl{\hat{K}}$($\mathsfsl{\check{K}}$
is a diagonal matrix).

In accordance with Proposition 1.1, a conservative mechanical system was found
associated with the dissipative system (11) in addition to its initial
conditions. Subject to these initial conditions, the dissipative system (11)
possesses a common phase curve $\gamma$ with the conservative system. As in
Eq. (7), we can consider that the components of the damping force
$\mathsfsl{C}\dot{\bm{q}}$ determine the components of a conservative force on
the phase curve $\gamma$

$\begin{array}[]{ccc}c_{11}\dot{q}_{1}=\varrho_{11}(q_{1})&\dots&c_{1n}\dot{q}_{n}=\varrho_{1n}(q_{1})\\\
\vdots&\ddots&\vdots\\\
c_{n1}\dot{q}_{1}=\varrho_{21}(q_{n})&\dots&c_{nn}\dot{q}_{n}=\varrho_{nn}(q_{n}).\end{array}$
(12)

For convenience, this conservative force is assumed to be an elastic restoring
force:

$\begin{array}[]{ccc}\varrho_{11}(q_{1})=\kappa_{11}(q_{1})q_{1}&\dots&\varrho_{1n}(q_{1})=\kappa_{1n}(q_{1})q_{1}\\\
\vdots&\ddots&\vdots\\\
\varrho_{n1}(q_{1})=\kappa_{n1}(q_{n})q_{n}&\dots&\varrho_{nn}(q_{n})=\kappa_{nn}(q_{n})q_{n}.\end{array}$
(13)

In a similar manner, the components of the non-conservative force
$\mathsfsl{\hat{K}}\bm{q}$ are equal to the components of a conservative force
on the phase curve $\gamma$

$\begin{array}[]{ccc}\hat{K}_{11}q_{1}=\chi_{11}(q_{1})&\dots&\hat{K}_{1n}q_{n}=\chi_{1n}(q_{1})\\\
\vdots&\ddots&\vdots\\\
\hat{K}_{n1}q_{1}=\chi_{21}(q_{n})&\dots&\hat{K}_{nn}q_{n}=\chi_{nn}(q_{n}).\end{array}$
(14)

The conservative force can likewise be assumed to an elastic restoring force:

$\begin{array}[]{ccc}\chi_{11}(q_{1})=\lambda_{11}(q_{1})q_{1}&\dots&\chi_{1n}(q_{1})=\lambda_{1n}(q_{1})q_{1}\\\
\vdots&\ddots&\vdots\\\
\chi_{n1}(q_{1})=\lambda_{n1}(q_{n})q_{n}&\dots&\chi_{nn}(q_{n})=\lambda_{nn}(q_{n})q_{n}.\end{array}$
(15)

By an appropriate transformation, an equivalent stiffness matrix
$\mathsfsl{\tilde{K}}$ that is diagonal in form can be obtained

$\mathsfsl{\tilde{K}}_{ii}=\sum_{l=1}^{n}\kappa_{il}(q_{l})+\lambda_{il}(q_{l}).$
(16)

Consequently, an $n$-dimensional conservative system is obtained

$\bm{\ddot{q}}+(\mathsfsl{\check{K}}+\mathsfsl{\tilde{K}})\bm{q}=0$ (17)

which shares the common phase curve $\gamma$ with the $n$-dimensional damping
system described by (11). In this paper, the conservative system is called the
’substitute’ conservative system. The Lagrangian of Eqs.(17) is

$\hat{L}=\frac{1}{2}\dot{\bm{q}}^{T}\dot{\bm{q}}-\int_{\bm{0}}^{\bm{q}}(\mathsfsl{\check{K}}\bm{q})^{T}\mathrm{d}\bm{q}-\int_{\bm{0}}^{\bm{q}}(\tilde{\mathsfsl{K}}\bm{q})^{T}\mathrm{d}\bm{q},$
(18)

with the Hamiltonian

$\hat{H}=\frac{1}{2}\bm{p}^{T}\bm{p}+\int_{\bm{0}}^{\bm{q}}(\mathsfsl{\check{K}}\bm{q})^{T}\mathrm{d}\bm{q}+\int_{\bm{0}}^{\bm{q}}(\tilde{\mathsfsl{K}}\bm{q})^{T}\mathrm{d}\bm{q},$
(19)

where $\bm{0}$ is the zero vector, and $\bm{p}=\dot{\bm{q}}$. Here $\hat{H}$
in Eq. (19) is the mechanical energy of the conservative system (17), because
$\int_{\bm{0}}^{\bm{q}}(\tilde{\mathsfsl{K}}\bm{q})^{T}\mathrm{d}\bm{q}$ is a
potential function such that $\hat{H}$ is independent of the path taken in
phase space.

Subject to a certain initial condition, one need merely to solve the
conservative system(17). But one must in advance obtain the numerical
approximation of the matrix $\mathsfsl{\tilde{K}}$ for a time step, such that
one can utilize the algorithm (10) to integrate the conservative system (17)
for a time step. One can repeat this process above up to the end. In this way
one obtains the numerical particular solution of the conservative system (17),
which is exactly the numerical particular solution of the damping one. The he
numerical approximation of the matrix $\mathsfsl{\check{K}}$ can be assumed
as:

$\displaystyle\mathsfsl{\tilde{K}}=\left[\begin{array}[]{ccc}\tilde{K}_{11}&\dots&0\\\
\vdots&\ddots&\vdots\\\ 0&\dots&\tilde{K}_{nn}\end{array}\right]$ (23)
$\displaystyle{\tilde{K}_{j}}({q_{j}}^{i})={c_{jl}}\dot{q}_{l}^{i}/{q_{j}}^{i}+\hat{K}_{jl}q_{l}^{i}/{q_{j}}^{i}$

Hence the explicit canonical scheme (10) can be modified into

$\begin{array}[]{l}\tilde{K}_{j}^{i}({q_{j}}^{i})={c_{jl}}\dot{q}_{l}^{i}/{q_{j}}^{i}+\hat{K}_{jl}q_{l}^{i}/{q_{j}}^{i}\\\
{p_{j}}^{i+1}={p_{j}}^{i}-\tau[{K_{j}}+\tilde{K}_{j}^{i}({q_{j}}^{i})]{q^{i}})\\\
{q_{j}}^{i+1}={q^{i}}+\tau{p_{j}}^{i+1}\\\ \end{array}$ (24)

The scheme above is a one order scheme. Furthermore one can construct higher
order explicit canonical schemes utilizing the method reported in the
literatures[6][7]. Now consider a map from $\bm{z}=\bm{z}(0)$ to
$\bm{z}^{\prime}=\bm{z}(\tau)$:

${\bm{z^{\prime}}}\approx(\prod\limits_{i=1}^{h}{{e^{{r_{i}}t{\mathsfsl{E}}}}}{e^{{s_{i}}\tau{\mathsfsl{F}}}}+O({\tau^{n+1}})){\mathsfsl{z}},$
(25)

where

$\displaystyle\bm{z}=\left[\begin{array}[]{l}\bm{p},\\\
\bm{q}\end{array}\right],\bm{z}^{\prime}=\left[\begin{array}[]{l}\bm{p}^{\prime},\\\
\bm{q}^{\prime}\end{array}\right],$
$\displaystyle{\mathsfsl{E}}=\left[{\begin{array}[]{*{20}{c}}0&0\\\ 1&0\\\
\end{array}}\right]{\mathsfsl{F}}=\left[{\begin{array}[]{*{20}{c}}0&{-({\mathsfsl{K}}+{\mathsfsl{\tilde{K}}})}\\\
0&0\\\ \end{array}}\right].$

In fact Eq.(25) is the succession of the following mappings,

$\begin{array}[]{l}{{\mathsfsl{p}}^{j+1}}={{\mathsfsl{p}}^{j}}-{s_{i}}\tau{V_{q}}({{\mathsfsl{q}}^{j}})\\\
{{\mathsfsl{q}}^{j+1}}={{\mathsfsl{q}}^{j}}+{r_{i}}\tau{U_{p}}({{\mathsfsl{p}}^{j+1}})\\\
\end{array}.$ (28)

In reality the difference between the equations above and Eq.(24) is that the
coefficients $s_{i},r_{i}$ before the time step $\tau$. In the literature [7]
a generalized method to determine $s_{i},r_{i}$ were given. Therefore, the
higher order explicit canonical scheme can be represented as:

$\begin{array}[]{l}{\mathsfsl{\tilde{K}}}({q^{j}})=\left[{\begin{array}[]{*{20}{c}}{{\tilde{K}_{1}}(q_{1}^{j})}&{}\hfil&0\\\
{}\hfil&\ddots&{}\hfil\\\ 0&{}\hfil&{{\tilde{K}_{n}}(q_{n}^{j})}\\\
\end{array}}\right]{\tilde{K}_{\alpha}}(q_{\alpha}^{j})=\sum\limits_{l=1}^{n}{{c_{\alpha
l}}\dot{q}_{l}^{j}/q_{\alpha}^{j}}+\hat{K}_{\alpha
l}q_{l}^{i}/q_{\alpha}^{i}\\\
{\mathsfsl{E}}=\left[{\begin{array}[]{*{20}{c}}0&0\\\ 1&0\\\
\end{array}}\right]\;\;\;\;{\mathsfsl{F}}=\left[{\begin{array}[]{*{20}{c}}0&{-({\mathsfsl{K}}+{\mathsfsl{\tilde{K}}})}\\\
0&0\\\ \end{array}}\right]\\\
{{\bm{z}}^{j+1}}=(\prod\limits_{i=1}^{h}{{e^{{s_{i}}\tau{\mathsfsl{F}}}}{e^{{r_{i}}\tau{\mathsfsl{E}}}}}){{\bm{z}}^{j}}\\\
\end{array}\ $ (29)

##  4 Numerical Examples

Two examples will be given to shown this numerical method29.

###  4.1 The First Example

To begin, we consider a Van Der Pol’s oscillator

$\ddot{x}+\mu\dot{x}({x^{2}}-1)+x=0,$ (30)

where $\mu=10$. The initial conditions are given by ${x_{0}}=1,\;\;\dot{x}=0$.
We employee the $4-$order explicit symplectic method (29) with coefficients

$\displaystyle s_{1}=s_{4}=[2+(\sqrt[3]{2}+1/\sqrt[3]{2})]/6,\ \
s_{2}=s_{3}=[1-(\sqrt[3]{2}+1/\sqrt[3]{2})]/6,$ $\displaystyle\ \
r_{1}=r_{3}=[2+(\sqrt[3]{2}+1/\sqrt[3]{2})]/3,\ \
r_{2}=-[2+(\sqrt[3]{2}+1/\sqrt[3]{2})]/3,\ \ r_{4}=0,$

and classical explicit $4-$order Runge-Kutta method to compute the resonance
of the Van Der Pol’s oscillator (31) respectively, then employ a same method
to integrate the results to the total energy, which is the sum of the
mechanical energy and the work done by damping forces in the system (30). The
both methods are run with the same step size $\tau=0.01$. The resonance is
shown in Fig. 2, and the total energy is shown in Fig. 3.

Figure 2: The resonance of the Van Der Pol’s oscillator

Figure 3: The total energy of the Van Der Pol’s oscillator

It is aparent from Fig. 3 that the explicit symplectic method (29) has
qualitatively different behavior to the Runge-Kutta method. The energy
divergence between the explicit symplectic method and the exact solution is
smaller than that between Runge-Kutta method and the exact solution. The
energy divergence between the explicit symplectic method and Runge-Kutta
method increases with the time evolution. Due to the increasement of the
energy, the phase difference between both the results in Fig. 2 increases also
with the time evolution.

###  4.2 The Second Example

In the second example, we consider a $2-$dimensional damped nonlinear Duffing
oscillator

$\begin{array}[]{l}2\ddot{q}_{1}+0.1\dot{q}_{1}+(2+0.1q_{1}^{2})q_{1}+q_{2}=0\\\
3\ddot{q}_{2}+0.2\dot{q}_{2}+q_{1}+(2+0.2q_{2}^{2})q_{2}=0,\end{array}$ (31)

with the initial conditions $q_{1}=0,\ \ q_{2}=0,\ \ \dot{q}_{1}=0,\ \
\dot{q}_{2}=1$. The program of the both methods with the step size $\tau=0.01$
are carried out to simulate Eq. (31). The resonance is shown in Fig. 2, the
numerical solution of the total energy is shown in Fig.5.

There is only tiny difference between resonance results of the two methods,
correspondingly, the difference among the total energy obtained by the
numerical methods and anlytical methods is very tiny. As numerical examples in
the other literatures[13], that explicit Runge-Kutta method must cause
numerical pseudo dissipation which might be positive or negative. The
difference between our numerical examples and the examples in the
literature[13] is the total energy in our examples and the mechanical energy
in their examples111Fig.6.1 in the literature[13].

Figure 4: The $1$-th displacement of the damped Duffing oscillator

Figure 5: Total energy of the damped dissipative oscillator

##  5 Conclusions

We have introduced a class of explicit symplectic algorithms to dissipative
mechanical systems successfully, by changing these algorithms into the
scheme.(29). Because the algorithms (29) are explicit and possess good energy
preserving characteristics, the explicit symplectic algorithms (29) is quite
suitable for long term integration of arbitrary dimensional nonlinear
dissipative mechanical systems.

## References

  * Feng [1985] K. Feng,  On difference schemes and symplectic geometry,  in: Ed. Feng Kang Proceeding of the 1984 Beijing Symposium on differential geometry and differential equations-computation of partial differential equations, Science Press, Beijing, 1985, pp. 42–58.
  * Wu et al. [1989] H. Wu, M. Qin, K. Feng,  Construction of canonical difference schemes for hamiltonian formalism via generating functions,  JCM 7 (1989) 71–96.
  * Wu et al. [1990] H. Wu, M. Qin, K. Feng,  Symplectic difference schemes for the linear hamiltonian canonical systems,  JCM 8 (1990) 371–380.
  * Feng [1991] K. Feng,  The hamiltonian way for computing hamiltonian dynamics,  Math. Appl. 56 (1991) 17–35.
  * Marsden et al. [1998] J. E. Marsden, G. W. Patrick, S. Shkoller,  Multisymplectic geometry, variational integrators, and nonlinear pdes,  Communications in Mathematical Physics 199 (1998) 351–395. Cited By (since 1996): 129.
  * Neri [1988] F. Neri, Lie algebras and canonical integration, Technical Report, Department of Physics,University of Maryland, 1988.
  * Yoshida [1990] H. Yoshida,  Construction of higher order symplectic integrators,  Physics Letters A 150 (1990) 262–268.
  * Uhlar and Betsch [2010] S. Uhlar, P. Betsch,  On the derivation of energy consistent time stepping schemes for friction afflicted multibody systems,  Computers & Structures 88 (2010) 737 – 754.
  * Leyendecker et al. [2004] S. Leyendecker, P. Betsch, P. Steinmann,  Energy-conserving integration of constrained hamiltonian systems – a comparison of approaches,  Computational Mechanics 33 (2004) 174–185. 10.1007/s00466-003-0516-2.
  * Betsch [2006] P. Betsch,  Energy-consistent numerical integration of mechanical systems with mixed holonomic and nonholonomic constraints,  Computer Methods in Applied Mechanics and Engineering 195 (2006) 7020 – 7035. Multibody Dynamics Analysis.
  * Luo and Guo [2009] T. Luo, Y. Guo,  Infinite-dimensional Hamiltonian description of a class of dissipative mechanical systems,  ArXiv e-prints (2009).
  * Feng and Qin [1987] K. Feng, M. Qin,  The symplectic methods for the computation of hamiltonian equations,  in: Numerical Methods for Partial Differential Equations, Springer, Berlin, 1987, pp. 17–35.
  * Kane et al. [2000] C. Kane, J. E. Marsden, M. Ortiz, M. West,  Variational integrators and the newmark algorithm for conservative and dissipative mechanical systems,  International Journal for Numerical Methods in Engineering 49 (2000) 1295–1325.