Functionally-fitted energy-preserving methods for solving oscillatory nonlinear Hamiltonian systems

In this paper, we consider nonlinear Hamiltonian systems:

where $d=2d_{1},f:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d},H:\mathbb{R}^{d}\rightarrow\mathbb{R}$ are sufficiently smooth functions and

is the canonical symplectic matrix. It is well known that the flow of (1) preserves the symplectic form $dy\wedge Jdy$ and the Hamiltonian or energy $H(y(t))$. In the spirit of geometric numerical integration, it is a natural idea to design schemes that preserve both the symplecticity of the flow and Hamiltonian function. Unfortunately, a numerical scheme cannot achieve this goal unless it generates the exact solution (see, e.g. [21], page 379). Hence researchers face a choice between preserving symplecticity or energy and many of them have given more weight on the former in the last decades, and readers are referred to [21] and references therein. Whereas investigations on energy-preserving (EP) methods are relatively insufficient (see, e.g. [1, 4, 6, 10, 11, 16, 18, 22, 26, 35]). Comparing to symplectic methods, EP methods are beneficial for improving nonlinear stability, easier to adapt the time step and more suitable for the integration of chaotic systems (see, e.g. [12, 19, 32, 34]).

On the other hand, in scientific computing and modelling, the design and analysis of methods for periodic or oscillatory systems have been considered by many authors (see, e.g. [2, 17, 20, 31, 38, 41]). Generally, these methods utilize a priori information of special problems and they are more efficient than general-purpose methods. A popular approach to constructing methods suitable for oscillatory problems is using the functionally-fitted (FF) condition, namely, deriving a suitable method by requiring it to integrate members of a given finite-dimensional function space $X$ exactly. If $X$ incorporates trigonometrical or exponential functions, the corresponding methods are also named by trigonometrically-fitted (TF) or exponentially-fitted (EF) methods (see, e.g. [14, 25, 29, 33]).

Therefore, combining the ideas of the EF/TF and structure-preserving methods is a promising approach to developing numerical methods which allow long-term computation of solutions to oscillatory Hamiltonian systems (1). Just as the research of symplectic and EP methods, EF/TF symplectic methods have been studied extensively by many authors (see, e.g. [7, 8, 9, 15, 36, 37, 40]). By contrast, as far as we know, only a few papers paid attention to the EF/TF EP methods (see, e.g. [27, 28, 39]). Usually the existing EF/TF EP methods are derived in the context of continuous-stage Runge–Kutta (RK) methods. The coefficients in these methods are determined by a system of equations resulting from EF/TF, EP and symmetry conditions. As mentioned at the end of [27], it is not easy to find such a system having a unique solution in the case of deriving high-order methods. What’s more, how to verify the algebraic order of such methods falls into a question. A common way is to check order conditions related to rooted trees. Again, this is inconvenient in the high-order setting since the number of trees increases extremely fast as the order grows. In this paper, we will construct FF EP methods based on the continuous finite element method, which is inherently energy-preserving (see, e.g. [1, 16, 35]). Intuitively, we are expected to increase the order of the method through enlarging the finite element space. By adding trigonometrical functions to the space, the corresponding method is naturally trigonometrically fitted. Thus we are hopeful of constructing FF EP methods, in particular TF EP methods, of arbitrarily high order.

The outline of this paper is as follows. In Section 2, we construct FF CFE methods and present important geometric properties of them. In Section 3, we interpret them as continuous-stage Runge–Kutta methods and analyse the algebraic order. We then discuss implementation details of these new methods in Section 4. Numerical results are shown in section 5, including the comparison between our new TF EP methods and other prominent structure-preserving methods in the literature. The last section is concerned with the conclusion and discussion.

Given two finite-dimensional function spaces $X$ and $Y$ whose members are $\mathbb{R}^{d}$-valued, the continuous finite element method for (1) is described as follows.

Find a continuous piecewise X-type function $U(t)$ on $I$ with $U(t_{0})=y_{0}$, such that

for any piecewise Y-type function $v(t)$, where $U(t)\approx y(t)$ on $I$ and $y(t)$ solves the IVP (1). The term ‘continuous finite element’(CFE) comes from the continuity of the finite element solution $U(t)$. Since (2) deals with an initial value problem, we need only to consider it on $[t_{0},t_{0}+h]$.

Find $u\in X_{h}(t_{0})$ with $u(0)=y_{0}$, such that

for any $v\in Y_{h}(t_{0})$, where

for $\tau\in[0,1]$. Since $U(t)$ is continuous, $y_{1}=u(1)$ is the initial value of the local problem on the next interval $[t_{1},t_{2}]$. Thus we can solve the global variational problem (2) on $I$ step by step.

In the special case of

(2) reduces to the classical continuous finite element method (see, e.g. [1, 24]) denoted by CFE$r$ in this paper. For the purpose of deriving functionally-fitted methods, we generalize $X$ and $Y$ a little:

Thus it is sufficient to give $X$ or $Y$ since they can be determined by each other. Besides, $Y$ is supposed to be invariant under translation and reflection , namely,

Clearly, $Y_{h}(t_{0})$ and $X_{h}(t_{0})$ are irrelevant to $t_{0}$ provided (5) holds. For convenience, we simplify $Y_{h}(t_{0})$ and $X_{h}(t_{0})$ by $Y_{h}$ and $X_{h}$ respectively. In the remainder of this paper, we denote the CFE method (2) or (3) based on the general function spaces (4) satisfying the condition (5) by FFCFE$r$.

It is noted that the FFCFE$r$ method (3) is defined by a variational problem and the well-definedness of this problem has not been confirmed yet. Here we presume the existence and uniqueness of the solution to (3). This assumption will be proved in the next section. With this premise, we are able to present three significant properties of the FFCFE$r$ method. At first, by the definition of the variational problem (2), the FFCFE$r$ method is functionally fitted with respect to the space $X$.

Moreover, the FFCFE$r$ method is inherently energy preserving method. The next theorem confirms this point.

∎

The FFCFE$r$ method can also be viewed as a one-step method $\Phi_{h}:y_{0}\rightarrow y_{1}=u(1)$. It is well known that reversible methods show a better long-term behaviour than nonsymmetric ones when applied to reversible differential systems such as (1) (see, e.g. [21]). This fact motivates the investigation of the symmetry of the FFCFE$r$ method. Since spaces $X$ and $Y$ satisfy the invariance (5), which is a kind of symmetry, the FFCFE$r$ method is expected to be symmetric.

∎

It is well known that polynomials cannot approximate oscillatory functions satisfactorily. If the problem (1) has a fixed frequency $\omega$ which can be evaluated effectively in advance, then the function space containing the pair $\left\{\sin(\omega t),\cos(\omega t)\right\}$ seems to be a more promising candidate for $X$ and $Y$ than the polynomial space. For example, possible $Y$ spaces for deriving the TF CFE method are

and

Correspondingly, by equipping the FF CFE method with the space $Y=Y_{1},Y_{2}$ or $Y_{3}$, we obtain three families of TF CFE methods. According to Theorem 2.3 and Theorem 2.4, all of them are symmetric energy-preserving methods. To exemplify this framework of the TF CFE method, in numerical experiments, we will test the TF CFE method denoted by TFCFE$r$ and TF2CFE$r$ based on the spaces (7) and (8). It is noted that TFCFE2 and TF2CFE2 coincide.

An interesting connection between CFE methods and RK-type methods has been shown in a few papers (see, e.g. [3, 24, 35]). Since the RK methods are dominant in the numerical integration of ODEs, it is meaningful and useful to transform the FFCFE$r$ method to the corresponding RK-type method which has been widely and conventionally used in applications. After the transformation, the FFCFE$r$ method can be analysed and implemented by standard techniques in ODEs conveniently. To this end, it is helpful to introduce the projection operation $P_{h}$. Given a continuous $\mathbb{R}^{d}$-valued function $w$ on $[0,1]$, its projection onto $Y_{h}$, denoted by $P_{h}w$, is defined by

Clearly, $P_{h}w(\tau)$ can be uniquely expressed as a linear combination of $\{\tilde{\varphi}_{i}(\tau)\}_{i=0}^{r-1}$:

Taking $v(\tau)=\tilde{\varphi}_{i}(\tau)e_{j}$ in (10) for $i=0,1,\ldots,r-1$ and $j=1,\ldots,d$, it can be observed that coefficients $U_{i}$ satisfy the equation

where

Since $\left\{\tilde{\varphi}_{i}\right\}_{i=0}^{r-1}$ are linearly independent for $h>0$, the stiffness matrix $M$ is nonsingular. Thus the projection can be explicitly expressed by

where

Clearly, $P_{\tau,\sigma}$ can be calculated by a basis other than $\left\{\tilde{\varphi}_{i}\right\}_{i=0}^{r-1}$ since they only differ in a linear transformation. If $\left\{\phi_{0},\ldots,\phi_{r-1}\right\}$ is an orthonormal basis of $X_{h}$ under the inner product $\langle\cdot,\cdot\rangle$, then $P_{\tau,\sigma}$ admits a simpler expression:

Now using (3) and the definition (10) of the operator $P_{h}$, we obtain that $u^{\prime}=hP_{h}(f\circ u)$ and

Integrating the above equation with respect to $\tau$, we transform the FFCFE$r$ method (3) into the continuous-stage RK method:

where

In particular,

for the CFE$r$ method for $i=0,1,\ldots,r-1$, where $\hat{p}_{i}(\tau)$ is the shifted Legendre polynomial of degree $i$ on $[0,1]$, scaled in order to be orthonormal. Thus the CFE$r$ method in the form (14) is identical to the energy-preserving collocation method of order $2r$ (see [22]) or the Hamiltonian boundary value method HBVM$(\infty,r)$ (see, e.g. [4]). For the FFCFE$r$ method, since $P_{\tau,\sigma},A_{\tau,\sigma}$ are functions of variable $h$ and $u(\tau)$ is implicitly determined by (14), it is necessary to analyse their smoothness with respect to $h$ before investigating the analytic property of the numerical solution $u(\tau)$. First of all, it can be observed from (11) that $P_{\tau,\sigma}=P_{\tau,\sigma}(h)$ is not defined at $h=0$ since the matrix $M$ is singular in this case. Fortunately, the following lemma shows that the singularity is removable.

∎

From (16) and (19), it can be observed that the FFCFE$r$ method (14) reduces to the CFE$r$ method when $h\to 0$, or equivalently, the energy-preserving collocation method of order $2r$ and HBVM$(\infty,r)$ method mentioned above. Since $A_{\tau,\sigma}=\int_{0}^{\tau}P_{\alpha,\sigma}d\alpha$ is also a smooth function of $h$, we can assume that

Furthermore, since the right function $f$ in (1) maps from $\mathbb{R}^{d}$ to $\mathbb{R}^{d}$, the $n$th-order derivative of $f$ at $y$ denoted by $f^{(n)}(y)$ is a multilinear map from $\underbrace{\mathbb{R}^{d}\times\ldots\times\mathbb{R}^{d}}_{n-fold}$ to $\mathbb{R}^{d}$ defined by

where $y=(y_{1},\ldots,y_{d})^{\intercal}$ and $z_{i}=(z_{i}^{1},\ldots,z_{i}^{d})^{\intercal}$, $i=1,\ldots,n$. With this background, we now can give the existence, uniqueness, especially the smoothness with respect to $h$ for the continuous finite element approximation $u(\tau)$ associated with the FFCFE$r$ method. The proof of the following theorem is based on a fixed-point iteration which is analogous to Picard iteration.

∎

Since our analysis on the algebraic order of the FFCFE$r$ method is mainly based on Taylor’s theorem, it is meaningful to investigate the expansion of $P_{\tau,\sigma}$.