#  Multiplicatively Repeated Non-Binary LDPC Codes

Kenta Kasai, David Declercq, Charly Poulliat, and Kohichi Sakaniwa, The
material in this paper was presented in part at 2010 IEEE International
Symposium on Information Theory (ISIT) [1]. K. Kasai and K. Sakaniwa are with
the Department of Communications and Integrated Systems, Graduate School of
Science and Engineering, Tokyo Institute of Technology, 2-12-1 Ookayama,
Meguro-ku, Tokyo 152-8550, Japan (e-mail: {kenta,
sakaniwa}@comm.ss.titech.ac.jp). D. Declercq and C. Poulliat are with the ETIS
Laboratory, UMR8051 ENSEA UCP CNRS, 95014 Cergy-Pontoise, France (e-mail:
{declercq, poulliat}@ensea.fr).

###### Abstract

We propose non-binary LDPC codes concatenated with multiplicative repetition
codes. By multiplicatively repeating the (2,3)-regular non-binary LDPC mother
code of rate 1/3, we construct rate-compatible codes of lower rates
$1/6,1/9,1/12,\dotsc$. Surprisingly, such simple low-rate non-binary LDPC
codes outperform the best low-rate binary LDPC codes so far. Moreover, we
propose the decoding algorithm for the proposed codes, which can be decoded
with almost the same computational complexity as that of the mother code.

###### Index Terms:

iterative decoding, low-rate code, non-binary low-density parity-check code,
rate compatible code, repetition code

##  I Introduction

In 1963, Gallager invented low-density parity-check (LDPC) codes [2]. Due to
sparsity of the code representation, LDPC codes are efficiently decoded by
belief propagation (BP) decoders. By a powerful optimization method density
evolution [3], developed by Richardson and Urbanke, messages of BP decoding
can be statistically evaluated. The optimized LDPC codes can approach very
close to Shannon limit [4].

Rate-adaptability is a desirable property of coding systems. Over time-varying
channels, the system adapts the coding rate according to the quality of the
channels. Using the different type of codes for different rates results in a
complex coding system. It is desirable to use a single encoder and decoder
pair compatible with different rates. Such a property of codes is referred to
as rate-compatibility. Moreover, rate-compatible codes allow us to transmit
bits gradually in conjunction with automatic repeat request (ARQ). By
puncturing a low rate code, we can construct rate-compatible codes of higher
rates.

In order to reliably transmit information over the very noisy communication
channels, one needs to encode the information at low coding rate. As described
in [5], one encounters a difficulty when designing low-rate LDPC codes. While,
for high rate codes, even binary regular LDPC codes have good thresholds. The
optimized low-rate structured LDPC codes, e.g. accumulate repeat accumulate
(ARA) code [6, Table. 1] of rate $1/6$ and multi-edge type LDPC code [5,
Table. X] of rate $1/10$ have good thresholds. However, the maximum row-
weights of those codes are as high as 11 and 28, respectively. Such high row
weights lead to dense parity-check matrices and degraded performance for short
code length. We note that, with very large code length, generalized LDPC codes
with Hadamard codes [7] perform very close to the ultimate Shannon limit [8].
However, the large code length leads to transmission latency. If two error
correcting codes with the same error-correcting capabilities and different
code length are given, the shorter code is preferred.

Another obstacle blocking the realization of the low-rate LDPC codes is the
large number of check node computations. For a fixed information length $K$,
it can be easily seen that the number $M$ of check nodes gets larger as the
coding rate $R$ gets lower. To be precise, $M=K(1-R)/R$. In the BP decoding,
computations of check nodes are usually more complex than those of variable
nodes. It is a desirable property for the low-rate LDPC codes to be decoded
with computational complexity comparable to that of the higher-rate LDPC
codes.

The problems for constructing low-rate LDPC codes are summarized as follows.

  * •

Problem 1: The Tanner graphs of low-rate LDPC codes tend to have many check
nodes that require more complex computations than variable nodes.

  * •

Problem 2: The Tanner graphs of optimized low-rate LDPC codes tend to have
check nodes of high degree, which results in the degraded decoding performance
for small code length.

  * •

Problem 3: The optimized low-rate LDPC codes need to be used with large code
length to exploit the potential decoding performance.

In this paper, we deal with all these issues.

In this paper, we consider non-binary LDPC codes defined by sparse parity-
check matrices over ${\mathrm{GF}}(2^{m})$ for $2^{m}>2$. Non-binary LDPC
codes were invented by Gallager [2]. Davey and MacKay [9] found non-binary
LDPC codes can outperform binary ones. Non-binary LDPC codes have captured
much attention recently due to their decoding performance [10, 11, 12, 13,
14].

It is known that irregularity of Tanner graphs help improve the decoding
performance of binary LDPC codes [4]. While, it is not the case for the non-
binary LDPC codes. The $(2,k)$-regular non-binary LDPC codes over
${\mathrm{GF}}(2^{m})$ are empirically known [15] as the best performing codes
for $2^{m}\geq 64$, especially for short code length. This means that, for
designing non-binary LDPC codes, one does not need to optimize the degree
distributions of Tanner graphs, since $(2,k)$-regular non-binary LDPC codes
are best. Furthermore, sparsity of $(2,k)$-regular Tanner graph helps
efficient decoding.

Sassatelli et al. proposed hybrid non-binary LDPC codes [16] whose symbols are
defined over the Galois fields of different sizes, e.g. over
${\mathrm{GF}}(2),{\mathrm{GF}}(8),$ and ${\mathrm{GF}}(16)$ and whose Tanner
graphs are irregular. In other words, the codes have two types of
irregularity, i.e. irregularity of the degree distributions of graphs and the
size distributions of Galois fields. To the best of the authors’ knowledge,
the decoding performance of the hybrid non-binary LDPC codes are best so far
among the low-rate codes of short code length.

In this paper, we investigate non-binary LDPC codes concatenated with
multiplicative repetition inner codes. We use a (2, 3)-regular non-binary LDPC
code of rate 1/3, as a mother code. By multiplicatively repeating the mother
code, we construct codes of lower rates 1/6,1/9,1/12,$\dotsc$. Furthermore, we
present a decoding algorithm for the proposed codes. And we show the
computational complexity of decoding is almost the same as that of the mother
code. The codes exhibit surprisingly better decoding performance than the best
codes so far for small and moderate code length.

The rest of this paper is organized as follows. Section II defines the
proposed codes. Section III describes the decoding algorithm for the proposed
codes. In Section IV, we investigate the thresholds for the proposed codes
transmitted over the binary erasure channels (BEC) by density evolution [4,
17]. In Section V, for the BEC and AWGN channels, we compare the decoding
performance of the proposed codes and the best known codes for short and
moderate code length.

##  II Concatenation of Non-Binary LDPC Codes and Multiplicative Repetition
Codes

We deal with elements of ${\mathrm{GF}}(2^{m})$ as non-binary symbols. For
transmitting over the binary input channels, each non-binary symbol in
${\mathrm{GF}}(2^{m})$ needs to be represented by a binary sequence of length
$m$. For each $m$, we fix a Galois field ${\mathrm{GF}}(2^{m})$ with a
primitive element $\alpha$ and its primitive polynomial $\pi$. Once a
primitive element $\alpha$ of ${\mathrm{GF}}(2^{m})$ is fixed, each symbol is
given a $m$-bit representation [18, pp. 110]. For example, with a primitive
element $\alpha\in{\mathrm{GF}}(2^{3})$ such that
$\pi(\alpha)=\alpha^{3}+\alpha+1=0$, each symbol is represented as
$0=(0,0,0)$, $1=(1,0,0)$, $\alpha=(0,1,0)$, $\alpha^{2}=(0,0,1)$,
$\alpha^{3}=(1,1,0)$, $\alpha^{4}=(0,1,1)$, $\alpha^{5}=(1,1,1)$ and
$\alpha^{6}=(1,0,1)$.

A non-binary LDPC code $C$ over ${\mathrm{GF}}(2^{m})$ is defined by the null
space of a sparse $M\times N$ parity-check matrix $H=\\{h_{ij}\\}$ defined
over ${\mathrm{GF}}(2^{m})$.

$\displaystyle C$ $\displaystyle=\\{x\in{\mathrm{GF}}(2^{m})^{N}\mid
Hx=0\in{\mathrm{GF}}(2^{m})^{M}\\}$

The $c$-th parity-check equation for $c=1,\dotsc,M$ is written as

$\displaystyle h_{c1}x_{1}+\cdots+h_{cN}x_{N}=0\in{\mathrm{GF}}(2^{m}),$

where $h_{c1},\dotsc,h_{cN}\in{\mathrm{GF}}(2^{m})$ and
$x_{1},\dotsc,x_{N}\in{\mathrm{GF}}(2^{m})$.

Figure 1: An example of a mother code $C_{1}$. A non-binary (2,3)-regular LDPC
code of rate 1/3 over ${\mathrm{GF}}(2^{m})$. Each variable node represents a
symbol in ${\mathrm{GF}}(2^{m})$. Each check node represents a parity-check
equation over ${\mathrm{GF}}(2^{m})$. The code length is 18 symbols in
${\mathrm{GF}}(2^{m})$ or equivalently $18m$ bits. Circle and square nodes
represent variable and check nodes, respectively. The lower-rate codes $C_{T}$
for $T=2,3\dotsc$ are constructed from $C_{1}$.

Binary LDPC codes are represented by Tanner graphs with variable and check
nodes [19, pp. 75]. The non-binary LDPC codes, in this paper, are also
represented by bipartite graphs with variable nodes and check nodes, which are
also referred to as Tanner graphs. For a given sparse parity-check matrix
$H=\\{h_{cv}\\}$ over ${\mathrm{GF}}(2^{m})$, the graph is defined as follows.
The $v$-th variable node and $c$-th check node are connected if $h_{cv}\neq
0$. By $v=1,\dotsc,N$ and $c=1,\dotsc,M$, we also denote the $v$-th variable
node and $c$-th check node, respectively.

A non-binary LDPC code with a parity-check matrix over ${\mathrm{GF}}(2^{m})$
is called $(d_{v},d_{c})$-regular if all the columns and all the rows of the
parity-check matrix have weight $d_{v}$ and $d_{c}$, respectively, or
equivalently all the variable and check nodes have degree $d_{v}$ and $d_{c}$,
respectively.

Let $C_{1}$ be a $(2,3)$-regular LDPC code defined over ${\mathrm{GF}}(2^{m})$
of length $N$ symbols or equivalently $mN$ bits and of rate $1/3$. The code
$C_{1}$ has a $2N/3\times N$ sparse parity-check matrix $H$ over
${\mathrm{GF}}(2^{m})$. The matrix $H$ has row weight 3 and column weight 2.
Fig. 1 shows the Tanner graph of an example $C_{1}$ of length $N=$18 symbols.

Figure 2: An example of $C_{2}$. A non-binary (2,3)-regular LDPC code over
${\mathrm{GF}}(2^{m})$ concatenated with inner multiplicative repetition codes
of length 2. The code length is 36 symbols or equivalently $36m$ bits. The
rate is 1/6.

By using $C_{1}$ as a mother code, we will construct codes
$C_{2},C_{3},\dotsc,C_{T}$ of lower rates in the following way. Choose $N$
coefficients $r_{N+1},\dotsc,r_{2N}$ uniformly at random from
${\mathrm{GF}}(2^{m})\setminus\\{0\\}$. The lower-rate code $C_{2}$ is
constructed as follows.

$\displaystyle C_{2}=\\{$
$\displaystyle(x_{1},\dotsc,x_{2N})|x_{N+v}=r_{N+v}x_{v},$
$\displaystyle\quad\text{ for }v=1,\dotsc,N,(x_{1},\dotsc,x_{N})\in C_{1}\\}.$

Since the resulting code $C_{2}$ has code length $2N$ and the same number of
codewords as $C_{1}$, then the rate is $1/6$. Fig. 2 shows the Tanner graph of
$C_{2}$ of length 2$N=$36 symbols. We say that $x_{N+v}=r_{N+v}x_{v}$ is a
multiplicative repetition symbol of $x_{v}$ for $v=1,\dotsc,N$. Each variable
node of degree one in Fig. 2 represents a multiplicative repetition symbol
$x_{N+v}$ for $v=1,\dotsc,N$. And each check node of degree two in Fig. 2
represents a parity-check constraint $x_{N+v}+r_{N+v}x_{v}=0$ for
$v=1,\dotsc,N$.

For $T\geq 3$, in a recursive fashion, by choosing $N$ coefficients
$r_{(T-1)N+1},\dotsc,r_{TN}$ randomly chosen from
${\mathrm{GF}}(2^{m})\setminus\\{0\\}$, the further low-rate code $C_{T}$ is
constructed from $C_{T-1}$ as follows.

$\displaystyle C_{T}=\\{$
$\displaystyle(x_{1},\dotsc,x_{TN})|x_{(T-1)N+v}=r_{(T-1)N+v}x_{v},$
$\displaystyle\text{ for }v=1,\dotsc,N,(x_{1},\dotsc,x_{(T-1)N})\in
C_{T-1}\\}.$

The code $C_{T}$ has length $TN$ and rate $1/(3T)$. Fig. 3 shows the Tanner
graph of $C_{3}$ of 3$N=$54 symbol code length. Fig. 4 shows the block diagram
of the encoding of $C_{3}.$ We refer to $T$ as the repetition parameter.

Concatenating a binary code with repetition codes is known as the worst coding
scheme. Indeed, repeating a binary code just doubles the number of channel use
without any improvement of the curve of the decoding error rate v.s.
$\mathrm{E_{b}/N_{0}}$. Note that the proposed code $C_{T}$ are not generated
by simple repetitions of the mother code but the random multiplicative
repetitions of non-binary symbols. Since the coefficient
$r_{N+v},\dotsc,r_{(T-1)N+v},\in{\mathrm{GF}}(2^{m})\setminus\\{0\\}$ is
randomly chosen, the multiplicative repetition
$x_{v}\mapsto(x_{v},r_{N+v}x_{v},\dotsc,r_{(T-1)N+v}x_{v},)$ can be viewed as
a random code of length $mT$ bits. In other words, the proposed codes can be
viewed as non-binary LDPC codes over ${\mathrm{GF}}(2^{m})$ serially
concatenated with $N$ random binary codes of length $mT$. Intuitively, this
explains why multiplicative repetition works better than simple repetition.

The construction of the proposed codes may remind some readers of Justesen
codes [20]. Note that the proposed construction chooses the multiplicative
coefficients uniformly at random. Note also that since the minimum distance of
$C_{1}$ is at most $O(\log(N))$ [15], the $C_{T}$ code has minimum distance is
at most $O(T\log(N))$.

Due to the repetition of symbols, the encoder are inherently rate-compatible.

##  III Decoding Scheme

The BP decoder for non-binary LDPC codes [21] exchanges probability vectors of
length $2^{m}$, called messages, between variable nodes and check nodes, at
each iteration round $\ell\geq 0$. The proposed codes $C_{T}$ for $T\geq 2$
also can be decoded by the BP decoding algorithm on the Tanner graphs of
$C_{T}$. In this section, instead of the immediate use of the BP decoding on
the Tanner graph of $C_{T}$, we propose a decoding algorithm which uses only
the Tanner graph of $C_{1}$ for decoding $C_{T}$ for $T\geq 2$.

Figure 3: An example of $C_{3}$. A non-binary (2,3)-regular LDPC code over
${\mathrm{GF}}(2^{m})$ concatenated with one inner multiplicative repetition
codes of length 3. The code length is 54 symbols or equivalently $54m$ bits.
The rate is 1/9.

The variable nodes of degree one in Fig. 2 and Fig. 3 represent multiplicative
repetition symbols of $C_{2}$ and $C_{3}$, respectively. If the BP decoding
algorithm is immediately applied to the proposed codes, all the variable nodes
and check nodes, including the variable nodes of those multiplicative
repetition symbols, are activated, i.e. exchage the messages. However, the
messages reaching the variable nodes of degree one do not change messages that
are sent back from the nodes. Therefore, the decoder does not need to pass the
messages all the way to those variable nodes of degree 1 and their adjacent
check nodes of degree 2. Consequently, after the variable nodes of degree 1
pass the initial messages to the upper part of the graph, the decoder uses
only the upper part of the graph, i.e. $C_{1}$.

SOURCE$C_{1}$$(x_{1},\dotsc,x_{N})$$(r_{N+1}x_{1},\dotsc,r_{2N}x_{N})$$(r_{2N+1}x_{1},\dotsc,r_{3N}x_{N})$$\begin{subarray}{c}\\!\\!\\!\\!\\!\\!\\!\\!{r_{N+v}}\\\
(v=1,\dotsc,N)\end{subarray}$$\begin{subarray}{c}\\!\\!\\!\\!\\!\\!\\!\\!{r_{2N+v}}\\\
(v=1,\dotsc,N)\end{subarray}$ Figure 4:  The block diagram of the encoder of
$C_{3}$. First, source of $N/3$ symbols in ${\mathrm{GF}}(2^{m})$ are encoded
with a (2,3)-regular LDPC code $C_{1}$ over ${\mathrm{GF}}(2^{m})$. Next, each
symbol in the codeword $x_{v}$, for $v=1,\dotsc,N$, is randomly multiplied by
$r_{N+v}$ and $r_{2N+v}$ from ${\mathrm{GF}}(2^{m})\setminus\\{0\\}$ to
generate $x_{N+v}$ and $x_{2N+v}$.

The computations of check nodes are more complex than those of variable nodes.
As posed in the Problem 1 in Section I, the number $M$ of the check nodes gets
higher as $R$ decreases. In general, LDPC codes of information length $K$ and
rate $R$ have $K(1-R)/R$ check nodes. In our setting, we have $K=N/3$
information symbols. The number of check nodes in the proposed code $C_{T}$
for $T\geq 2$ is also given by $K(1-R)/R$. However, $(T-1)N$ check nodes of
degree 2 adjacent to the $(T-1)N$ variable nodes of degree 1 do not need to
participate in the BP decoding iterations. The only $2N/3$ active check nodes
in the mother code $C_{1}$ participate in the BP decoding algorithm for
decoding $C_{T}$ for $T\geq 2$. Note that the number of active check nodes
$2N/3$ remains unchanged for any $T\geq 1$. This is highly preferable property
for low-rate LDPC codes, which relieves the Problem 1. The Problem 2 is also
relieved, since the maximum degree of check nodes in the mother code $C_{1}$
is as small as 3.

The BP decoding involves mainly 4 parts, i.e. the initialization, the check to
variable computation, the variable to check computation, and the tentative
decision parts. For $v=1,\dotsc,NT$, let $X_{v}$ be the random variables with
realizations $x_{v}$. Let $Y_{v}$ be the random variables with realizations
$y_{v}$ which is received value from the channel $\Pr(Y_{v}|X_{v})$ and the
probability of transmitted symbol $\Pr(X_{v})$ is assumed to be uniform.

We assume the decoder knows the channel transition probability

$\displaystyle\Pr(X_{v}=x|Y_{v}=y_{v}),v=1,\dotsc,NT$ (1)

for $x\in{\mathrm{GF}}(2^{m})$. When the transmissions take place over the
memoryless binary-input output-symmetric channels, we can rewrite (1) as

$\displaystyle\Pr(X_{v}=x|Y_{v}=y_{v})=\prod_{i=1}^{m}\Pr(X_{v,i}=x_{i}|Y_{v,i}=y_{v,i}),$

where $(x_{1},\dotsc,x_{m})\in{\mathrm{GF}}(2)^{m}$ is the binary
representation of $x\in{\mathrm{GF}}(2^{m})$ and $X_{v,i}$ is the random
variable of the transmitted bit, and the corresponding channel output
$y_{v,i}$ and its random variable $Y_{v,i}$.

###  III-A Decoding Algorithm

initialization :
For each variable node $v$ in $C_{1}$ for $v=1,\dotsc,N$, compute
$p_{v}^{(0)}(x)$ as follows.

$\displaystyle\begin{split}p_{v}^{(0)}(x)&=\xi\Pr(X_{v}=x|Y_{v}=y_{v})\\\
&\prod_{t=1}^{T-1}{\Pr(X_{tN+v}=r_{tN+v}x|Y_{tN+v}=y_{tN+v})},\end{split}$ (2)

for $x\in{\mathrm{GF}}(2^{m})$, where $\xi$ is the normalization factor so
that $\sum_{x\in{\mathrm{GF}}(2^{m})}p^{(0)}_{v}(x)=1$. Each variable node
$v=1,\dotsc,N$ in $C_{1}$ sends the initial message
$p_{vc}^{(0)}=p_{v}^{(0)}\in{\mathbb{R}}^{2^{m}}$ to each adjacent check node
$c$. Set the iteration round as $\ell:=0$.

check to variable :
For each check node $c=1,\dotsc,M$ in $C_{1}$, let $\partial c$ be the set of
the adjacent variable nodes of $c$. It holds that $\\#\partial c=3$, since the
mother code $C_{1}$ is $(2,3)$-regular. Each $c$ has 3 incoming messages
$p_{vc}^{(\ell)}$ for $v\in\partial c$ from the 3 adjacent variable nodes. The
check node $c$ sends the following message
${p}^{(\ell+1)}_{cv}\in{\mathbb{R}}^{2^{m}}$ to each adjacent variable node
$v\in\partial c$.

$\displaystyle\tilde{p}^{(\ell)}_{vc}(x)={p}^{(\ell)}_{vc}(h_{cv}^{-1}x)\text{
for $x\in{\mathrm{GF}}(2^{m})$},$
$\displaystyle\tilde{p}^{(\ell+1)}_{cv}=\otimes_{v^{\prime}\in\partial
c\setminus{\\{v\\}}}\tilde{p}^{(\ell)}_{v^{\prime}c},$
$\displaystyle{p}^{(\ell+1)}_{cv}(x)=\tilde{p}^{(\ell+1)}_{cv}(h_{cv}x)\text{
for $x\in{\mathrm{GF}}(2^{m})$}.$

where $p_{1}\otimes p_{2}\in{\mathbb{R}}^{2^{m}}$ is a convolution of
$p_{1}\in{\mathbb{R}}^{2^{m}}$ and $p_{2}\in{\mathbb{R}}^{2^{m}}$. To be
precise,

$(p_{1}\otimes
p_{2})(x)=\sum_{\begin{subarray}{c}y,z\in{\mathrm{GF}}(2^{m})\\\
x=y+z\end{subarray}}{p_{1}(y)p_{2}(z)}\text{ for $x\in{\mathrm{GF}}(2^{m})$}.$

The convolution seems the most complex part of the decoding algorithm. Indeed,
the convolutions are efficiently calculated via FFT and IFFT [22], [17].
Increment the iteration round as $\ell:=\ell+1$.

variable to check :
Each variable node $v=1,\dotsc,N$ in $C_{1}$ has two adjacent check nodes
since the mother code $C_{1}$ is $(2,3)$-regular. Let $c$ and $c^{\prime}$ be
the two adjacent check nodes of $v$. The message
$p^{(\ell)}_{vc}\in{\mathbb{R}}^{2^{m}}$ sent from $v$ to $c$ is given by

$\displaystyle p^{(\ell)}_{vc}(x)=\xi
p_{v}^{(0)}(x)p^{(\ell)}_{c^{\prime}v}(x)\text{ for
$x\in{\mathrm{GF}}(2^{m})$},$

where $\xi$ is the normalization factor so that
$\sum_{x\in{\mathrm{GF}}(2^{m})}p^{(\ell)}_{vc}(x)=1$.

tentative decision :
For each $v=1,\dotsc,N$, the tentatively estimated $v$-th transmitted symbol
is given as

$\displaystyle\hat{x}_{v}^{(\ell)}=\mathop{\mathrm{argmax}}_{x\in{\mathrm{GF}}(2^{m})}p_{v}^{(0)}(x)p^{(\ell)}_{cv}(x)p^{(\ell)}_{c^{\prime}v}(x),$

where $c$ and $c^{\prime}$ are the two adjacent check nodes of $v$. If
$\underline{\hat{x}}^{(\ell)}:=(\hat{x}_{1}^{(\ell)},\dotsc,\hat{x}_{N}^{(\ell)})$
forms a codeword of $C_{1}$, in other words, $\underline{\hat{x}}^{(\ell)}$
satisfies every parity-check equations

$\displaystyle\sum_{v\in\partial
c}h_{cv}\hat{x}_{v}^{(\ell)}=0\in{\mathrm{GF}}(2^{m})$

for all $c=1,\dotsc,M$, the decoder outputs $\underline{\hat{x}}^{(\ell)}$ as
the estimated codeword. Otherwise repeat the latter 3 decoding steps. If the
iteration round $\ell$ reaches a pre-determined number, the decoder outputs
FAIL.

The decoder is inherently rate-compatible. Indeed, for decoding the different
$C_{T}$ of rate $1/(3T)$ for $T=1,2,\dotsc$, the decoder only needs the Tanner
graph of the mother code $C_{1}$.

##  IV Erasure Channel Analysis

In the binary case, we can predict the asymptotic decoding performance of LDPC
codes transmitted over the general memoryless binary-input output-symmetric
channels in the large code length limit by density evolution [4]. Density
evolution also can be used to analyze non-binary LDPC codes [23, 24]. However,
for large field size, it becomes computationally intensive and tractable only
for the BEC.

Rathi and Urbanke developed the density evolution which enables the prediction
of the decoding performance of the non-binary LDPC codes over the BEC in the
limit of large code length. For a given code ensemble, density evolution gives
the maximum channel erasure probability at which the decoding erasure
probability, averaged over all the LDPC codes in the ensemble goes to zero.
The maximum channel erasure probability given by the density evolution is
referred to as the threshold.

It is shown in [17] that for the transmissions over the BEC with non-binary
LDPC codes defined over ${\mathrm{GF}}(2^{m})$, the decoding results depend on
the binary representation, i.e. the primitive element. In other words, two
isomorphic fields do not, in general, yield the identical decoding results.
Rathi and Urbanke also observed that the difference of the threshold is of the
order of $10^{-4}$ for the different fields. The density evolution [17] is
developed for the non-binary LDPC code ensembles with parity-check matrices
defined over the general linear group $\mathrm{GL}({\mathrm{GF}}(2),m)$. In
this section, we will use the density evolution to evaluate the thresholds of
non-binary LDPC codes defined over ${\mathrm{GF}}(2^{m})$. This is a fair
approximation, since in [17], it is reported that the threshold for the code
ensemble with parity-check matrices defined over ${\mathrm{GF}}(2^{m})$ and
$\mathrm{GL}({\mathrm{GF}}(2),m)$ have almost the same thresholds within the
order of $10^{-4}$.

When the transmission takes place over the BEC and all-zero codeword is
assumed to be sent, the messages, described by probability vectors
$(p(x))_{x\in{\mathrm{GF}}(2^{m})}$ of length $2^{m}$ in general, can be
reduced to linear subspaces [17] of ${\mathrm{GF}}(2)^{m}$. To be precise, for
each message in the BP decoding algorithm, a subset of ${\mathrm{GF}}(2)^{m}$

$\displaystyle\\{\bm{x}\in{\mathrm{GF}}(2)^{m}\mid p(x)\neq 0\\},$

forms a linear subspace of ${\mathrm{GF}}(2)^{m}$, where $\bm{x}$ is the
binary representation of $x\in{\mathrm{GF}}(2^{m})$.

Define $P^{(\ell)}=(P^{(\ell)}_{0},\dotsc,P^{(\ell)}_{m})$ and
$Q^{(\ell)}=(Q^{(\ell)}_{0},\dotsc,Q^{(\ell)}_{m})$ as the probability vectors
of length $m+1$ such that $P_{i}^{(\ell)}$ (resp. $Q_{i}^{(\ell)}$) is the
probability that a message sent from variable (resp. check) nodes has
dimension $i$ at the $\ell$-th iteration round of the BP decoding algorithm.
The density evolution gives us the update equations of $P^{(\ell)}$ and
$Q^{(\ell)}$ for $\ell\geq 0$.

Figure 5: Thresholds $\epsilon^{\ast}$ of $C_{T}$ over ${\mathrm{GF}}(2^{m})$
for the BEC and $T$=1,2,3 and 5 from below. The rate is $1/(3T)$. The straight
lines show the Shannon limits $1-1/(3T)$.

Rathi and Urbanke[17] developed the density evolution for the BEC that tracks
probability mass functions of the dimension of the linear subspaces. For
$\ell\geq 0$, the density evolution tracks the probability vectors
$P^{(\ell)}$ and $Q^{(\ell)}$ which are referred to as densities. The initial
messages in (2) can be seen as the intersection of $T$ subspaces of the
messages received as the channel outputs. The density of the initial messages
is given by $P^{(0)}$ as follows,

$\displaystyle P^{(0)}=\overbrace{E\boxdot\cdots\boxdot E}^{\text{$T$
times}},$ $\displaystyle E:=(E_{0},\dotsc,E_{m}),$ $\displaystyle
E_{i}:=\binom{m}{i}\epsilon^{i}(1-\epsilon)^{m-i},$

where $\epsilon$ is the channel erasure probability of the BEC. The operator
$\boxdot$ is defined as follows.

$\displaystyle\left[P\boxdot Q\right]_{k}$
$\displaystyle=\sum_{i=k}^{m}\sum_{j=k}^{k+m-i}C_{\boxdot}(m,k,i,j)P_{i}Q_{j},$
$\displaystyle{C_{\boxdot}}(m,k,i,j)$
$\displaystyle:=2^{(i-k)(j-k)}\frac{\left[\begin{array}[]{c}i\\\
k\end{array}\right]\left[\begin{array}[]{c}m-i\\\
j-k\end{array}\right]}{\left[\begin{array}[]{c}m\\\ j\end{array}\right]},$

where $\displaystyle\left[\begin{array}[]{c}m\\\
k\end{array}\right]=\prod_{l=0}^{k-1}\frac{2^{m}-2^{l}}{2^{k}-2^{l}}$ is a
2-Gaussian binomial.

Since the mother code is (2,3)-regular, the update equations of density
evolution is given by

$\displaystyle Q^{(\ell+1)}=P^{(\ell)}\boxtimes P^{(\ell)},$ $\displaystyle
P^{(\ell+1)}=P^{(0)}\boxdot Q^{(\ell+1)},$

where the operator $\boxtimes$ is defined as follows.

$\displaystyle\left[P\boxtimes Q\right]_{k}$
$\displaystyle=\sum_{i=0}^{k}\sum_{j=k-i}^{k}C_{\boxtimes}(m,k,i,j)P_{i}Q_{j},$
$\displaystyle{C_{\boxtimes}}(m,k,i,j)$
$\displaystyle:=2^{(k-i)(k-j)}\frac{\left[\begin{array}[]{c}m-i\\\
m-k\end{array}\right]\left[\begin{array}[]{c}i\\\
k-j\end{array}\right]}{\left[\begin{array}[]{c}m\\\ m-j\end{array}\right]}.$

Since the messages of dimension 0 corresponds to the successful decoding, the
threshold is defined as follows.

$\displaystyle\epsilon^{\ast}:=\sup_{\epsilon\in[0,1]}\\{\epsilon\in[0,1]\mid\lim_{\ell\to\infty}P^{(\ell)}_{0}=1\\}.$

In the large code length limit, if $\epsilon<\epsilon^{\ast}$ the reliable
transmissions are possible with the proposed $C_{T}$.

Fig. 5 draws the thresholds of $C_{T}$ defined with parity-check matrices over
$\mathrm{GL}_{m}({\mathrm{GF}}(2))$ for repetition parameter $T=1,\dotsc,5$
and $m=1,\dotsc,10$. The threshold $\epsilon^{\ast}=1/\sqrt[T]{2}$ for the
binary case $m=1$ is decided by the stability condition [19]. It can be seen
that the thresholds are not monotonic with respect to $m$. For repetition
parameter $T=1$, i.e., the mother code has the maximal threshold at $m=6$. For
$T>2$, the maximal threshold is attained around at $m=8$.

Fig. 6 compares the proposed codes and the best existing low-rate LDPC codes
respect to the thresholds for the BEC. It can be seen that the proposed codes
have better thresholds especially for lower rates.

Figure 6:  The asymptotic decoding performance over the BEC of the proposed
codes and the best known low-rate codes. One curve corresponds to the proposed
code $C_{T}$ over ${\mathrm{GF}}(2^{6})$ of coding rates $R=1/(3T)$ with
repetition parameter $T=1,\dotsc,10$. The punctured $C_{1}$ of rate 1/2 is
also plotted. The other curve corresponds to the bit-wise shortened non-binary
LDPC code over ${\mathrm{GF}}(2^{6})$ proposed by Klinc et al. [25, Fig. 1].
The vertical axis indicates $(1-\epsilon^{\ast}-R)/R$ which is the normalized
gap between the capacity $1-\epsilon^{\ast}$ and the rate $R$, where
$\epsilon^{\ast}$ is the threshold.  Figure 7:  The solid curve shows the
decoding performance of the proposed code $C_{2}$ whose mother code is a
$(2,4)$-regular non-binary LDPC code defined over ${\mathrm{GF}}(2^{8})$. The
coding rate is 1/4. The transmission takes place over the BEC. The arrow
indicates the threshold 0.72898 of $C_{2}$. The code length is of length 1024,
8192 and 65536. For comparison, the decoding performance of accumulated LDPC
(ALDPC) codes [26, Fig. 15] is shown. It is known that the ALDPC codes achieve
the capacity of the BEC in the limit of large code length and exhibit good
decoding performance with finite code length. The frame error rate of
corresponding random codes of rate 1/4 under maximum-likelihood decoding are
calculated by [27, Eq. (3.2)]. It can be seen that the proposed codes exhibit
better decoding performance than the ALDPC codes with code length up to 8192.
The error floors of the proposed codes can not be observed down to FER
$10^{-5}$ while the ALDPC codes have high error floors even with code length
as long as 65536 bits.  Figure 8:  The thresholds for the AWGN channels of the
proposed codes $C_{T}$ of rate $1/(3T)$ for $T=1,\dotsc,20$. The codes are
defined over ${\mathrm{GF}}(2^{8}).$ Figure 9:  Frame error rate versus
parameter for the proposed codes $C_{2}$ and $C_{7}$ over
${\mathrm{GF}}(2^{8})$ transmitting over the AWGN channel. The rates of
$C_{2}$ and $C_{7}$ are 1/6 and 1/21, respectively. The information length are
set to 1024, 4096, 16834, and 65536. The arrows indicate the corresponding
threshold values. Observe how the curves move closer to these threshold values
for increasing codeword lengths.  Figure 10: The frame error rate of the
proposed codes $C_{T}$ for $T=1,2,6$ and hybrid non-binary LDPC codes. It also
shows the performance of rate half punctured mother code $C_{1}$. All these
codes have 192 information bits. The curves labeled SP59 are the corresponding
Shannon’s 1959 sphere-packing bound [28, 29, 30] for rate 1/3, 1/6 and 1/18.
Figure 11: The frame error rate of the proposed codes $C_{T}$ for $T=1,2,6$,
multi-edge type (MET) LDPC code of rate 1/2 [5] and 1/6, and accumulate repeat
accumulate (ARA) code [6] of rate 1/6. All these codes have 1024 information
bits except that the MET LDPC code of rate 1/2 has 1280 information bits. It
also shows the performance of rate half punctured mother code $C_{1}$. The
curves labeled SP59 are the corresponding Shannon’s 1959 sphere-packing bounds
for rate 1/3, 1/6 and 1/18.

##  V Numerical Results

In this section, we give some numerical results of the proposed codes. Fig. 7
shows the decoding performance of the proposed code $C_{2}$ whose mother code
is a $(2,4)$-regular non-binary LDPC code defined over ${\mathrm{GF}}(2^{8})$.
The transmission takes place over the BEC. The compared accumulated LDPC
(ALDPC) codes are designed to achieve the capacity in the limit of large code
length. It can be seen that the proposed codes exhibit better decoding
performance than the ALDPC codes with code length up to 8192. The error floors
of the proposed codes can not be observed down to frame error rate $10^{-5}$
while the ALDPC codes have high error floors even with code length as long as
65536 bits.

Fig. 8 shows the thresholds of the proposed codes $C_{T}$ of rate $1/(3T)$ for
$T=1,\dotsc,20$. The codes are defined over ${\mathrm{GF}}(2^{8}).$ The
threshold values are calculated by the Monte Carlo simulation method. The
method was originally suggested in [31, p. 22] and an efficient calculation
was developed in [32, Section VII]. It can be observed that the proposed codes
leave a gap to the ultimate Shannon limit
$\mathrm{E_{b}/N_{0}}=\log_{10}(\ln(2))\approx-1.59$ [dB] [8] even in the
limit of large repetition parameter $T$. Fig. 9 depicts the simulation results
for the AWGN channels of $C_{2}$ and $C_{7}$ of very long code length. We
observe the convergence to a sharp threshold effect at the predicted threshold
values as the information length $k$ increases.

We demonstrate the decoding performance of the short and moderate-length
proposed codes $C_{T}$ for $T=1,2,3,4,6$ over the binary-input AWGN channels.
The mother code $C_{1}$ is constructed by the optimization method in [15]. The
coefficients $r_{N+1},\dotsc,r_{TN}$ are chosen uniformly at random from
${\mathrm{GF}}(2^{m})\setminus\\{0,1\\}$, where $1\in{\mathrm{GF}}(2^{m})$ is
the multiplicative identity. We fix $m=8$ for its good performance and the
computer-friendly representation of byte.

Fig. 10 shows the decoding performance of $C_{T}$ for $T=1,2,3,6$ of rates
${1}/{(3T)}$. It also shows a hybrid non-binary LDPC code [16] of rate 1/6 and
punctured $C_{1}$ of rate $1/2$. All these codes have 192 information bits.
The proposed code $C_{2}$ outperforms the hybrid non-binary LDPC code which is
the best code so far for that rate and code length. The code $C_{3}$ of rate
1/9 has about 0.5 [dB] coding gain from $C_{2}$ of rate 1/6. As we show on
these curves, the proposed construction, although simple, allows to design
codes with very low rates without large loss of gap to the Shannon limits.

The same property can be seen for the proposed codes with larger information
bits. Fig. 11 shows the decoding performance of $C_{T}$ for $T=1,2,3,6$, the
binary multi-edge type LDPC code of rate 1/2 and 1/6, and the binary ARA code
[6] of rate 1/6. All these codes have 1024 information bits except that the
MET LDPC code of rate 1/2 has 1280 information bits. It also shows the
performance of a punctured mother code $C_{1}$ of rate 1/2. Among the codes of
rate 1/6, the proposed code $C_{2}$ has the best performance both at water-
fall and error-floor regions.

As posed in Problem 3, conventional low-rate codes required large code length
to exploit the potential performance. It can be seen that the proposed codes
exhibit better decoding performance both at small and moderate code length.

##  VI Conclusions

We propose non-binary LDPC codes concatenated with inner multiplicative
repetition codes. The performance of the proposed codes exceeds the hybrid
non-binary codes, multi-edge type LDPC codes and ARA codes both at the water-
fall and error-floor regions. The encoder and decoder are inherently rate-
compatible, and especially the decoder complexity is almost the same as the
mother code.

## Acknowledgments

The authors are grateful to T. Richardson for providing the data of the multi-
edge type LDPC code of rate 1/6 in Fig. 11. The authors would like to thank
anonymous reviewers of ISIT2010 and Transactions on Information Theory, and I.
Sason for their suggestions and comments. K. K. wishes to thank T. Uyematsu
for valuable comments.

## References

  * [1] K. Kasai, D. Declercq, C. Poulliat, and K. Sakaniwa, “Rate-compatible non-binary LDPC codes concatenated with multiplicative repetition codes,” in _Proc. 2010 IEEE Int. Symp. Inf. Theory (ISIT)_ , Jun. 2010, pp. 844 –848.
  * [2] R. G. Gallager, _Low Density Parity Check Codes_. in Research Monograph series, MIT Press, Cambridge, 1963\.
  * [3] T. Richardson and R. Urbanke, “The capacity of low-density parity-check codes under message-passing decoding,” _IEEE Trans. Inf. Theory_ , vol. 47, no. 2, pp. 599–618, Feb. 2001.
  * [4] T. J. Richardson, M. A. Shokrollahi, and R. L. Urbanke, “Design of capacity-approaching irregular low-density parity-check codes,” _IEEE Trans. Inf. Theory_ , vol. 47, no. 2, pp. 619–637, Feb. 2001.
  * [5] T. Richardson and R. Urbanke, “Multi-edge type LDPC codes,” 2003.
  * [6] D. Divsalar, S. Dolinar, and C. Jones, “Low-rate LDPC codes with simple protograph structure,” in _Proc. 2005 IEEE Int. Symp. Inf. Theory (ISIT)_ , Sep. 2005, pp. 1622–1626.
  * [7] G. Yue, L. Ping, and X. Wang, “Generalized low-density parity-check codes based on Hadamard constraints,” _IEEE Trans. Inf. Theory_ , vol. 53, no. 3, pp. 1058–1079, Mar. 2007.
  * [8] S. Haykin, _Communication Systems_ , 4th ed. John Wiley & Sons, 2001.
  * [9] M. Davey and D. MacKay, “Low-density parity check codes over GF($q$),” _IEEE Commun. Lett._ , vol. 2, no. 6, pp. 165–167, Jun. 1998.
  * [10] W. Chang and J. Cruz, “Nonbinary LDPC codes for 4-kB sectors,” _IEEE Trans. Magn._ , vol. 44, no. 11, pp. 3781–3784, Nov. 2008.
  * [11] I. Djordjevic and B. Vasic, “Nonbinary LDPC codes for optical communication systems,” _Photonics Technology Letters, IEEE_ , vol. 17, no. 10, pp. 2224–2226, Oct. 2005.
  * [12] B. Zhou, J. Kang, S. Song, S. Lin, K. Abdel-Ghaffar, and M. Xu, “Construction of non-binary quasi-cyclic LDPC codes by arrays and array dispersions,” _IEEE Trans. Commun._ , vol. 57, no. 6, pp. 1652–1662, Jun. 2009.
  * [13] M. Arabaci, I. Djordjevic, R. Saunders, and R. Marcoccia, “High-rate nonbinary regular quasi-cyclic LDPC codes for optical communications,” _Lightwave Technology, Journal of_ , vol. 27, no. 23, pp. 5261–5267, Dec. 2009.
  * [14] B. Zhou, J. Kang, Y. Tai, S. Lin, and Z. Ding, “High performance non-binary quasi-cyclic LDPC codes on euclidean geometries LDPC codes on euclidean geometries,” _IEEE Trans. Commun._ , vol. 57, no. 5, pp. 1298–1311, May 2009.
  * [15] C. Poulliat, M. Fossorier, and D. Declercq, “Design of regular (2,$d_{c}$)-LDPC codes over GF($q$) using their binary images,” _IEEE Trans. Commun._ , vol. 56, no. 10, pp. 1626–1635, Oct. 2008.
  * [16] L. Sassatelli, D. Declercq, and C. Poulliat, “Low-rate non-binary hybrid LDPC codes,” in _Proc. 5th Int. Symp. on Turbo Codes and Related Topics_ , Sep. 2008, pp. 225–230.
  * [17] V. Rathi and R. Urbanke, “Density Evolution, Threshold and the Stability Condition for non-binary LDPC Codes,” _IEE Proceedings - Communications_ , vol. 152, no. 6, pp. 1069–1074, 2005.
  * [18] F. J. MacWilliams and N. J. A. Sloane, _The Theory of Error-Correcting Codes_. Amsterdam: Elsevier, 1977.
  * [19] T. Richardson and R. Urbanke, _Modern Coding Theory_. Cambridge University Press, Mar. 2008.
  * [20] J. Justesen, “Class of constructive asymptotically good algebraic codes,” _IEEE Trans. Inf. Theory_ , vol. 18, no. 5, pp. 652–656, Sep. 1972.
  * [21] M. Davey and D. MacKay, “Low density parity check codes over GF($q$),” in _Information Theory Workshop, 1998_ , Jun. 1998, pp. 70–71.
  * [22] D. Declercq and M. Fossorier, “Decoding algorithms for nonbinary LDPC codes over GF($q$),” _IEEE Trans. Commun._ , vol. 55, no. 4, pp. 633–643, Apr. 2007.
  * [23] A. Bennatan and D. Burshtein, “On the application of LDPC codes to arbitrary discrete-memoryless channels,” _IEEE Trans. Inf. Theory_ , vol. 50, no. 3, pp. 417–438, Mar. 2004.
  * [24] G. Li, I. Fair, and W. Krzymien, “Density evolution for nonbinary LDPC codes under Gaussian approximation,” _IEEE Trans. Inf. Theory_ , vol. 55, no. 3, pp. 997 –1015, Mar. 2009.
  * [25] D. Klinc, J. Ha, and S. McLaughlin, “Optimized puncturing and shortening distributions for nonbinary LDPC codes over the binary erasure channel,” in _Proc. 46th Annual Allerton Conf. on Commun., Control and Computing_ , Sep. 2008, pp. 1053–1058.
  * [26] H. Pfister and I. Sason, “Accumulate repeat accumulate codes: Capacity-achieving ensembles of systematic codes for the erasure channel with bounded complexity,” _IEEE Trans. Inf. Theory_ , vol. 53, no. 6, pp. 2088 –2115, Jun. 2007.
  * [27] C. Di, D. Proietti, I. Telatar, T. Richardson, and R. Urbanke, “Finite-length analysis of low-density parity-check codes on the binary erasure channel,” _IEEE Trans. Inf. Theory_ , vol. 48, no. 6, pp. 1570 –1579, Jun. 2002.
  * [28] C. E. Shannon, “Capacity of the band-limited Gaussian channel,” _Bell System Tech. J._ , vol. 38, pp. 611–656, May 1959.
  * [29] A. Valembois and M. Fossorier, “Sphere-packing bounds revisited for moderate block lengths,” _IEEE Trans. Inf. Theory_ , vol. 50, no. 12, pp. 2998 – 3014, Dec. 2004.
  * [30] G. Wiechman and I. Sason, “An improved sphere-packing bound for finite-length codes over symmetric memoryless channels,” _IEEE Trans. Inf. Theory_ , vol. 54, no. 5, pp. 1962 –1990, May 2008.
  * [31] M. Davey, “Error-correction using low-density parity-check codes,” Ph.D. dissertation, Univ. Cambridge, Cambridge, U.K, Dec. 1999.
  * [32] L. Sassatelli and D. Declercq, “Nonbinary hybrid LDPC codes,” _IEEE Trans. Inf. Theory_ , vol. 56, no. 10, pp. 5314 –5334, Oct. 2010.


In this paper, we present two Hermite polynomial based approaches to derive one-step numerical integrators for mechanical systems. These methods are based on discretizing the configuration using Hermite polynomials which leads to numerical trajectories continuous in both configuration and velocity. First, we incorporate Hermite polynomials for time-discretization and derive one-step variational methods by discretizing the Lagrange-d'Alembert principle over a single time step. Second, we present the Galerkin approach to derive one-step numerical integrators by setting the weighted average of the residual of the equations of motion over a time step to zero.

We consider three numerical examples to understand the numerical performance of the one-step variational and Galerkin methods. We first study a particle in a double-well potential and compare the variational approach results with the corresponding results for the Galerkin approach. We then study the Duffing oscillator to understand the numerical behavior in presence of dissipative forces. Finally, we apply the proposed methods to a nonlinear aeroelastic system with two degrees of freedom. Both variational and Galerkin one-step methods capture conservative and nonconservative dynamics accurately with excellent energy behavior. The one-step Galerkin methods exhibit better trajectory and energy performance than the one-step variational methods and the variational integrators.

Keywords: Structure-preserving numerical methods; One-step time-integrators; Variational methods; Galerkin methods; Conservative and dissipative mechanical systems.

§ INTRODUCTION

The field of numerical integration of  ordinary differential equations (ODEs) reached scientific maturity towards the end of the 20th century when general purpose codes based on Runge-Kutta or linear multistep methods became widely available. In fact these off-the-shelf numerical integration packages have become an essential tool in studying complex dynamical behavior of nonlinear systems. Over the past decades, with growing interest in long-time simulation of dynamical systems from fields such as astronomy or molecular dynamics, the geometric properties of the governing differential equations became crucial for numerical simulation. This need for special numerical methods for certain classes of problems led to the advent of the field of geometric numerical integration (GNI) methods. Unlike the traditional methods that are designed to give computationally inexpensive solutions with small overall error, GNI methods are designed to preserve the underlying geometric structure of the governing differential equations. Motivation for applying these structure-preserving numerical methods to engineering problems ranges from preservation of qualitative features, such as invariants of motion or configuration space structure, to long-time stability - accurate numerical simulation for exponentially long times.  The essential aspects of GNI methods for ODEs are presented in [1, 2].

Many mechanical systems of interest to engineers possess physically meaningful invariants such as momentum, energy, vorticity or rotational symmetry. Among such mechanical systems, Lagrangian/Hamiltonian systems form the most important class of ODEs in the context of structure-preserving numerical integration. One of the key features of Lagrangian/Hamiltonian mechanical systems is that the governing equations can be derived from variational principles. Wendlandt and Marsden [3] derived mechanical integrators for conservative dynamical systems by discretizing Hamilton's principle, instead of discretizing the governing equations. These numerical integrators are both symplectic and momentum-preserving. These numerical integrators, because of the variational nature of their derivation, are known as variational integrators. Building on the ideas of Maeda, Moser and Veselov [4, 5, 6, 7, 8], Marsden and his colleagues [9] developed the discrete mechanics framework and extended these methods to nonconservative [10] and/or constrained mechanical systems [11]. Leok and Shingel [12] used Hermite polynomials in the discrete mechanics setting to derive variational integrators that approximate time derivatives of the trajectory with accuracy. Kane et al. [13] developed symplectic-energy-momentum integrators for conservative Lagrangian systems by imposing an additional discrete energy equation to compute the adaptive time step. Sharma et al. [14] extended these energy-preserving, adaptive time step variational integrators to Lagrangian systems with external forcing.  A recent review paper by Sharma et al. [15] provides an overview of variational integrators with focus on their applications in different engineering fields.

Alternatively, the structural dynamics community has used time finite elements to derive numerical integrators from the variational principle since the late 1960s. Argryis and Scharpf [16] used cubic Hermite shape functions to formulate time finite elements and discretize Hamilton's principle to obtain numerical integrators for initial value problems. Baruch and Riff [17] presented six different formulations of time finite elements based on Hamilton's law of varying action, and later on, Riff and Baruch [18] implemented numerical integrators based on these formulations. Recently, Mergel et al. [19] developed $\mathcal{C}^1-$continuous time integrators based on Hamilton's law of varying action for conservative systems. Unlike the discrete mechanics approach used for deriving variational integrators, the time finite elements approach only considers the variational principle over one time step and uses Hermite polynomials for discretizing the configuration which leads to $\mathcal{C}^1$–continuous trajectories. Apart from the accuracy and computational stability, discretization uses Hermite polynomials also leads to numerical methods that are amenable to feedback control implementation and analysis since they naturally yield configuration states and velocities.

Most of these time finite elements approaches based on Hamilton's principle or Hamilton's law of varying action fit into the framework of continuous Galerkin  or discontinuous Galerkin method. Hulme [20] employed the so-called continuous Galerkin method to derive one-step methods for numerically solving systems of nonlinear first-order ODEs. Betsch and Steinmann [21, 22]  applied Petrov-Galerkin based time finite elements to the Hamiltonian formulation of the $N$-Body problem and semidiscrete nonlinear elastodynamics.

The purpose of this paper is twofold. First, we derive Hermite polynomial based $\mathcal{C}^1$–continuous numerical integrators for mechanical systems from two different approaches. For the variational methods, we discretize the variational principle over a single fixed time step and then use the discrete mechanics framework to develop one-step variational methods. We also use Galerkin's method of weighted residuals of the governing equations over a single fixed time step to develop one-step Galerkin methods. Second, we study the numerical performance of the developed methods for three different classes of mechanical systems. We  also investigate the linear stability and the symplecticity of both one-step methods.

The remainder of the paper is organized as follows. In Section <ref> we review the basic concepts from variational integrators and Galerkin methods used later in this work. In Section <ref> we present Hermite polynomial based one-step variational and Galerkin methods.  First, we introduce cubic Hermite polynomials used for time-discretization in this work. Then we utilize the discrete mechanics approach to derive the one-step variational methods for mechanical systems with external forcing. Finally, we use the Galerkin approach to derive one-step Galerkin methods. In Section <ref> we investigate the linear stability and symplectic nature of the proposed one-step methods. In Section <ref> we study three numerical examples to understand the numerical performance of the proposed one-step methods. Finally, in Section <ref> we provide concluding remarks and suggest future directions for this work.

§ BACKGROUND

The numerical integration of mechanical systems can be approached in different ways. Traditional methods apply time-discretization directly to the governing equations of motion to obtain time-integration algorithms. Unfortunately, this approach does not account for the qualitative properties of the dynamical system. In this section, we review the basic concepts from variational mechanics and Galerkin methods used in the development of our one-step variational and Galerkin methods.
Consider a time-invariant Lagrangian mechanical system with a finite-dimensional, smooth configuration manifold $\Q$, state space $\T\Q$, and Lagrangian $L: \T\Q \to \mathbb{R}$.
For such an autonomous Lagrangian system with time-independent external forcing $\textbf{f}_L (\q(t),\dot{\q}(t))$, the Lagrange-d'Alembert principle characterizes trajectories $\q(t)$ as those satisfying

\begin{equation}
\delta \ \int^{t_{\rm f}}_{t_{\rm i}} L(\q(t),\dot{\q}(t)) \ \text{d}t + \int ^{t_{\rm f}}_{t_{\rm i}} \textbf{f}_L (\q(t),\dot{\q}(t)) \cdot \delta \q \ \text{d} t  =0,
\label{eq:lda}
\end{equation}

where the first term considers the variation of the action integral and the second term accounts for the virtual work done by the external forces when the path $\q(t)$ is varied by $\delta \q(t)$. Using integration by parts and setting the variations at the endpoints equal to zero gives the forced Euler-Lagrange equations

\begin{equation}
\textbf{M}_{\rm eq}(\q(t),\dot{\q}(t),\ddot{\q}(t))=\frac{\partial L(\q(t),\dot{\q}(t))}{\partial \q}   - \frac{\text{d}}{\text{d}t} \left( \frac{\partial L(\q(t),\dot{\q}(t))}{\partial \dot{\q}} \right) + \textbf{f}_L(\q(t),\dot{\q}(t)) =\textbf{0},
\label{eq:fg}
\end{equation}

where $\textbf{M}_{\rm eq}$ denotes equations of motion written in a residual form. For a general mechanical system with a separable Lagrangian of the form $L(\q(t),\dot{\q}(t))=\frac{1}{2}\dot{\q}(t)^\top \mathbb{M}(\q)\dot{\q}(t)-U(\q)$, the equations of motion are given by

\begin{equation}
\textbf{M}_{\rm eq}(\q(t),\dot{\q}(t),\ddot{\q}(t))=\mathbb{M}(\q)\frac{\text{d}^2\q}{\text{d}t^2} + \frac{\partial U(\q)}{\partial \q} - \textbf{f}_L(\q(t),\dot{\q}(t)) = \textbf{0},
\end{equation}

where $\mathbb{M}(\q)$ is the mass matrix and $U(\q)$ is the potential energy of the mechanical system.

§.§ Discrete Variational Mechanics and Variational Integrators

Variational integrators are time-integration methods that utilize concepts from discrete mechanics,
a discrete analogue to continuous-time variational mechanics. Although these methods were
originally developed for conservative dynamical systems, due to the variational nature of their
construction, these methods can be easily extended to nonconservative mechanical systems by
discretizing the corresponding Lagrange-d’Alembert principle. The basic idea is to first construct
discrete-time approximations of both the action integral and virtual work terms in (<ref>) and then apply a discrete variational principle to obtain discrete-time trajectories of the mechanical system.

Consider a discrete Lagrangian system with configuration manifold $\Q$ and discrete state space $\Q \times \Q$. For a fixed time step $\dt=\frac{t_{\rm f}-t_{\rm i}}{N}$, the discrete trajectory $\{\ \q_k  \  \}_{k=0}^N\ $ is defined by the configuration of the system at the sequence of times $\{\ t_k= t_{\rm i} + k\dt~|~k=0,...,N\}\ $.
We introduce the discrete Lagrangian function $L_{\rm d}(\q_k,\q_{k+1})$, an approximation of the action integral along the curve from $\q_k$ to $\q_{k+1}$, which approximates the integral of the Lagrangian in the following sense

\begin{equation}
L_{\rm d}(\q_k,\q_{k+1}) \approx \int^{t_{k+1}}_{t_k} L(\q(t),\dot{\q}(t)) \  \text{d}t.
\label{eq:dl}
\end{equation}

To discretize external forcing, we define two discrete forces $\textbf{f}_{\rm d}^{\pm} :  \Q \times \Q \to \T^*\Q$ which approximate the continuous-time force integral that appears in (<ref>) over one time step in the following sense

\begin{equation}
\textbf{f}_{\rm d}^+(\q_k,\q_{k+1})\cdot\delta \q_{k+1} +  \textbf{f}_{\rm d}^-(\q_k,\q_{k+1})\cdot\delta \q_{k}  \approx \int ^{t_{k+1}}_{t_k} \textbf{f}_L (\q(t), \dot{\q}(t)) \cdot \delta \q \ \text{d}t.
\end{equation}

The discrete Lagrange-d'Alembert principle seeks $\{\ \q_k  \  \}_{k=0}^N\ $ that satisfy

\begin{equation}
\delta \sum_{k=0}^{N-1} L_{\rm d}(\q_k,\q_{k+1}) + \sum_{k=0}^{N-1}[\textbf{f}_{\rm d}^+(\q_k,\q_{k+1})\cdot\delta \q_{k+1} +  \textbf{f}_{\rm d}^-(\q_k,\q_{k+1})\cdot\delta \q_{k}] = 0,
\end{equation}

which yields the following forced discrete Euler-Lagrange equations

\begin{equation}
\frac{\partial L_{\rm d}(\q_{k-1},\q_{k})}{\partial \q_k} + \frac{\partial L_{\rm d}(\q_k,\q_{k+1})}{\partial \q_{k}} + \textbf{f}_{\rm d}^+(\q_{k-1},\q_k) + \textbf{f}_{\rm d}^-(\q_{k},\q_{k+1}) = \textbf{0}   \ \ \ \ k=1,..., N-1.
\end{equation}

These discrete equations can be recast in a standard time-marching form as follows

\begin{align}
\label{eqn:vif_imp}
-\frac{\partial L_{\rm d}(\q_k,\q_{k+1})}{\partial \q_k} - \textbf{f}_{\rm d}^-(\q_k,\q_{k+1}) &=\textbf{p}_k, \\
\label{eqn:vif_exp}
\textbf{p}_{k+1}&=\frac{\partial L_{\rm d}(\q_{k},\q_{k+1})}{\partial \q_{k+1}} + \textbf{f}_{\rm d}^+(\q_k,\q_{k+1}),
\end{align}

where $\textbf{p}_k$ is the discrete momentum corresponding to the discrete configuration $\q_k$.
Due to the variational nature of their derivation, discrete trajectories from variational integrators inherit geometric properties from their continuous-time counterpart. These algorithms are ideal for long-time simulation because of their numerical stability and excellent energy behavior over exponentially long times for conservative as well as nonconservative systems.

§.§ Galerkin Methods for Numerical Integration

Galerkin methods are a class of methods used for converting continuous differential equation problems to discrete problems. The basic idea behind these methods is to seek approximate solutions to the differential equation in a finite-dimensional space spanned by a set of basis functions.
Consider a Hilbert space $U$, and a bilinear form $a : U \times V \to \Real$ which is both bounded and V-eliptic. We consider the following abstract problem posed as a weak formulation on the Hilbert spaces $U$ and $V$, namely,

\begin{equation}
Find\quad u \in U\quad  s.t \quad \quad a(u,v)=\ell(v)  \quad \forall v \in V,
\label{eq:weak}
\end{equation}

where $\ell$ is a bounded linear functional on $V$. In general, it is very rare to find an exact solution of (<ref>) because $U$ and $V$ are infinite-dimensional. A natural approach to construct an approximate solution is to solve a finite-dimensional analogue of (<ref>). Let $U_N \subseteq U$ and $V_N \subseteq V$ be $N-$dimensional subspaces. We project the original problem onto $U_n$, i.e.,

\begin{equation}
Find\quad u_n \in U_n\quad  s.t \quad \quad a(u_n,v_n)=\ell(v_n)  \quad \forall v_n \in V_n.
\label{eq:galerkin}
\end{equation}

Reducing the original problem to a $N-$dimensional subspace allows us to numerically solve for $u_n$ as a finite linear combination of basis vectors in $U_n$. For a finite number of basis functions, the Galerkin methods lead to a system of equations with a finite number of unknowns. Depending on the choice of test functions, these methods can be classified into different methods such as Bubnov-Galerkin, Petrov-Galerkin, collocation methods, or the least squares method.
Although the Galerkin methods can be used for solving a wide variety of problems, we restrict our attention to using ideas from Galerkin methods to develop time integrators for dynamical systems. Consider an autonomous system of first-order ODEs given by

\begin{equation}
\frac{\text{d}\y}{\text{d}t} = \f(\y), \quad \quad \y(t_0)=\y_0,
\end{equation}

on a finite time interval $t \in [t_0,t_f]$ where $\y(t_0)=\y_0$ is the initial condition for the given initial value problem. We rewrite the governing first-order equations in the residual form by defining

\begin{equation}
\R\left(\y,\frac{\text{d}\y}{\text{d}t}\right)= \frac{\text{d}\y}{\text{d}t} - \f(\y).
\end{equation}

We discretize the total time interval into $N$ subintervals of fixed time step size $\dt=\frac{t_f-t_0}{N}$. On each subinterval $[t_k, t_{k+1}]$, we consider piecewise smooth polynomial approximations of $\y$

\begin{equation}
\y(t) \approx \y_{\rm{d},k}(t) = \sum_{i=1}^{n+1} b^{(k)}_i \Phi_{i,k}(t),
\end{equation}

where $\y_{\rm{d},k}(t)$ is the approximation over the $k$th subinterval and  $\Phi_{i,k}(t)$ are basis functions that are $n$th degree polynomials on $k$th subinterval.  Since the coefficients $b_i^{(k)}$ may change from one subinterval to the next, the global approximation $\y_{\rm d}(t)$ need not be as smooth as the trial basis functions $\Phi_{i,k}(t)$. Using this piecewise polynomial approximation, we can derive a $\mathcal{C}^0-$continuous numerical integrator by requiring that $\y_{\rm{d},k}(t)$ is a local Galerkin approximation on the $k$th subinterval  and $\y_{\rm d}$ is continuous on $[t_0,t_f]$. For example, if we use $w_i=\Phi_{i,k}$ as test functions then that leads to solving the following $n+1$ equations for $k=0,\cdots,N-1$

\begin{equation}
\y_{\rm{d},k}(t_{k})=\begin{cases}\y_{\rm{d},k-1}(t_{k}), \quad & k\geq 1 \\ \y_0,\quad & k=0 \end{cases} \end{equation}

\begin{equation}
\bigg \langle \textbf{R}\left(\y_{\rm d},\frac{\text{d}\y_{\rm d}}{\text{d}t}\right), \Phi_{i,k} \bigg \rangle_k = 0, \quad i=2, \cdots,n+1,
\end{equation}

where $\langle v,w\rangle_k=\int_{t_k}^{t_{k+1}}v(t)w(t)\text{d}t$ is the inner product on the $k$th subinterval.

§ ONE-STEP TIME-INTEGRATION METHODS

§.§ Hermite Polynomials

We discretize the continuous trajectory over one time step using cubic Hermite polynomials

\begin{equation}
\q(t)\approx \q_{\rm d}(t) = \q_0N_1(t) + \textbf{v}_0N_2(t) +  \q_{1}N_3(t) + \textbf{v}_{1}N_4(t),
\label{eq:hermite}
\end{equation}

where $\q_{\rm d}(t)$ is the discrete approximation and the Hermite polynomials are given by

\begin{multline*}
N_1(t)=2 \left( \frac{t}{\dt} \right)^3 - 3\left( \frac{t}{\dt} \right)^2  + 1, \hfill  N_3(t)=-2 \left( \frac{t}{\dt} \right)^3 + 3\left( \frac{t}{\dt} \right)^2, \\  N_2(t) = \dt \left[  \left( \frac{t}{\dt} \right)^3 - 2\left( \frac{t}{\dt}   \right)^2 + \left( \frac{t}{\dt}   \right)  \right], \hfill  N_4(t) = \dt \left[  \left( \frac{t}{\dt} \right)^3 - \left( \frac{t}{\dt}   \right)^2  \right],
\end{multline*}

where $\dt$ is the fixed time step. As shown in Figure <ref>, at initial time $t=0$, we have $N_1(0)=1$ along with $N_2(0)=N_3(0)=N_4(0)=0$ which leads to $\q_{\rm d}(0)=\q_0$. Similarly, we also have $\dot{\q}_{\rm d}(0)=\textbf{v}_0$, $\q_{\rm d}(\dt)=\q_1$ and $\dot{\q}_{\rm d}(\dt)=\textbf{v}_1$. Using piecewise Hermite polynomials for discretization in the one-step approach leads to numerical solutions that are $\mathcal{C}^1$– continuous.

Cubic Hermite polynomials plotted over one time step.

§.§ One-step Variational Methods

Unlike the discrete mechanics framework developed by Marsden and West [9], we only consider the variational principle over one time step $\dt$ and derive one-step numerical integrators from the discretized variational principle.

We use the discrete approximation $\q_{\rm d}(t)$ from (<ref>) to obtain a discrete action $S_{\rm d}$ which approximates the action integral over one time step in the following sense

\begin{equation}
S_{\rm d}(\q_0,\textbf{v}_0,\q_{1},\textbf{v}_{1}) = \int_{0}^{\dt} L(\q_{\rm d},\dot{\q}_{\rm d})  \  \text{d}t \approx \int_{0}^{\dt} L(\q,\dot{\q})  \  \text{d}t.
\end{equation}

Thus, we have used a cubic Hermite polynomial to obtain the discrete action $S_{\rm d}$ in terms of the configuration and velocity at the endpoints. We introduce discrete forces corresponding to displacement and velocity variations to approximate the virtual work

\begin{multline}
\int ^{\dt}_{0} \textbf{f}_L (\q(t), \dot{\q}(t)) \cdot\delta \q \ \text{d}t \approx \int  ^{\dt}_{0} \textbf{f}_L (\q_{\rm d}(t), \dot{\q}_{\rm d}(t)) \cdot\delta \q_{\rm d} \ \text{d}t \\
=\textbf{f}_{\rm d}^+(\q_0,\textbf{v}_0,\q_{1},\textbf{v}_{1})\cdot\delta \q_{1} +  \textbf{f}_{\rm d}^-(\q_0,\textbf{v}_0,\q_{1},\textbf{v}_{1})\cdot\delta \q_{0} \ \ \ \ \ \ \ \ \  \ \ \ \ \   \\
+ \textbf{g}_{\rm d}^+(\q_0,\textbf{v}_0,\q_{1},\textbf{v}_{1})\cdot\delta \textbf{v}_{1} +  \textbf{g}_{\rm d}^-(\q_0,\textbf{v}_0,\q_{1},\textbf{v}_{1})\cdot\delta \textbf{v}_{0}.
\end{multline}

Using $\delta \q_{\rm d}(t) = \delta \q_0N_1(t) + \delta \textbf{v}_0N_2(t) +  \delta \q_{1}N_3(t) + \delta \textbf{v}_{1}N_4(t)$, we can obtain the discrete forces $\textbf{f}_{\rm d}^{\pm}$ corresponding to displacement variations

\begin{equation}
\textbf{f}_{\rm d}^+(\q_0,\textbf{v}_0,\q_{1},\textbf{v}_{1}) = \int  ^{\dt}_{0} \textbf{f}_L (\q_{\rm d}(t), \dot{\q}_{\rm d}(t)) \ N_3(t) \ \text{d}t,
\end{equation}

\begin{equation}
\textbf{f}_{\rm d}^-(\q_0,\textbf{v}_0,\q_{1},\textbf{v}_{1}) = \int  ^{\dt}_{0} \textbf{f}_L (\q_{\rm d}(t), \dot{\q}_{\rm d}(t)) \ N_1(t) \ \text{d}t,
\end{equation}

and discrete forces $\textbf{g}_{\rm d}^{\pm}$ corresponding to velocity variations

\begin{equation}
\textbf{g}_{\rm d}^+(\q_0,\textbf{v}_0,\q_{1},\textbf{v}_{1}) = \int  ^{\dt}_{0} \textbf{f}_L (\q_{\rm d}(t), \dot{\q}_{\rm d}(t)) \ N_4(t) \ \text{d}t,
\end{equation}

\begin{equation}
\textbf{g}_{\rm d}^-(\q_0,\textbf{v}_0,\q_{1},\textbf{v}_{1}) = \int  ^{\dt}_{0} \textbf{f}_L (\q_{\rm d}(t), \dot{\q}_{\rm d}(t)) \ N_2(t) \ \text{d}t.
\end{equation}

The discrete Lagrange-d'Alembert principle  using the one-step variational approach seeks $ \q_{\rm d}(t)$  (<ref>) that satisfy

\begin{multline*}
\delta S_{\rm d}(\q_0,\textbf{v}_0, \q_{1},\textbf{v}_{1}) + [\textbf{f}_{\rm d}^+(\q_0,\textbf{v}_0, \q_{1},\textbf{v}_{1})\cdot\delta \q_{1} +  \textbf{f}_{\rm d}^-(\q_0,\textbf{v}_0, \q_{1},\textbf{v}_{1})\cdot\delta \q_{0}] \\ + [\textbf{g}_{\rm d}^+(\q_0,\textbf{v}_0, \q_{1},\textbf{v}_{1})\cdot\delta \textbf{v}_{1} +  \textbf{g}_{\rm d}^-(\q_0,\textbf{v}_0, \q_{1},\textbf{v}_{1})\cdot\delta \textbf{v}_{0}]  = 0,
\end{multline*}

which gives

\begin{multline*}
\left(  \frac{\partial S_{\rm d}}{\partial \q_0} + \textbf{f}_{\rm d}^-(\q_0,\textbf{v}_0, \q_{1},\textbf{v}_{1}) \right) \cdot \delta \q_0 + \left(  \frac{\partial S_{\rm d}}{\partial \textbf{v}_0} + \textbf{g}_{\rm d}^-(\q_0,\textbf{v}_0, \q_{1},\textbf{v}_{1}) \right) \cdot \delta \textbf{v}_0 \\ +  \left(  \frac{\partial S_{\rm d}}{\partial \q_1} + \textbf{f}_{\rm d}^+(\q_0,\textbf{v}_0, \q_{1},\textbf{v}_{1}) \right) \cdot \delta \q_1 + \left(  \frac{\partial S_{\rm d}}{\partial \textbf{v}_1} + \textbf{g}_{\rm d}^+(\q_0,\textbf{v}_0, \q_{1},\textbf{v}_{1}) \right) \cdot \delta \textbf{v}_1  =0,
\end{multline*}

where $S_{\rm d}:=S_{\rm d}(\q_0,\textbf{v}_0, \q_{1},\textbf{v}_{1}).$ Setting variations at endpoints to zero, i.e. $\delta \q_0= \delta \q_1=\textbf{0}$, gives

\begin{equation}
\frac{\partial S_{\rm d}}{\partial \textbf{v}_0}(\q_0,\textbf{v}_0, \q_{1},\textbf{v}_{1}) + \textbf{g}_{\rm d}^-(\q_0,\textbf{v}_0, \q_{1},\textbf{v}_{1})=\textbf{0},
\end{equation}

\begin{equation}
\frac{\partial S_{\rm d}}{\partial \textbf{v}_1}(\q_0,\textbf{v}_0, \q_{1},\textbf{v}_{1}) + \textbf{g}_{\rm d}^+(\q_0,\textbf{v}_0, \q_{1},\textbf{v}_{1}) =\textbf{0}.
\end{equation}

Given $(\q_0,\textbf{v}_0)$, these coupled nonlinear equations can be solved to obtain $(\q_1,\textbf{v}_1)$. Thus, for Lagrangian systems with external forcing the one-step variational approach can be used to derive numerical integrators that are continuous in both configuration and velocities.

It is important to note that we have used the time finite elements approach [16] where considering the discrete Lagrange-d'Alembert principle over a single time step plays a crucial role in ensuring that the numerical integrators are $\mathcal{C}^1-$continuous and stable. Riff and Baruch [23] have shown that numerical integrators derived by considering the discrete action sum over the entire time interval are unconditionally unstable.

In [12], Hermite-based prolongation-collocation variational integrators (PCVIs) are constructed from discrete Lagrangian $L_{\rm d}(\q_0,\q_1)$ by means of expressing every parameter in $\q_{\rm d}(t)$ as function of $(\q_0,\q_1)$ and with the help of some extra equations based on the prolongation-collocation approach. PCVIs are fundamentally different from the one-step variational methods developed here in two ways. First, one-step variational methods proposed in this work are $\mathcal{C}^1-$continuous whereas the PCVIs in [12] are $\mathcal{C}^0-$continuous. Second, PCVIs are constructed by summing up the discrete action over the entire time interval $T$ whereas the one-step variational methods are constructed by considering the discrete action $S_{\rm d}(\q_0,\textbf{v}_0, \q_{1},\textbf{v}_{1})$ over a single time step $\Delta t$.

§.§ One-step Galerkin Methods

In this subsection, we consider the Petrov-Galerkin method or the weighted residual method to obtain numerical integrators for mechanical systems with external forcing. For a given system of equations we first write it in the residual form $\textbf{R}\left(t,\q_{\rm d},\frac{\text{d}\q_{\rm d}}{\text{d}t},...,\frac{\text{d}^n\q_{\rm d}}{\text{d}t^n}\right)$ and this continuous system of ODEs is transformed into following discrete equations

\begin{equation}
\bigg \langle \textbf{R}\left(t,\q_{\rm d},\frac{\text{d}\q_{\rm d}}{\text{d}t},...,\frac{\text{d}^n\q_{\rm d}}{\text{d}t^n}\right), w_i \bigg \rangle = 0, \quad i=1,2,...,N,
\end{equation}

where $\q_{\rm d}(t)$ is the assumed solution form and  $w_i$ are test functions. For our study, we focus on time-integration of mechanical problems with $\textbf{R}=\textbf{M}_{\rm eq}(\q,\dot{\q},\ddot{\q})$ where $\textbf{M}_{\rm eq}$ are the equations of motion for the mechanical system. Similar to the variational approach discussed in the previous section, we use cubic Hermite polynomials as solution functions for approximating the continuous solution over one time step. Our goal is to use the Petrov-Galerkin approach to derive one-step methods so we consider the following shifted Legendre polynomials as test functions

\begin{equation*}
P_0(t) =1, \quad P_1(t)= \frac{1}{\dt} (2t-\dt),
\end{equation*}

Given $(\q_0,\textbf{v}_0)$, the one-step Galerkin method yields

\begin{equation*}
\int_0^{\dt} \textbf{M}_{\rm eq}(\q_{\rm d}(t),\dot{\q}_{\rm d}(t),\ddot{\q}_{\rm d}(t))\left( 1 \right) \text{d}t=\textbf{0},
\end{equation*}

\begin{equation*}
\int_0^{\dt} \textbf{M}_{\rm eq}(\q_{\rm d}(t),\dot{\q}_{\rm d}(t),\ddot{\q}_{\rm d}(t))\left(\frac{1}{\dt}(2t-\dt)\right)\text{d}t=\textbf{0}.
\end{equation*}

Given  $(\q_0,\textbf{v}_0)$, these two coupled nonlinear equations can be solved to obtain the configuration $\q_1$ and velocity $\textbf{v}_1$ at the next time step. Just like the variational approach, this system of nonlinear equations can be used as a one-step numerical integrator.

§.§ One-step Methods based on Higher-order Hermite Polynomials

In this subsection, we demonstrate how to derive one-step methods with higher-order Hermite polynomials. The one-step methods presented so far have been based on cubic Hermite polynomials which lead to numerical integrators that are continuous in both configuration and velocities. For discretization using higher-order Hermite polynomials, we consider the following discrete trajectory over one time step

\begin{equation}
\q_{\rm d}(t)=\sum_{j=0}^{n-1} \left( \q^{(j)}_0 H_{n,j}(t) + \q^{(j)}_1 H_{n,j}(\dt-t) \right),
\end{equation}

where we have written the discrete trajectory in terms of Hermite basis functions and values of $\q(t)$ and its derivatives $\q^{(j)}(t)$ at endpoints of each interval. The Hermite basis functions are

\begin{equation}
H_{n,j}(t)= \frac{t^j}{j!} \left(1-\frac{t}{\dt}\right)^n \sum_{s=0}^{n-j-1} \begin{pmatrix}
n+s-1 \\
\end{pmatrix}\left( \frac{t}{\dt} \right)^s.
\end{equation}

Thus, the discrete trajectory $\q_{\rm d}$ is represented by a $2n-1$ degree polynomial which satisfies

\begin{equation}
\q^{(j)}(0)=\q_{\rm d}^{(j)}(0), \quad \quad \quad  \q^{(j)}(\dt)=\q_{\rm d}^{(j)}(\dt),   \quad \quad j=0,...,n-1.
\end{equation}

It is clear from the above expression that for  $n=1$, the discrete trajectory simply reduces to a linear interpolant between endpoints $\q_0$ and $\q_1$. For $n=2$, the discrete trajectory is the cubic interpolant discussed in Section <ref>. For discretization using higher-order Hermite polynomials with $n \geq 3$, the discrete trajectory $\q_{\rm d}$ over the fixed time step is represented by $2n-1$ degree polynomials with $2n$ unknown coefficients. For an initial condition of the form $(\q_0,\textbf{v}_0)$, the numerical integration problem reduces to solving for the remaining $2n-2$ coefficients. The first $n-2$ coefficients are $\q^{j}(0)$ for $j=2,\cdots, n-1$ and the other $n$ coefficients are $\q^{j}(\dt)$ for $j=0,\cdots, n-1$.

Variational Approach: For conservative Lagrangian systems, the discrete Hamilton's principle after setting the configuration variations at the endpoints to zero (i.e. $\delta \q_0=\delta \q_1=\textbf{0}$) leads to the following discrete equations

\begin{equation}
\frac{\partial S_{\rm d}}{\partial \q^{(j)}(0)}= \frac{\partial S_{\rm d}}{\partial \q^{(j)}(\dt)} = \textbf{0}, \quad \quad j=1,\cdots, n-1.
\end{equation}

Thus, solving these $2n-2$ coupled nonlinear equations gives the $2n-2$ coefficients and this one-step method can be seen as a numerical integrator from $(\q_0,\textbf{v}_0)$ to $(\q_1,\textbf{v}_1)$ with a $2n-1$ degree Hermite piecewise polynomial interpolating the configuration over every fixed time step. This approach can be extended to Lagrangian systems with forcing by discretizing the Lagrange-d'Alembert principle.

Galerkin Approach: One-step Galerkin methods with discretization using higher-order Hermite polynomials involve the use of shifted Legendre polynomials up to order $2n-2$ as test functions. Given $(\q_0,\textbf{v}_0)$, the governing discrete equations are given by

\begin{equation}
\int_0^{\dt} \textbf{M}_{\rm eq}(\q_{\rm d}(t),\dot{\q}_{\rm d}(t),\ddot{\q}_{\rm d}(t)) P_j(t)  \ \text{d}t = \textbf{0}, \quad \quad j=0,\cdots, 2n-3,
\end{equation}

where $P_j(t)$ are shifted Legendre polynomials of degree $j$.

§ NUMERICAL PROPERTIES

In order to understand the numerical properties of the proposed one-step methods, we investigate the linear stability and symplectic nature of these algorithms. We consider the simple harmonic oscillator with a single degree of freedom for both studies.

§.§ Linear Stability

We closely follow Leimkuhler and Reich [2] to study the linear stability of the proposed one-step methods for the simple harmonic oscillator. We consider the following Lagrangian system

\begin{equation}
L(q,\dot{q})=\frac{1}{2}\dot{q}^2 - \frac{1}{2}\omega^2q^2,
\end{equation}

where $\omega$ is the natural frequency. The governing equation is

\begin{equation}
\ddot{q} + \omega^2q^2 = 0.
\end{equation}

The analytical solution for this Lagrangian system is given by

\begin{equation}
\begin{bmatrix}
v(t) \\ \omega q(t) \end{bmatrix} =
\underbrace{ \begin{bmatrix}
\cos(\omega t) & -\sin(\omega t) \\
\sin(\omega t) & \cos(\omega t) \end{bmatrix}}_{=A_{\omega}}
\begin{bmatrix}
v(0) \\ \omega q(0) \end{bmatrix} .
\end{equation}

Since $A_{\omega}A_{\omega}^\top =\mathbb{I}$, $A_{\omega}$ is orthogonal and thus, the eigenvalues are $\lambda_{1,2}=e^{\pm i \omega t}$ with $|\lambda_{1,2}|=1$.
For this linear dynamical system, time-marching equations for both one-step methods can be written in the following form

\begin{equation}
\begin{bmatrix}
v_{k+1} \\ \omega q_{k+1}
\end{bmatrix} = {\Ra{A_{z} }} \begin{bmatrix}
v_{k} \\ \omega q_k
\end{bmatrix},
\end{equation}

where $A_{z}$ is the amplification matrix with $z:=\dt\omega$.  A sufficient condition for the asymptotic stability of a numerical method is that the eigenvalues of the amplification matrix $A_{z}$ are in the unit disk of the complex plane and are simple if they lie on the unit circle. We investigate the linear stability of both Hermite-based one-step methods and compare the results with the variational integrators based on the discrete mechanics.

* For the one-step variational method we have

\begin{equation*}
A_{z} = \begin{bmatrix}
\frac{7z^4 - 192z^2 + 420 }{2z^4 + 18z^2 + 420}&    \frac{15z(3z^2 - 28)}{2z^4 + 18z^2 + 420} \\   \frac{z(z^4 - 52z^2 + 420) }{2z^4 + 18z^2 + 420} &    \frac{3z^4 - 104z^2 + 240 }{2z^4 + 18z^2 + 420} \end{bmatrix},
\end{equation*}

with eigenvalues $\lambda_{1,2}=\frac{17z^4-192z^2 + 420 \pm z \sqrt{15(z^2-10)(3z^2-28)(z^2-42)}}{2z^4 + 18z^2 + 420}$. The stability region for the one-step variational method is shown in Figure <ref>. We observe a small region of instability for $\sqrt{\frac{28}{3}}< z< \sqrt{10}$. Thus, the one-step variational method is stable for $z<\sqrt{\frac{28}{3}}$.

* For the one-step Galerkin method we have

\begin{equation*}
A_{z} = \begin{bmatrix}
\frac{3z^4 - 104z^2 + 240 }{3z^4 + 48z^2 + 720}&    \frac{24z(z^2 - 10)}{3z^4 + 48z^2 + 720} \\   \frac{z(z^4 - 72z^2 + 720) }{3z^4 + 48z^2 + 720} &    \frac{3z^4 - 104z^2 + 240 }{3z^4 + 48z^2 + 720} \end{bmatrix},
\end{equation*}

with eigenvalues $\lambda_{1,2}=\frac{3z^4-104z^2 + 240 \pm 2z \sqrt{2(z^2-10)(z^2-12)(z^2-60)}}{z^4 + 16z^2 + 240}$. The stability region for the one-step Galerkin method is shown in Figure <ref>. Simillar to the one-step variational method, we observe a small region of instability for $\sqrt{10}< z< \sqrt{12}$. Thus, the one-step Galerkin method is stable for $z<\sqrt{10}$.

* The amplification matrix $A_{z}$ for the variational integrator based on the midpoint rule is

\begin{equation*}
A_{z} = \begin{bmatrix}
\frac{4-z^2 }{ z^2 + 4}&    \frac{4z}{z^2 + 4} \\   \frac{-4z }{z^2 + 4} &    \frac{4-z^2  }{z^2 + 4} \end{bmatrix},
\end{equation*}

The amplification matrix for this method is orthogonal and hence, the method is stable for all $z \in \mathbb{R}$.

The eigenvalues for the proposed one-step methods are plotted for different $z$ values in Figure <ref> where both one-step methods have similar stability regions. The one-step Galerkin method is stable for $\dt<\frac{\sqrt{10}}{\omega}$ whereas the one-step variational method is stable for $\dt < \frac{\sqrt{\frac{28}{3}}}{\omega}$.

[subfigure]oneside,margin=1.8cm,0 cm

One-step Galerkin

One-step variational

Linear stability analysis of both one-step methods.

§.§ Symplectic Nature

As mentioned in Section <ref>, the variational approach to one-step methods only considers the action integral over one time step whereas the variational integrators consider the action integral over a finite number of fixed time steps. We know that variational integrators derived from the latter approach yield numerical algorithms that automatically preserve the canonical symplectic form. For a Hamiltonian system with Hamiltonian $H(\textbf{p},\q)$, the  symplectic flow map $\phi_t(\textbf{p}_0,\q_0)=(\textbf{p}(t),\q(t))$ satisfies the following condition

\begin{equation}
\left( \frac{\partial \phi_t}{\partial \textbf{y}_0} \right)^\top J \left( \frac{\partial \phi_t}{\partial \textbf{y}_0} \right) =J,
\end{equation}

where $\textbf{y}_0=(\textbf{p}_0,\q_0)$ and $J=\begin{bmatrix}
0 & I \\ -I & 0
\end{bmatrix}$ is the symplectic matrix. Similar to this condition, a given one-step method $\phi_{\dt}:(\textbf{p}_k,\q_k) \to (\textbf{p}_{k+1},\q_{k+1})$ is symplectic if it satisfies $ \left( \frac{\partial \phi_{\dt}}{\partial \textbf{y}_k} \right)^\top J \left( \frac{\partial \phi_{\dt}}{\partial \textbf{y}_k} \right) =J$ for $\textbf{y}_k=(\textbf{p}_k,\q_k)$. The key step in this process is to compute the following Jacobian matrix

\begin{equation}
\frac{\partial \phi_{\dt}}{\partial \textbf{y}_k} = \begin{bmatrix}
\frac{\partial \textbf{p}_{k+1}(\textbf{p}_k,\q_k) }{\partial \textbf{p}_k}&  \frac{\partial \textbf{p}_{k+1}(\textbf{p}_k,\q_k) }{\partial \q_k} \\ \frac{\partial \q_{k+1}(\textbf{p}_k,\q_k) }{\partial \textbf{p}_k}&  \frac{\partial \q_{k+1}(\textbf{p}_k,\q_k) }{\partial \q_k}
\end{bmatrix}.
\end{equation}

Since both proposed methods are generally implicit, the computation for the Jacobian matrix involves differentiating the governing discrete equations and then solving a system of linear equations for the entries in the Jacobian matrix. We study the symplectic nature of the one-step algorithms for both linear and nonlinear conservative systems. It is important to note that the one-step methods developed in this paper are formulated on the state space. In order to check the condition for symplecticity, we need to define $\textbf{p}_k=m\textbf{v}_k$ to write these algorithms on phase space. Instead of writing the algorithms on phase space, we pick $m=1$ to simplify the expressions.
First, we check the condition for the simple harmonic oscillator with $L(q,\dot{q})=\frac{1}{2}\dot{q}^2 -\frac{1}{2}q^2$ where $q(t)$ is the displacement. For a fixed time step $\dt$, the Jacobians for one-step methods are

\begin{equation*}
\frac{\partial \phi_{\dt,V}}{\partial \textbf{y}_k} = \begin{bmatrix}
\frac{(44\dt^5 + 143\dt^4 - 700\dt^3 + 189\dt^2 + 1176\dt - 882) }{(7\dt(\dt^4 + 9\dt^2 + 210))}&  \frac{(8\dt^5 + 33\dt^4 - 224\dt^3 - 217\dt^2 + 1568\dt + 294) }{(7(\dt^4 + 9\dt^2 + 210))} \\ \frac{-(66\dt^5 + 169\dt^4 - 434\dt^3 + 1092\dt^2 + 588\dt + 1764) }{(14\dt(\dt^4 + 9\dt^2 + 210))}&  \frac{-(12\dt^5 + 39\dt^4 - 224\dt^3 - 56\dt^2 + 784\dt - 588) }{(14\dt(\dt^4 + 9\dt^2 + 210))} \end{bmatrix},
\end{equation*}

\begin{equation*}
\frac{\partial \phi_{\dt,G}}{\partial \textbf{y}_k} = \begin{bmatrix}
\frac{ (2( \dt^4 - 52\dt^2 + 120))}{\dt^4 + 16\dt^2 + 240}&  \frac{(\dt(\dt^4 - 132\dt^2 + 1440)) }{6(\dt^4 + 16\dt^2 + 240)} \\ \frac{ -(\dt(\dt^4 - 130\dt^2 + 1200 )) }{5(\dt^4 + 16\dt^2 + 240)}&  \frac{-(\dt^6 - 240\dt^4 + 6240\dt^2  - 14400) }{60(\dt^4 + 16\dt^2 + 240)} \end{bmatrix}.
\end{equation*}

For this linear  dynamical system, Jacobians from both variational and Galerkin one-step methods satisfy the condition for symplecticity. For a general mechanical system with a nonlinear potential energy $U(q)$ and Lagrangian $L=\frac{1}{2}m\dot{q}^2 -U(q)$, we find that none of the one-step methods satisfy the required condition for symplecticity. It is important to note that the above condition is only to check whether the algorithms preserve the canonical symplectic form. In fact, at present, one can only check whether a given integration scheme exhibits a specific symplectic structure; one cannot determine whether any such symplectic structure exists, in general. In the past, some of the well-known methods such as Newmark methods [10] have been shown to preserve a noncanonical symplectic form via nonlinear transformations but it is not generally known how to test for the existence of a noncanonical symplectic form for a given algorithm.

§ NUMERICAL RESULTS

In this section, we consider three examples to study the numerical performance of the proposed Galerkin and variational one-step methods. We first consider a nonlinear conservative system to demonstrate the good energy performance of the proposed methods. We also present an order analysis study to understand the convergence behavior of one-step variational and Galerkin methods. We then consider the Duffing oscillator to investigate the numerical behavior of the proposed one-step methods in the presence of dissipative forces. Finally, we consider a nonlinear aeroelastic system to study how the proposed methods perform for a coupled nonlinear dynamical system.

§.§ Particle in a Double-well Potential

In this subsection, we apply the proposed methods to a particle in a double-well potential with Lagrangian

\begin{equation}
L(q,\dot{q})= \frac{1}{2}m\dot{q}^2 - \frac{1}{2} \left( q^4 -q^2 \right),
\end{equation}

The Euler-Lagrange equation for this conservative system is given by

\begin{equation}
m\ddot{q} - q + 2q^3 = 0 .
\end{equation}

We have compared the numerical results for $m=1$ and fixed time step $\dt=0.1$ for two initial conditions. The phase portrait comparisons in Figure <ref> and Figure <ref>  show how both variational and Galerkin methods agree with the benchmark solution for both initial conditions. The corresponding energy error plots in Figure <ref> and Figure <ref> demonstrate the bounded energy error for the one-step methods. For both cases, the energy performance for the Galerkin method is substantially better than for the variational method. The energy error comparison in Figure <ref> shows that the one-step variational method has energy error magnitude around $10^{-7}$ whereas the Galerkin method has energy error around $10^{-10}$. Similarly, in Figure <ref> the variational method has energy error magnitude around $10^{-5}$ and the Galerkin method has energy error around $10^{-8}$.

[subfigure]oneside,margin=1.8cm,0 cm

Phase portrait

Energy error

Comparison between the two one-step methods for $q(0)=0.74, \ \dot{q}(0)=0$.

[subfigure]oneside,margin=1.8cm,0 cm

Phase portrait

Energy error

Comparison between the two one-step methods for $q(0)=0.995, \ \dot{q}(0)=0$.

In Figure <ref>, we have studied the numerical behavior of these algorithms for different fixed time step values to understand how the configuration, velocity and energy error values decrease with decrease in the step size. We have also performed the convergence analysis for variational integrators derived from the discrete mechanics framework to understand how they compare to the proposed one-step methods. Since both one-step methods proposed in this work require the solution of two coupled nonlinear equations at each time step, we have implemented variational integrators that require solving two nonlinear equations. Thus, we have considered a quadratic trajectory over one time step in the discrete mechanics framework by introducing an interior point [9] and the resulting variational integrator leads to two coupled nonlinear implicit equations at each time step.

The configuration and velocity convergence plots in Figure <ref> and Figure <ref> show that the one-step variational method has second order convergence whereas both the variational integrator and the one-step Galerkin method have fourth order convergence. For the energy error, the one-step variational and the variational integrator show second order convergence, whereas the one-step Galerkin method shows fourth order convergence. Thus, the one-step Galerkin approach gives better trajectory and energy performance than both the one-step variational method and the variational integrator.

[subfigure]oneside,margin=2cm,0 cm

Trajectory error

Velocity error

Energy error

Convergence analysis of maximum temporal error in configuration, velocity and energy for the nonlinear conservative system.

§.§ Duffing Oscillator

In this subsection, we study the numerical performance of these algorithms in the presence of dissipation. The governing second-order, nonlinear differential equation for the Duffing oscillator is

\begin{equation}
\ddot{x} + \delta \dot{x} + \alpha x + \beta x^3 = 0,
\end{equation}

where $x(t)$ is the displacement at time $t$, $\delta$ is the linear damping, $\alpha$ is the linear stiffness, and $\beta$ is the nonlinear stiffness coefficient. The Lagrangian and external forcing for this system are

\begin{equation}
L(x,\dot{x})=\frac{1}{2}\dot{x}^2  - \frac{1}{2}\alpha x^2 -\frac{1}{4}\beta x^4, \quad \quad f(\dot{x},t) = -\delta \dot{x} .
\end{equation}

We have fixed the stiffness parameters to $\alpha=1$ and $\beta=0.5$ and studied this dissipative nonlinear dynamical system for three cases, i.e. $\delta \in \{0.025,0.05,0.1\}$.  With these specific parameter values, the Duffing oscillator can be thought of as the double-well potential system with dissipation. The numerically computed trajectories from the one-step methods are compared with the benchmark solution in Figure <ref>. The plots in Figure <ref> demonstrate that both methods are able to capture the dissipation effect accurately. In fact, the discrete trajectories are indistinguishable from the benchmark solution for all three cases.

The energy error plots in Figure <ref> compare the energy performance for the one-step methods and the Galerkin approach outperforms the variational approach in all three cases. For all three cases, the energy error for the one-step variational method starts around $10^{-4}$ whereas the one-step Galerkin methods exhibit energy error around $10^{-6}$. The energy error for both methods decreases over the time due to the presence of dissipative forces and the rate of decrease in energy error increases with increasing values of the damping parameter $\delta$. This decrease in energy error is seen clearly in Figure <ref> and Figure <ref> for higher $\delta$ values.

[subfigure]oneside,margin=2cm,0 cm

Duffing oscillator numerical simulation with $\alpha=-1,\beta=2$, and time step $\dt=0.1$.

[subfigure]oneside,margin=2cm,0 cm

Energy error behavior for Duffing oscillator numerical simulation with $\alpha=-1,\beta=2$, and time step $\dt=0.1$.

§.§ Aeroelastic System

In this subsection, we consider the open-loop behavior of the nonlinear aeroelastic system studied by Shukla and Patil [24]. As shown in Figure <ref>, the model contains a flat plate supported by a linear spring in the plunge degree of freedom and the cubic nonlinear spring in the pitch degree of freedom. The flat plate is free to move up and down along the plunge degree of freedom and to rotate about the pitch degree of freedom. The Lagrangian for this system is given by

\begin{equation}
L(h, \dot{h},\alpha,\dot{\alpha})= m_{\rm T} \dot{h}^2 + I_{\alpha}\dot{\alpha}^2 + m_{\rm W}x_{\alpha}\dot{h}\dot{\alpha} - \frac{1}{2}k_hh^2 -\frac{1}{2}k_{\alpha_0}\alpha^2 -  \frac{1}{3}k_{\alpha_1}\alpha^3 - \frac{1}{4}k_{\alpha_2}\alpha^4  ,
\end{equation}

where $m_{\rm W}$ is the mass of the wing and $m_{\rm T}$ is the total mass of the aeroelastic system. The parameter $I_{\alpha}$ represents the moment of inertia about the elastic axis. The terms $k_h$ and $k_{\alpha_{\{0,1,2\}}}$ are the stiffness functions along the plunge and pitch degrees of freedom respectively.
The external nonconservative forces are

\begin{equation}
f_{h}= - c_h \dot{h} + \rho U^2bC_{L_{\alpha}} \alpha_{\rm eff},
\end{equation}

\begin{equation}
f_{\alpha} = -c_{\alpha}\dot{\alpha} + \rho U^2b^2C_{M_{\alpha}} \alpha_{\rm eff},
\end{equation}

where $c_h$ and $c_{\alpha}$ are damping coefficients, $C_{L_{\alpha}}$ and $C_{M_{\alpha}}$ are the derivatives of the lift and moment coefficients, and  $\alpha_{\rm eff}=\left( \alpha + \frac{\dot{h}}{U} + \left(\frac{1}{2}-a \right)b \frac{\dot{\alpha}}{U} \right)$ is the effective angle of attack. The equations of motion for this aeroelastic system are

\begin{equation}
m_{\rm T}\ddot{h} + m_{\rm W}x_{\alpha}b\ddot{\alpha} + c_h \dot{h} + k_hh -  \rho U^2bC_{L_{\alpha}} \alpha_{\rm eff} = 0,
\end{equation}

\begin{equation}
I_{\alpha}\ddot{\alpha} + m_{\rm W}x_{\alpha}b\ddot{h} + c_{\alpha} \dot{\alpha} + k_{\alpha}(\alpha)\alpha + \rho U^2b^2C_{M_{\alpha}} \alpha_{\rm eff}   = 0,
\end{equation}

where $k(\alpha)=k_{\alpha_0} + k_{\alpha_1}\alpha + k_{\alpha_2}\alpha^2$.

Sketch of a two degrees of freedom aeroelastic section model.

We have studied the dynamic behavior of the nonlinear aeroelastic system for initial conditions $(h_0,\dot{h}_0,\alpha_0,\dot{\alpha}_0)=(0.01,0,0.1,0)$ at freestream velocity $U=0.9U_f$ where $U_f$ is the linear flutter velocity. The phase space plots given in Figure <ref> and Figure <ref> clearly demonstrate how both one-step methods capture the subcritical limit cycle oscillations (LCOs) accurately.

[subfigure]oneside,margin=2cm,0 cm

Plunge response

Pitch response

Subcritical LCO simulation using the one-step variational method.

[subfigure]oneside,margin=2cm,0 cm

Plunge response

Pitch response

Subcritical LCO simulation using the one-step Galerkin method.

The energy plot in Figure <ref> shows how the total energy of the aeroelastic system evolves over time. Initially there is a sharp decrease in energy followed by an increase, and eventually when the system exhibits periodic motion with constant amplitude the total energy oscillates around a fixed value. As shown in Figure <ref>, both one-step methods track the change in energy accurately for the nonlinear aeroelastic system. The energy error comparison in Figure <ref> demonstrates how the Galerkin approach has better energy behavior than the one-step variational approach. The one-step variational method has energy error magnitude around $10^{-4}$ whereas the Galerkin method has energy error around $10^{-6}$.

Total energy of the aeroelastic system.

Energy error comparison between one-step Galerkin and variational methods.

§ CONCLUSIONS

In this paper we have developed Hermite polynomial based one-step variational and Galerkin methods for mechanical systems with external forcing. We have utilized cubic Hermite polynomials over one time step for discretization and the resulting numerical algorithms are continuous in both configuration and velocity. We also demonstrated an approach to obtain one-step methods using higher-order Hermite polynomials. We showed that both one-step methods are symplectic for linear dynamical systems but they do not preserve the canonical symplectic form for general nonlinear dynamical systems. We also investigated the linear stability of the proposed one-step methods and both methods exhibit excellent stability for large time steps.

We have studied the numerical behavior of these algorithms through three different numerical examples. The energy performance and convergence analysis results for the conservative example showed how both one-step methods achieve good numerical performance by obtaining $\mathcal{C}^1$–continuous trajectories. We have also presented results for a dissipative system and the numerical plots show that both one-step methods capture the effect of the dissipative forces accurately over long time intervals. Finally, numerical studies for the coupled aeroelastic system show how both one-step methods capture the limit cycle oscillations accurately. The numerical results from all three examples showed that the Galerkin approach has significantly better energy behavior. In fact, the one-step Galerkin method is also better than the variational integrator, with same nonlinear equations per time step, in terms of the energy accuracy.

Future research directions motivated by this work are: obtaining theoretical results about the convergence behavior and geometric properties of the proposed one-step methods; investigating the connection between the one-step Galerkin methods and energy-momentum integrators; and applying these one-step methods to discretizations of infinite-dimensional systems.

Declaration of interest: None.

Funding: This material is based upon work supported by the National Science Foundation under Grant No. 1826152.

[1]
Hairer E, Lubich C, Wanner G.
Geometric Numerical Integration: Structure-preserving Algorithms for
Ordinary Differential Equations. vol. 31.
Springer Science & Business Media; 2006.

[2]
Leimkuhler B, Reich S.
Simulating Hamiltonian Dynamics.
Cambridge Monographs on Applied and Computational Mathematics.
Cambridge University Press; 2005.

[3]
Wendlandt J, Marsden J.
Mechanical integrators derived from a discrete variational principle.
Physica D: Nonlinear Phenomena. 1997;106(3-4):223–246.

[4]
Maeda S.
Extension of discrete Noether theorem.
Math Japonica. 1981;26(1):85–90.

[5]
Maeda S.
Lagrangian formulation of discrete systems and concept of
difference space.
Math Japonica. 1982;27:345–356.

[6]
Moser J, Veselov AP.
Discrete versions of some classical integrable systems and
factorization of matrix polynomials.
Communications in Mathematical Physics. 1991;139(2):217–243.

[7]
Veselov AP.
Integrable Lagrangian correspondences and the factorization of
matrix polynomials.
Functional Analysis and Its Applications. 1991;25(2):112–122.

[8]
Veselov AP.
Integrable discrete-time systems and difference operators.
Functional Analysis and its Applications. 1988;22(2):83–93.

[9]
Marsden J, West M.
Discrete mechanics and variational integrators.
Acta Numerica. 2001;10:357–514.

[10]
Kane C, Marsden JE, Ortiz M, West M.
Variational integrators and the Newmark algorithm for conservative
and dissipative mechanical systems.
International Journal for Numerical Methods in Engineering.

[11]
Leyendecker S, Marsden JE, Ortiz M.
Variational integrators for constrained dynamical systems.
ZAMM-Journal of Applied Mathematics and Mechanics.

[12]
Leok M, Shingel T.
Prolongation–collocation variational integrators.
IMA Journal of Numerical Analysis. 2012;32(3):1194–1216.

[13]
Kane C, Marsden JE, Ortiz M.
Symplectic-energy-momentum preserving variational integrators.
Journal of mathematical physics. 1999;40(7):3353–3371.

[14]
Sharma H, Patil M, Woolsey C.
Energy-preserving variational integrators for forced Lagrangian
Communications in Nonlinear Science and Numerical Simulation.

[15]
Sharma H, Patil M, Woolsey C.
A review of structure-preserving numerical methods for engineering
Computer Methods in Applied Mechanics and Engineering.

[16]
Argyris J, Scharpf D.
Finite elements in time and space.
The Aeronautical Journal. 1969;73(708):1041–1044.

[17]
Baruch M, Riff R.
Hamilton's principle, Hamilton's law-6 to the n power correct
AIAA Journal. 1982;20(5):687–692.

[18]
Riff R, Baruch M.
Time finite element discretization of Hamilton's law of varying
AIAA journal. 1984;22(9):1310–1318.

[19]
Mergel JC, Sauer RA, Ober-Blöbaum S.
C1-continuous space-time discretization based on Hamilton's law of
varying action.
ZAMM-Journal of Applied Mathematics and Mechanics/Zeitschrift für
Angewandte Mathematik und Mechanik. 2017;97(4):433–457.

[20]
Hulme BL.
One-step piecewise polynomial Galerkin methods for initial value
Mathematics of Computation. 1972;26(118):415–426.

[21]
Betsch P, Steinmann P.
Conservation properties of a time FE method. Part I: time-stepping
schemes for N-body problems.
International Journal for Numerical Methods in Engineering.

[22]
Betsch P, Steinmann P.
Conservation properties of a time FE method—part II:
Time-stepping schemes for non-linear elastodynamics.
International Journal for Numerical Methods in Engineering.

[23]
Riff R, Baruch M.
Stability of time finite elements.
AIAA journal. 1984;22(8):1171–1173.

[24]
Shukla H, Patil M.
Nonlinear state feedback control design to eliminate subcritical
limit cycle oscillations in aeroelastic systems.
Nonlinear Dynamics. 2017;88(3):1599–1614.