High-order generalized-$\alpha$ methods

TL;DR

This paper extends the generalized-$ extalpha$ method to third-order accuracy, maintaining stability and controllable numerical dissipation, thus enhancing time integration techniques for computational simulations.

Contribution

The paper introduces a third-order generalized-$ extalpha$ method, proving its unconditional stability and discussing potential for higher-order generalizations.

Findings

Third-order method is unconditionally stable.

Numerical dissipation can be controlled similarly to second-order.

Easy integration into existing second-order generalized-$ extalpha$ frameworks.

Abstract

The generalized-$\alpha$ method encompasses a wide range of time integrators. The method possesses high-frequency dissipation while minimizing unwanted low-frequency dissipation and the numerical dissipation can be controlled by the user. The method is unconditionally stable and is of second-order accuracy in time. We extend the second-order generalized-$\alpha$ method to third-order in time while the numerical dissipation can be controlled in a similar fashion. We establish that the third-order method is unconditionally stable. We discuss a possible path to the generalization to higher order schemes. All these high-order schemes can be easily implemented into programs that already contain the second-order generalized-$\alpha$ method.

Tables

Table 1. Table 1: C ​ ( p ) 𝐶 𝑝 C(p) for various values of p 𝑝 p .

$k$	2	3	4	5	6	7	8	9	10	11
$C\left(p\right)$	$\frac{1}{2}$	$\frac{5}{12}$	$\frac{1}{3}$	$\frac{31}{120}$	$\frac{1}{5}$	$\frac{41}{252}$	$\frac{1}{7}$	$\frac{31}{240}$	$\frac{1}{9}$	$\frac{61}{660}$

Equations

$$\dot{u}+\lambda u$$

$$u(0)$$

$$V_{n+{\alpha }_m}$$

$$U_{n+1}$$

$$V_{n+{\alpha }_m}$$

$$U_{n+{\alpha }_f}$$

$$V_n^{{\alpha }_m}$$

$$V_{n+1}$$

$$U_{n+1}$$

$$V_n^{{\alpha }_m}$$

$$U_n^{{\alpha }_f}$$

$$U_{n}^{\alpha_{f}}=U_{n}+\tau\big{(}V_{n}+\alpha_{f}(V_{n+1}-V_{n})\big{)}=U_{n}+\tau V_{n+\alpha_{f}}.$$

$$V_n^{{\alpha }_m}$$

$$U_{n+1}^{(k)}$$

$$⋮$$

$$U_{n+1}^{(1)}$$

$$U_{n+1}^{(0)}$$

$$V_n^{{\alpha }_m}$$

$$U_n^{{\alpha }_f}$$

$${\gamma }_1={\gamma }_2=\frac{5}{12}+{\alpha }_m-{\alpha }_f.$$

$$U_{n+1}-\frac{{\tau }^2}{2}{\gamma }_2{A}_{n+1}$$

$$V_{n+1}-\tau {\gamma }_1{A}_{n+1}$$

$$\tau {\alpha }_f\lambda V_{n+1}+\tau {\alpha }_m{A}_{n+1}$$

$$\left[\begin{array}{ccc}1& 0& -\frac{{\gamma }_2}{2}\\ 0& 1& -{\gamma }_1\\ 0& {\alpha }_f\lambda \tau & {\alpha }_m\end{array}\right]\left[\begin{array}{c}{U}_{n+1}\\ \tau V_{n+1}\\ {\tau }^2{A}_{n+1}\end{array}\right]=\left[\begin{array}{ccc}1& 1& \frac{1-{\gamma }_2}{2}\\ 0& 1& 1-{\gamma }_1\\ -\lambda \tau & ({\alpha }_f-1)\lambda \tau -1& ({\alpha }_m-1)\end{array}\right]\left[\begin{array}{c}{U}_n\\ \tau V_n\\ {\tau }^2{A}_n\end{array}\right].$$

$$G$$

$$=\left[\begin{array}{ccc}\frac{2{\alpha }_m-({\gamma }_2-2{\gamma }_1{\alpha }_f)\lambda \tau }{2{\alpha }_m+2{\gamma }_1{\alpha }_f\lambda \tau }& \frac{2{\alpha }_m+2{\gamma }_1{\alpha }_f\lambda \tau -{\gamma }_2(1+\lambda \tau )}{2{\alpha }_m+2{\gamma }_1{\alpha }_f\lambda \tau }& \frac{{\alpha }_m+{\gamma }_1{\alpha }_f\lambda \tau -{\gamma }_2(1+{\alpha }_f\lambda \tau )}{2({\alpha }_m+{\gamma }_1{\alpha }_f\lambda \tau )}\\ -\frac{{\gamma }_1\lambda \tau }{{\alpha }_m+{\gamma }_1{\alpha }_f\lambda \tau }& \frac{{\alpha }_m-{\gamma }_1+{\gamma }_1({\alpha }_f-1)\lambda \tau }{{\alpha }_m+{\gamma }_1{\alpha }_f\lambda \tau }& \frac{-{\gamma }_1+{\alpha }_m}{{\alpha }_m+{\gamma }_1{\alpha }_f\lambda \tau }\\ -\frac{\lambda \tau }{{\alpha }_m+{\gamma }_1{\alpha }_f\lambda \tau }& -\frac{1+\lambda \tau }{{\alpha }_m+{\gamma }_1{\alpha }_f\lambda \tau }& \frac{-1+{\alpha }_m+(-1+{\gamma }_1){\alpha }_f\lambda \tau }{{\alpha }_m+{\gamma }_1{\alpha }_f\lambda \tau }\end{array}\right].$$

$$G_0{U}_{n+1}-G_1{U}_n+G_2{U}_{n-1}-G_3{U}_{n-2}=0,$$

$$U_{n+1}$$

$$U_{n-1}$$

$$U_{n-2}$$

$$(G_0-G_1+G_2-G_3)U_n+\tau (G_0-G_2+2G_3)V_n+\frac{{\tau }^2}{2}(G_0+G_2-4G_3)A_n+\frac{{\tau }^3}{6}(G_0-G_2+8G_3){\dot{A}}_n=\mathcal{O}({\tau }^4),$$

$$V_n$$

$$A_n$$

$${\dot{A}}_n$$

$$\displaystyle\frac{\lambda^{3}\tau^{3}}{12(\alpha_{m}+\gamma_{1}\alpha_{f}\lambda\tau)}\cdot\Big{[}$$

$$\displaystyle+\lambda\tau(-5-2\gamma_{1}+6\gamma_{2}+12\alpha_{f}-12\gamma_{1}\alpha_{f})\Big{]}=\mathcal{O}(\tau^{4})$$

$$-5+6{\gamma }_1+6{\gamma }_2+12{\alpha }_f-12{\alpha }_m$$

$${\gamma }_2-{\gamma }_1$$

$${\gamma }_1={\gamma }_2=\frac{5}{12}+{\alpha }_m-{\alpha }_f$$

$${\gamma }_1=\frac{3{\alpha }_m}{2+3{\alpha }_f},{\gamma }_2=\frac{10-9{\alpha }_f-36{\alpha }_f^2+6{\alpha }_m+36{\alpha }_m{\alpha }_f}{12+18{\alpha }_f},$$

$$-5+6{\gamma }_1+6{\gamma }_2+12{\alpha }_f-12{\alpha }_m$$

$$-5-2{\gamma }_1+6{\gamma }_2+12{\alpha }_f-12{\gamma }_1{\alpha }_f$$

Full text

††journal: arXiv

1 Introduction

The generalized-$\alpha$ method was introduced by Chung and Hulbert in (4) for solving hyperbolic equations arising in structural dynamics. The method was then applied to solve the parabolic equations such as the Reynolds-averaged Navier-Stokes equations in the computational fluid dynamics; see (7). Since then, the method has been widely used in engineering and sciences due to the following three attractive features: second-order accuracy in time, unconditional stability, and user-control on the high-frequency numerical dissipation.

The generalized-$\alpha$ method produces an algorithm which provides an optimal combination of high-frequency and low-frequency dissipation in the sense that for a given value of high-frequency dissipation, the algorithm minimizes the low-frequency dissipation; see (4). The robustness of the generalized-$\alpha$ method for parabolic systems has been successfully used to simulate a wide range of engineering applications (8, 1, 2, 5, 6).

To our best knowledge, the generalized-$\alpha$ method are limited to be second-order accurate in time while the high-order Runge-Kutta schemes, the Adams-Moulton schemes, and backward differentiation formulas (see (3)) do not control numerical dissipation explicitly. Thus, the goal of this work is to devise and analyze a third-order generalized-$\alpha$ method. The generalized-$\alpha$ method involves a parameter $\alpha_{m}$ for the time-derivative term and a parameter $\alpha_{f}$ for other terms. These parameters are then used for the discrete representation of the ordinary differential equation. The main idea of our generalization is to view the solution representations as Taylor expansions and then assign extra parameters to the higher-order terms in this expansion to achieve higher-order accuracy. More precisely, in the generalized-$\alpha$ method, the term $U_{n+\alpha_{f}}$ in the discrete PDE is written as $U_{n}+\alpha_{f}(U_{n+1}-U_{n})$. This representation limits the method from obtaining higher-order accuracy. Thus, we view $U_{n}+\alpha_{f}(U_{n+1}-U_{n})$ as a Taylor expansion of $U_{n+\alpha_{f}}$ and expand the representation $U_{n}^{\alpha_{f}}=U_{n}+\tau\dot{U}_{n}+\tau\alpha_{f}(\dot{U}_{n+1}-\dot{U}_{n})$ to seek for third-order accuracy. Herein, $\tau$ is the time-step size and the dot specifies a time derivative. The parameter $\alpha_{f}$ in the notation $U_{n+\alpha_{f}}$ in the generalized-$\alpha$ method behaves like a sub-time-step. This notation $U_{n+\alpha_{f}}$ is no longer valid for $U_{n}+\tau\dot{U}_{n}+\tau\alpha_{f}(\dot{U}_{n+1}-\dot{U}_{n})$. Thus, we use $U_{n}^{\alpha_{f}}$. To a certain extent, the sub-time-step is performed on the higher-order derivatives. We represent the other terms in the discrete PDE using a similar construction. We determine the free parameters $\alpha_{m}$ and $\alpha_{f}$ by Taylor series analysis. We then study the spectral properties of the resulting amplification matrix to determine the unconditional stability region. The numerical dissipation is then user-controlled by the values of $\alpha_{m}$ and $\alpha_{f}$ in the unconditional stability region.

The rest of this paper is organized as follows. Section 2 describes the main idea of the high-order generalized-$\alpha$ methods. We prove the formal third-order accuracy in time. Section 3 establishes the unconditional stability region. We also discuss the control on the eigenvalues of the resulting amplification matrix. Concluding remarks are given in Section 4.

2 High-order generalized-$\alpha$ methods

We consider the first-order ordinary differential equation (ODE)

[TABLE]

where $u_{0}$ is the initial solution. We partition the time interval $[0,\mathcal{T}]$ as $0=t_{0}<t_{1}<\cdots<t_{N}=\mathcal{T}$ and let $\tau_{n}=t_{n+1}-t_{n}$ be the time step-size. We assume a uniform partitioning and denote the time step-size as $\tau$. We denote by $U_{n},V_{n},A_{n}$ the approximations to $U(t_{n}),\dot{U}(t_{n}),\ddot{U}(t_{n})$, respectively. The time-marching scheme of the generalized-$\alpha$ method for solving (2.1) is given by:

[TABLE]

with the initial solution $U_{0}=u_{0}$ and initial velocity $-\lambda U_{0}$.

The equations in (2.2c) and (2.2d) represent the approximations of $V_{n+\alpha_{m}}$ and $U_{n+\alpha_{f}}$, respectively, while the equation (2.2a) concerns the consistency of the discrete approximation to the ODE. When $\gamma_{1}=\frac{1}{2}+\alpha_{m}-\alpha_{f}$, this is equivalent to the second-order generalized-$\alpha$ method written in a different way; see (7).

2.1 Main idea

Equations in (2.2c)-(2.2d) resemble a sub-step time-marching. The limited accuracy of the sub-step time-marching restricts the method from obtaining higher-order accuracy. The reason is that the linear combination of two Taylor expansions can remove part of the error terms.

To overcome this limitation, we view these representations as low-order accurate Taylor expansions. We then seek higher-order schemes by applying higher-order accurate Taylor expansions. Thus, for example, we seek a third-order generalized-$\alpha$ method in the form

[TABLE]

with the initial solution $U_{0}=u_{0}$, initial velocity $-\lambda U_{0}$, and initial acceleration $\lambda^{2}U_{0}$.

We determine the coefficients $\gamma_{1}$ and $\gamma_{2}$ to achieve third-order accuracy. Herein, we change the notation $V_{n+\alpha_{m}}$ and $U_{n+\alpha_{f}}$ to $V_{n}^{\alpha_{m}}$ and $U_{n}^{\alpha_{f}}$ as the terms on the right-hand sides of the equations (2.3d) and (2.3e) are no longer sub-step time-marching on the primary variables. They are approximations in terms of the parameters $\alpha_{m}$ and $\alpha_{f}$. The sub-step time-marching occurs on the first derivative. More precisely, we can rewrite

[TABLE]

Thus, for $U_{n}^{\alpha_{f}}$, the sub-step time-marching is on its derivative $V_{n}$. By considering sub-step time-marching on higher-order derivatives, we obtain higher-order generalized-$\alpha$ schemes. In general, for $k\geq 0$, we seek for $(k+2)$-th order generalized-$\alpha$ method in the form

[TABLE]

where the superscripts $(k)$ represent the $k$-th order derivative in time. The initial conditions are given by $U_{0}^{(j)}=(-\lambda)^{j}U_{0},j=0,1,\cdots,k+1.$ For $k=0,1$, this reduces to the second- and third-order generalized-$\alpha$ methods.

2.2 Third-order accuracy in time

To construct a third-order accurate scheme in the form of (2.3), we need to find the conditions on the parameters. We have the following result.

Theorem 1**.**

Assume that the solution is sufficiently smooth with respect to time. The scheme in (2.3) is third-order accurate in time given

[TABLE]

Proof.

Plugging the last two equations in (2.3) into the first equation, we obtain

[TABLE]

which can be rewritten as a matrix system

[TABLE]

Thus, the amplification matrix becomes

[TABLE]

For the amplification matrix with arbitrary entries, a symbolic calculation verifies that

[TABLE]

where $G_{0}=1$, $G_{1}$ is the trace of $G$ (first invariant), $G_{2}$ is the sum of principal minors of $G$ (second invariant), and $G_{3}$ is the determinant of $G$ (third invariant). We apply the Taylor expansions to obtain

[TABLE]

Substituting (2.11) into (2.10), we obtain

[TABLE]

Assuming sufficient regularity of the solution in time, we apply first and second derivatives to the discrete scheme and obtain

[TABLE]

Substituting these equations into (2.12), we obtain

[TABLE]

The second term in the bracket is of order $\mathcal{O}(\tau^{4})$. Thus, this can be viewed as a part of the local truncation error and we move it to the right hand side. Assuming the same representation on $A_{n}$ and $A_{n+1}$ in (2.3), then $\gamma_{2}=\gamma_{1}$. Thus, we obtain third-order scheme (4th-order local truncation error) when

[TABLE]

which has the solution

[TABLE]

A scheme which has a fourth-order local truncation error leads to a third-order accurate scheme in time. ∎

Remark 1**.**

Similarly, one can show that the scheme in (2.3) has third-order accuracy in time given that

[TABLE]

which is the solution of

[TABLE]

in (2.14). The second equation in either (2.15) or (2.18) is required for obtaining third-order accuracy in time. The second equation in either (2.15) or (2.18) is an auxiliary equation which contributes to a higher order (higher than or equal to $\mathcal{O}(\tau^{4})$) error term in (2.14). A different equation, which is in terms both $\gamma_{1}$ and $\gamma_{2}$, leads to a different third-order scheme. Each scheme has a different unconditionally stable region, which we will discuss in detail in Section 3.

2.3 Higher-order accuracy in time

We follow the derivations in the proof of Theorem 1 to seek higher-order schemes in the form of (2.5). To seek $p$-th order ($p\geq 4$) scheme, we substitute the last two equations in (2.5) into the first equation and obtain a system written in a matrix form

[TABLE]

where

[TABLE]

Similarly, with a slight abuse of notation, the amplification matrix is defined as

[TABLE]

For the amplification matrix with arbitrary entries, one has

[TABLE]

where $G_{j}$ is the sum of the principal minors of order $j$ and this reduces to (2.10) for the third-order ($p=3$) scheme. Applying Taylor expansions of $U_{n+1-j}$ around $U_{n}$ and following the procedure for the third-order scheme, we obtain the following condition for $p$-th order scheme.

[TABLE]

where $C(p)$ is a function of $p$ and some values are given in Table 1. The general pattern as well as the number of possible solutions for each order are open problems and subject to future work.

Remark 2**.**

As for the third-order case, different auxiliary conditions lead to different stability regions. We apply $\gamma_{j}=\gamma_{1},j=2,\cdots,p-1$ so that the last column of the amplification matrix when $\tau\to\infty$ is zero except the last entry, which is an eigenvalue bounded by 1.

3 Stability analysis and eigenvalue control

Section 2 studies the high-order accurate implicit schemes. However, these schemes have conditional and unconditional stability regions. In general, for implicit schemes, unconditionally stable schemes are preferred to the conditionally stable schemes. Thus, we focus on finding the unconditionally stable regions for high-order schemes.

Theorem 2**.**

The third-order scheme (2.3) with $\gamma_{j},j=1,2$ defined in Theorem 1 is unconditionally stable for

[TABLE]

Proof.

To show unconditionally stability, we show the equivalent condition that all the eigenvalues of the amplification matrix $G$ defined in (2.9) are bounded by 1 (strictly less than 1 for repeated real roots) for arbitrary $\lambda$ and $\tau$. The multiplication of $\lambda$ and $\tau$ appears as a single factor in $G$. In general, $\lambda$ is a complex number with positive real part. We denote $T=\lambda\tau$ and let $Re(T)$ vary in $\mathbb{R}^{+}$.

Following closely the analysis on the second-order generalized-$\alpha$ method in (4, 7), firstly, let $Re(T)\to 0,$ the amplification matrix reduces to

[TABLE]

which has eigenvalues

[TABLE]

For a complex eigenvalue, we bound its modulus by 1. Using symbolic calculation, we obtain

[TABLE]

Now, let $Re(T)\to\infty,$ with $\gamma_{1}=\gamma_{2}$, the amplification matrix reduces to

[TABLE]

The eigenvalues are

[TABLE]

Similarly, if $\alpha_{f}\geq 9/16$, then we obtain real eigenvalues, which leads to

[TABLE]

or

[TABLE]

If $\alpha_{f}<9/16$, we obtain complex eigenvalues, which leads to

[TABLE]

Thus, for $Re(T)\to\infty,$ we obtain

[TABLE]

Combining both (3.4) and (3.10) and taking their intersection (as $Re(T)\to\infty$, the stability region is reduced) give the desired results. For finite value of $Re(T)$, we verify symbolically that all the eigenvalues are bounded by 1 in the region defined by (3.1). ∎

The unconditionally stability region is different when we apply $\gamma_{j},j=1,2$ defined in Remark 1 to scheme (2.3). The symbolic analysis is more involved. We show numerically the unconditional stability regions in Figure 1. Figure 1 shows that the third-order scheme (2.3) with (2.17) has a larger stability region.

To control the high-frequency numerical dissipation, following closely the analysis on the second-order generalized-$\alpha$ method in (4, 7), we set that all the eigenvalues at high-frequency, that is when $Re(T)\to\infty$, to be a user-controlled parameter $\rho_{\infty}$ as

[TABLE]

which has a solution

[TABLE]

for $0\leq\rho_{\infty}\leq 1$. One controls the eigenvalues of the amplification matrix by setting $\rho_{\infty}$ and the high-frequency damping is happening when setting $\rho_{\infty}$ close to zero. Other solutions to (3.11) are possible and we refer to A for details.

4 Concluding remarks

The generalized-$\alpha$ method unifies the description of several families of second-order time integrators with the attractive feature of controlling the high-frequency damping. We extend the method to higher-order schemes while maintaining all the attractive features. In particular, at each time step, the second-order generalized-$\alpha$ method solves implicitly one matrix system and then updates the other variables explicitly. This feature is also maintained for the higher-order schemes. Additionally, the third-order method is still a single step method, allowing it to be easily introduced in a time adaptive loop. The generalization of the stability analysis to higher-order schemes as well as the generalization of the schemes for hyperbolic equations will be the subject of future work.

Acknowledgement

This publication was made possible in part by the CSIRO Professorial Chair in Computational Geoscience at Curtin University and the Deep Earth Imaging Enterprise Future Science Platforms of the Commonwealth Scientific Industrial Research Organization, CSIRO, of Australia. Additional support was provided by the European Union’s Horizon 2020 Research and Innovation Program of the Marie Skłodowska-Curie grant agreement No. 777778, the Mega-grant of the Russian Federation Government (N 14.Y26.31.0013), the Institute for Geoscience Research (TIGeR), and the Curtin Institute for Computation. The J. Tinsley Oden Faculty Fellowship Research Program at the Institute for Computational Engineering and Sciences (ICES) of the University of Texas at Austin has partially supported the visits of VMC to ICES. The authors also would like to acknowledge the contribution of an Australian Government Research Training Program Scholarship in supporting this research.