MGPS: Midpoint-Series Group Preserving Scheme for Discretizing Nonlinear Dynamics

Abstract

In this article, we propose a new computational method for initial value problems in ordinary differential equations. The algorithm combines the merits of the group preserving scheme (GPS), which has the ability of avoiding possible spurious solutions utilizing the inherent symmetry group, the cone structure of the nonlinear dynamical system, and the classic midpoint rule. The error and stability analysis are included to demonstrate the convergence properties of the presented method. From the numerical experimental results we obtained, the algorithm can be said to be computationally effective and possesses better simulation ability generally. Meanwhile, it works well with the periodic Hamiltonian system.

1. Introduction

Systems of ordinary differential equations (ODEs) are most versatile for characterizing the behavior of dynamical systems. They are ubiquitous in fields such as science and engineering as well as physical processes, biology, etc. ODEs, partial differential equations (PDEs), as well as other types of equations can then be used to model systems in continuous time (e.g., integro-differential or delay equations). In discrete time, they can be depicted by state machines, Petri nets, finite automata, and so on. As an example, consider the most important hybrid system. Figure 1 illustrates a simple hybrid system: the self-regulating thermostat, which is well characterized by a hybrid automaton. There are two states, $on$ and $off$, denoted by $q_0$ and $q_1$, respectively. The temperature is symbolized by the real-valued variable x. The two ODEs are used to portray the evolution of the continuous operation in the state $q_0$ cooling and $q_1$ state warming, separately. Initially, the system is in the $q_0$ state, i.e., $off$, and the temperature is $x=20$, which drops depending on the flow condition $\dot{x}=-0.1x$. As the temperature falls below 19, the system state switches from $q_0$ to the $on$ state $q_1$. The heater opens and the temperature increases depending on the flow condition $\dot{x}=5-0.1x$. The state switches from $q_1$ to the $off$ $q_0$ if the temperature exceeds 21. Subject to the invariant term $x\ge 18$, the heater is triggered when the temperature falls to 18 degrees at the least.

From this example, one crucial step for various tasks is solving the flow condition that characterizes the continuously dynamic evolution. Formally, ${\varphi }_q:X\times \mathbb{R}\to X$ denotes a flow on $X\subset {\mathbb{R}}^n$ in state q, which is $\dot{x}=f\left(x\right)$ about variable X. $In{v}_q\subseteq X$ is the domain of permitted evolution in q. ${\mathrm{\Delta }}_i\in {\mathbb{R}}^{+}$ is the operation time. ${\gamma }_i:[0,{\mathrm{\Delta }}_i]\to X$ denotes the curve. For all $t\in [0,{\mathrm{\Delta }}_i]$, ${\gamma }_i(t)={\varphi }_{q_i}({\gamma }_i(0),t)$ and ${\gamma }_i(t)\in In{v}_{q_i}$. The subsequent sections cover how our proposed approach can be used to obtain the discretized flow pipes.

Many successes have been obtained in linear systems with various mathematical models such as convex polyhedra [1], griddy polyhedra [2], zonetope [3], ellipsoidal [4], and support function [5] used for approximate simulations expressing the reachable set of hybrid systems. Many traditional tools have emerged for verifying linear hybrid systems, such as HYTECH [6], D/DT [7], CheckMate [8], PHAVER [9], SpaceEx [10], etc.

Nevertheless, for a large number of nonlinear continuous systems, no analytical solution exists and, even if it did, it would be prohibitively expensive to work out and understand its properties. When dealing with continuous and hybrid systems observed over real numbers, the term “exact” is insufficient because there may be numerical errors in observation, noise, and other nonidealities [9,11] over-approximating complex dynamics with simple rectangular flow constraints on the dotted variables, which is operative under some relatively simple nonlinear dynamics. However, the mass of the approximation may be too rough to meet the needs of security verification with comparatively high accuracy requirements.

Applying numerical analysis methods to discrete nonlinear flow constraints is also a proven method. The numerical approximation ”captures“ the behavior of the analytical solution, which is one of the most desirable properties of numerical schemes. Many classical numerical methods have been proposed for this task in the past, such as Euler and Heun, Runge–Kutta, and some others. Normally, the first-order numerical method is easy to implement and low in complexity, but the global error is large and unsuitable for systems with high accuracy requirements. Conversely, higher-order methods such as the Runge–Kutta method can produce numerical solutions with higher accuracy but at the expense of complexity. Furthermore, the spurious solution problem persists in these traditional numerical methods.

As a valuable class of numerical analysis schemes, geometric numerical integration methods serve a vital role since they respect the structure of the problem, which is of a geometric nature. In connection with the geometric characteristics of specific systems, mathematicians and physicists have formulated various numerical algorithms for structure preservation. In [12], a nonstandard finite-difference approach was developed for the Lotka–Volterra system, an essential predator–prey interaction mathematical model. Hamiltonian systems are dominated by the Hamiltonian equations and are prevalent in modern physics in an extensive variety of problems. To preserve the contact transformation property of the Hamilton equations, Ref. [13] devised methods of integration of Hamilton equations and presented many important findings on numerical symplectic integration for the first time, including the symplectic Euler method. Ref. [14] later concludes that the leapfrog method is symplectic. Refs. [15,16] systematically proposed the Hamiltonian equations and the Hamiltonian algorithm (i.e., the symplectic geometric algorithm), and opened up a new field of Hamiltonian algorithms. Ref. [17] discussed the discrete analogue of the gradient of a function and showed how discrete gradients can be used in the numerical integration of ODEs. Ref. [18] implemented a modified discrete gradient method for the purpose of preserving (almost completely) small oscillation periods of any time step. Ref. [19] proposed locally exact discrete gradient numerical integrators, which substitute a function for its step size and exactly preserve the energy integral. Ref. [20] goes further to develop the results of the above literature to the Nth order. Refs. [21,22] presented a new geometric integrator, which combines the implicit midpoint rule (IMR, i.e., second-order implicit Runge–Kutta method) [23] with an appropriate spatial discretization in the PDE cases. It can retain the energy properties of the solution [24], and the function $\mathbf{f}$ needs to be evaluated only once at each step (single-step method, unlike many other Runge–Kutta methods).

Unfortunately, the foregoing numerical approaches either suffer from so-called spurious solutions and ghost fixed points, or are feasible solely for certain systems. In order to avoid the flaws, ref. [25] firstly devised the group preserving scheme (GPS) by utilizing inherent symmetry group and the (null) cone structure of the nonlinear dynamical system. For the aim of improving the accuracy of GPS, an enhanced GPS is proposed in [26]. In this paper, combining the merits of the implicit midpoint rule and GPS, we construct a novelty midpoint group preserving scheme (MGPS). Compared with enhanced GPS, which is also a predictor–corrector method, MGPS overcomes some weaknesses of enhanced GPS.

• The predictor. In enhanced GPS, the Euler method with the same step size is used to calculate $x_{n+1}$, while MGPS uses GPS to calculate $x_{n+1}$, which can obtain the same benefits as the advantages of GPS over the Euler method.
• The corrector. When computing ${\eta }_n^e$, the enhanced GPS replaces $x_n$ with ${\overline{x}}_n$; the MGPS remains unchanged and is still used to calculate ${\eta }_n^e$. This improves the accuracy while ensuring that the cone structure is held.

The remaining parts of this paper are structured below. Section 2 firstly introduces GPS and analyzes its stability property. Then, we devise the midpoint group preserving scheme (MGPS) and details of error and stability properties. Some numerical tests are conducted to assess the effectiveness of the proposed solver in Section 3. The numerical results are very encouraging. This demonstrates improved performance and higher accuracy compared with other well-known second-order methods present in the literature. Finally, Section 4 sums up this work and addresses the path for the future.

2. Discretization

Discretization of the nonlinear continuous dynamic evolution is a vital measure to understand the properties of the system while differential equations have no explicit solutions. Consider the following ordinary differential equation with the initial value

$$\dot{\mathbf{x}}(t)=\mathbf{f}(t,\mathbf{x}(t)),\mathbf{x}(t_0)={\mathbf{x}}_0,t\in \mathbb{R},\mathbf{x}\in {\mathbb{R}}^k$$

(1)

It is commonly utilized to demonstrate the attributes and behaviors of dynamic systems. $\mathbf{x}$ represents k-dimensional state vector, and $\mathbf{f}$ denotes a vector function about $\mathbf{x}$ and t. The Lipschitz condition is assumed to ensure the existence and uniqueness of the solution:

$$\begin{array}{c}\hfill |\mathbf{f}(t,{\mathbf{x}}_1)-\mathbf{f}(t,{\mathbf{x}}_2)|\le L(|{\mathbf{x}}_1-{\mathbf{x}}_2|)\forall (t,{\mathbf{x}}_1),(t,{\mathbf{x}}_2)\in \mathbb{R}\times {\mathbb{R}}^k\end{array}$$

(2)

Suppose that the running time delay ${\Delta }_i$ is uniformly partitioned into n steps. A sequence of time-value pairs $\{(t_j,u_j),j\in \{0,\dots ,n-1\}\}$ are generated along a solution curve. Take the Euler method as an example, we can obtain curve ${\gamma }_i(\frac{j}{n}·{\mathrm{\Delta }}_i)\in 〚In{v}_{q_i}〛$ at each time point $t_j=\frac{j}{n}·{\mathrm{\Delta }}_i$, where $〚In{v}_{q_i}〛$ denotes the collection that meets $In{v}_{q_i}$. The continuous trajectory is thus converted into a finite discrete-time sequence in this manner: $\eta =\langle q_1,{\mathrm{\Delta }}_1·0,{\gamma }_1(0)\rangle$$\stackrel{t}{\to }\langle q_1,{\mathrm{\Delta }}_1·\frac{1}{n},{\gamma }_1({\mathrm{\Delta }}_1·\frac{1}{n})\rangle \stackrel{t}{\to }\cdots \stackrel{t}{\to }\langle q_1,{\mathrm{\Delta }}_1·\frac{n-1}{n},{\gamma }_1({\mathrm{\Delta }}_1·\frac{n-1}{n})\rangle \stackrel{r_{q_1,q_2}}{⟶}\langle q_2,{\mathrm{\Delta }}_2·0,{\gamma }_2(0)\rangle \stackrel{t}{\to }\cdots$.

Unfortunately, spurious solutions and ghost fixed points are broadly present in classical numerical methods and may yield unpredictable outcomes. The problem of such spurious solutions schemes is studied in [27] when first-order differential equations are computed employing the Runge–Kutta. Ref. [25] devised a group preserving scheme in the Minkowski space utilizing homogeneous coordinates to avoid the shortcomings, which serve as the foundation of our approach.

2.1. Group Preserving Scheme

The idea of GPS is that embedding Equation (1) into $k+1$ dimensional augmented dynamical system yields the equation below.

$$\frac{d}{dt}\left[\begin{array}{c}\mathbf{x}\\ \parallel \mathbf{x}\parallel \end{array}\right]=\left[\begin{array}{cc}{\mathbf{0}}_{k\times k}& \frac{\mathbf{f}(t,\mathbf{x})}{\parallel \mathbf{x}\parallel }\\ \frac{{\mathbf{f}}^{\mathbf{T}}(t,\mathbf{x})}{\parallel \mathbf{x}\parallel }& 0\end{array}\right]\left[\begin{array}{c}\mathbf{x}\\ \parallel \mathbf{x}\parallel \end{array}\right]$$

(3)

The first equation in (3) is obviously identical to the original Equation (1), and with the addition of the second equation, we obtain an enhancement of $\mathbf{X}:={({\mathbf{x}}^{\mathbf{T}},\parallel \mathbf{x}\parallel )}^{\mathbf{T}}$ Minkowskian structure of state variables meeting the cone condition

$${\mathbf{X}}^{\mathbf{T}}\mathbf{gX}=0$$

(4)

Herein, $\mathbf{g}$ denotes Minkowski metrics and is expressed by

$$\mathbf{g}=\left[\begin{array}{cc}{\mathbf{I}}_{\mathbf{k}}& {\mathbf{0}}_{k\times 1}\\ {\mathbf{0}}_{1\times k}& -1\end{array}\right]$$

(5)

where ${\mathbf{I}}_{\mathbf{k}}$ signifies the k-order identity matrix. The inference below leads to this conclusion:

$${\mathbf{X}}^{\mathbf{T}}\mathbf{gX}=\mathbf{x}·\mathbf{x}-{\parallel \mathbf{x}\parallel }^2={\parallel \mathbf{x}\parallel }^2-{\parallel \mathbf{x}\parallel }^2=0$$

(6)

Then, as a result of the previously stated cone condition, Equation (3) has now become

$$\dot{\mathbf{X}}=\mathbf{A}(t,\mathbf{X})\mathbf{X}$$

(7)

where

$$\mathbf{A}(t,\mathbf{X}):=\left[\begin{array}{cc}{\mathbf{0}}_{k\times k}& \frac{\mathbf{f}(t,\mathbf{x})}{\parallel \mathbf{x}\parallel }\\ \frac{{\mathbf{f}}^{\mathbf{T}}(t,\mathbf{x})}{\parallel \mathbf{x}\parallel }& 0\end{array}\right]$$

(8)

It is a (local) Lie algebra of the proper orthochronous Lorentz group ${\mathbf{SO}}_o(k,1)$ that fulfills

$${\mathbf{A}}^{\mathbf{T}}\mathbf{g}+\mathbf{gA}=0$$

(9)

This motivates the development of the so-called group preserving scheme with discretized mapping $\mathbf{G}$ that preserves exactly the following attributes [25]:

$$\left\{\begin{array}{cc}\hfill & {\mathbf{G}}^{\mathbf{T}}\mathbf{gG}=\mathbf{g}\hfill \\ \hfill & \det \mathbf{G}=1\hfill \\ \hfill & G_0^0>0\hfill \end{array}\right.$$

(10)

$G_0^0$ denotes the 00-th element of $\mathbf{G}$, whereas det signifies the abbreviation for determinant. Provided that ${\mathbf{X}}_n$ is the value of $\mathbf{X}$ at $t=t_n$, it is required to work out ${\mathbf{X}}_{n+1}$ at $t=t_{n+1}$. It can be known from Equation (8) that $\mathbf{A}$ is not a constant matrix. Let $h=t_{n+1}-t_n$. Then, the GPS is formulated by

$${\mathbf{X}}_{n+1}=\mathbf{G}(n){\mathbf{X}}_n$$

(11)

$\mathbf{G}(n)\in {\mathbf{SO}}_o(k,1)$ denotes the group value at $t_n$. The rest of the issue is determining how to obtain the expression of $\mathbf{G}(n)$. The exponential mapping is adopted to calculate $\mathbf{G}(n)$ with the next expression:

$$\mathbf{G}(n)=exp\left[h\mathbf{A}(n)\right]=\left[\begin{array}{cc}{\mathbf{I}}_k+\frac{(a_n-1)}{{\parallel {\mathbf{f}}_n\parallel }^2}{\mathbf{f}}_n{\mathbf{f}}_n^{\mathbf{T}}& \frac{b_n{\mathbf{f}}_n}{\parallel {\mathbf{f}}_n\parallel }\\ \frac{b_n{\mathbf{f}}_n^{\mathbf{T}}}{\parallel {\mathbf{f}}_n\parallel }& a_n\end{array}\right]$$

(12)

where

$$a_n:=cosh(\frac{h\parallel {\mathbf{f}}_n\parallel }{\parallel {\mathbf{x}}_n\parallel }),b_n:=sinh(\frac{h\parallel {\mathbf{f}}_n\parallel }{\parallel {\mathbf{x}}_n\parallel })$$

(13)

The first line of the formula (11) is picked up and one obtains

$${\mathbf{x}}_{n+1}={\mathbf{x}}_n+{\eta }^e{\mathbf{f}}_n$$

(14)

where

$${\eta }^e:=\frac{(a_n-1){\mathbf{f}}_n·{\mathbf{x}}_n+b_n\parallel {\mathbf{x}}_n\parallel \parallel {\mathbf{f}}_n\parallel }{{\parallel {\mathbf{f}}_n\parallel }^2}$$

(15)

From $a_n>1,\forall \mathbf{h}>0$ and $\parallel {\mathbf{f}}_n\parallel \parallel {\mathbf{x}}_n\parallel \ge {\mathbf{f}}_n·{\mathbf{x}}_n\ge -\parallel {\mathbf{f}}_n\parallel \parallel {\mathbf{x}}_n\parallel$, the inequalities can be concluded:

$$\left\{\begin{array}{c}{\eta }^e\le \left[exp(\frac{h\parallel {\mathbf{f}}_n\parallel }{\parallel {\mathbf{x}}_n\parallel })-1\right]\frac{\parallel {\mathbf{x}}_n\parallel }{\parallel {\mathbf{f}}_n\parallel }\hfill \\ {\eta }^e\ge \left[1-exp(-\frac{h\parallel {\mathbf{f}}_n\parallel }{\parallel {\mathbf{x}}_n\parallel })\right]\frac{\parallel {\mathbf{x}}_n\parallel }{\parallel {\mathbf{f}}_n\parallel }\hfill \end{array}\right.$$

(16)

Stable Analysis

This section will analyze the stability of the GPS (14). Regularly, the following test equation is used for analyzing numerical schemes’s stability.

$$\begin{array}{c}\hfill {\mathbf{x}}^{\prime }=\lambda \mathbf{x},\mathbf{x}(0)={\mathbf{x}}_0\end{array}$$

(17)

whereby $\lambda$ is a complex number. It is also assumed that ${\mathbf{x}}_0\ne 0$, which would otherwise infer trivial zero solutions. Putting (17) into (14), one obtains

$$\begin{array}{cc}\hfill & {\mathbf{x}}_1={\mathbf{x}}_0+{\eta }_0^e{\mathbf{f}}_0={\mathbf{x}}_0+{\eta }_0^e(\lambda {\mathbf{x}}_0)=(1+\lambda {\eta }_0^e){\mathbf{x}}_0\hfill \\ \hfill & {\mathbf{x}}_2={\mathbf{x}}_1+{\eta }_1^e{\mathbf{f}}_1={\mathbf{x}}_1+{\eta }_1^e(\lambda {\mathbf{x}}_1)=(1+\lambda {\eta }_0^e)(1+\lambda {\eta }_1^e){\mathbf{x}}_0\hfill \end{array}$$

Inductively,

$$\begin{array}{cc}\hfill {\mathbf{x}}_{n+1}& ={\mathbf{x}}_0\prod _{i=1}^n(1+\lambda {\eta }_i^e)\hfill \end{array}$$

(18)

is derived. The analytical solution $\mathbf{x}={\mathbf{x}}_0{e}^{\lambda t}$ ultimately converges to 0 as t goes to infinity under the restriction $Re(\lambda )<0$. In numerical terms, it is corresponding to the situation where ${\mathbf{x}}_n$ becomes 0 when n becomes infinite in Equation (18). In order to meet this condition, it needs

$$\begin{array}{c}\hfill \prod _{i=1}^n|(1+\lambda {\eta }_i^e)|<1\end{array}$$

(19)

Besides, the following Equations (20)–(25) could be clearly deduced.

$$\begin{array}{cc}\hfill \mathbf{x}·\mathbf{x}& ={\parallel \mathbf{x}\parallel }^2\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill \parallel \mathbf{f}\parallel & =\sqrt{(\lambda \mathbf{x})(\lambda \mathbf{x})}=|\lambda |\parallel \mathbf{x}\parallel \hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill \mathbf{f}·\mathbf{f}& ={\parallel \mathbf{f}\parallel }^2={|\lambda |}^2{\parallel \mathbf{x}\parallel }^2\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill \mathbf{f}·\mathbf{x}& =\lambda \mathbf{x}·\mathbf{x}=\lambda {\parallel \mathbf{x}\parallel }^2\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill \parallel \mathbf{x}\parallel ·\parallel \mathbf{f}\parallel & =\parallel \mathbf{x}\parallel ·\parallel \lambda \mathbf{x}\parallel ={|\lambda |\parallel \mathbf{x}\parallel }^2\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill b_n-a_n& =-e^{\frac{h\parallel {\mathbf{f}}_n\parallel }{\parallel {\mathbf{x}}_n\parallel }}=-e^{-h|\lambda |}\hfill \end{array}$$

(25)

Hence, Equation (15) can be transformed into

$$\begin{array}{cc}\hfill {\eta }_n^e& =\frac{(a_n-1)\lambda {\parallel {\mathbf{x}}_n\parallel }^2+b_n|\lambda |{\parallel {\mathbf{x}}_n\parallel }^2}{{|\lambda |}^2{\parallel {\mathbf{x}}_n\parallel }^2}\hfill \\ \hfill & =\frac{(a_n-1)\lambda +b_n|\lambda |}{{|\lambda |}^2}\hfill \end{array}$$

(26)

In practice, $\lambda \in \mathbb{R}$ is often restricted; then, when $\lambda <0$, Equation (26) could be transformed into

$$\begin{array}{cc}\hfill {\eta }_n^e& =\frac{(a_n-1)\lambda -b_n\lambda }{(-\lambda )(-\lambda )}\hfill \\ \hfill & =\frac{b_n-(a_n-1)}{-\lambda }\hfill \\ \hfill & =\frac{1-e^{-h|\lambda |}}{-\lambda }\hfill \\ \hfill & =\frac{1-e^{h\lambda }}{-\lambda }\hfill \end{array}$$

(27)

Consequently, Equation (18) is deduced

$${\mathbf{x}}_{n+1}={\mathbf{x}}_0\prod _{i=1}^n(1+e^{h\lambda }-1)={\mathbf{x}}_0{(e^{h\lambda })}^n$$

(28)

Thus, $e^{h\lambda }<1$ for $\forall h>0$, which implies that the interval of absolute stability contains $(-\infty ,0)$.

2.2. Midpoint Group Preserving Scheme

As the well-known second-order Runge–Kutta method, the implicit midpoint rule is a predictor–corrector numerical algorithm. The predicted value at $t_{n+\frac{1}{2}}$ is approximated by Euler’s method with a half step [23]

$$x_{n+1}=x_n+h\mathbf{f}(t_{n+\frac{1}{2}},\frac{{\mathbf{x}}_n+{\mathbf{x}}_{n+1}}{2}),t_{n+\frac{1}{2}}=t_n+\frac{h}{2}$$

(29)

The improvement of the midpoint method is to as accurately as possible obtain some of the curvature that will occur in the solution before calculating ${\mathbf{x}}_{n+1}$. Integrating the notion of GPS (14), step size h is replaced with ${\eta }_n^e$. Thus, the midpoint group preserving scheme can be formulated naturally

$$\left\{\begin{array}{cc}\hfill {\mathbf{x}}_{n+1}& ={\mathbf{x}}_n+{\eta }_n^e{\mathbf{f}}_n\hfill \\ \hfill {\mathbf{x}}_{n+1}& ={\mathbf{x}}_n+{\eta }_n^e\mathbf{f}(t_{n+\frac{1}{2}},\frac{{\mathbf{x}}_n+{\mathbf{x}}_{n+1}}{2})\hfill \end{array}\right.$$

(30)

where ${\eta }_n^e$ is the same as the definition of (15). In the following subsections, we will discuss its error estimator and stability analysis.

2.2.1. Error Analysis

Firstly, we consider the error associated with the proposed algorithm (30). Let the local truncation error be $e_n$; then,

$$\mathbf{x}(t_{n+1})=\mathbf{x}(t_n)+{\eta }_n^e\mathbf{f}(t_{n+\frac{1}{2}},\frac{\mathbf{x}(t_n)+\mathbf{x}(t_{n+1})}{2})+e_n$$

(31)

We slightly modify the right-hand side and write

$$\mathbf{x}(t_{n+1})=\mathbf{x}(t_n)+{\eta }_n^e\mathbf{f}(t_{n+\frac{1}{2}},\mathbf{x}(t_{n+\frac{1}{2}}))+e_{n_1}+e_n$$

(32)

Then, by Lipschitz condition (2),

$$\begin{array}{cc}\hfill |e_{n_1}|& ={\eta }_n^e|\mathbf{f}(t_{n+\frac{1}{2}},\frac{\mathbf{x}(t_n)+\mathbf{x}(t_{n+1})}{2})-\mathbf{f}(t_{n+\frac{1}{2}},\mathbf{x}(t_{n+\frac{1}{2}}))|\hfill \\ \hfill & \le \frac{1}{2}{\eta }_n^{e}L|\mathbf{x}(t_n)+\mathbf{x}(t_{n+1})-2\mathbf{x}(t_{n+\frac{1}{2}})|\hfill \end{array}$$

(33)

and expanding $\mathbf{x}(t_{n+1})$ and $\mathbf{x}(t_{n+\frac{1}{2}})$ in Taylor series about point $t_n$ with $\mathbf{x}=\mathbf{x}(t_n)$, we have

$$|e_{n_1}|\le \frac{1}{2}{\eta }_n^{e}L|\mathbf{x}+(\mathbf{x}+h{\mathbf{x}}^{\prime })-2(\mathbf{x}+\frac{1}{2}h{\mathbf{x}}^{\prime })+\mathcal{O}(h^2)|$$

(34)

In practice, since $|{\eta }_n^e-h|=\mathcal{O}(h)$, the right side of (34) is $\mathcal{O}(h^3)$. Now, from (32), since $\mathbf{f}(t,x(t))=x^{\prime }(t)$, we obtain

$$\begin{array}{cc}\hfill e_{n_1}+e_n& =\mathbf{x}(t_{n+1})-\mathbf{x}(t_n)-{\eta }_n^e{\mathbf{x}}^{\prime }(t_{n+\frac{1}{2}})\hfill \\ \hfill & =h{\mathbf{x}}^{\prime }+\frac{1}{2}{h}^2{\mathbf{x}}^{″}-{\eta }_n^e({\mathbf{x}}^{\prime }+\frac{1}{2}h{\mathbf{x}}^{″})+\mathcal{O}(h^3)\hfill \\ \hfill & =\mathcal{O}(h^3)\hfill \end{array}$$

(35)

Hence,

$$|e_n|\le |e_{n_1}|+\mathcal{O}(h^3)\le c{h}^3$$

(36)

Thus, local truncation error $e_n$ is bounded.

Next, let the global truncation error $E_n=\mathbf{x}(t_n)-{\mathbf{x}}_n$. The recurrence relation for $En$ is obtained by subtracting (30) from (31)

$$E_{n+1}=E_n+{\eta }_n^e(\mathbf{f}(t_{n+\frac{1}{2}},\frac{\mathbf{x}(t_n)+\mathbf{x}(t_{n+1})}{2})-\mathbf{f}(t_{n+\frac{1}{2}},\frac{{\mathbf{x}}_n+{\mathbf{x}}_{n+1}}{2}))+e_n$$

(37)

The Lipschitz condition implies

$$\begin{array}{cc}\hfill |E_{n+1}|& \le |E_n|+\frac{1}{2}{\eta }_n^{e}L|\frac{\mathbf{x}(t_n)+\mathbf{x}(t_{n+1})}{2}-\frac{{\mathbf{x}}_n+{\mathbf{x}}_{n+1}}{2}|+|e_n|\hfill \\ \hfill & \le |E_n|+\frac{1}{2}{\eta }_n^{e}L(|E_n|+|E_{n+1}|)+c{h}^3\hfill \end{array}$$

(38)

Thereby,

$$|E_{n+1}|\le \frac{1+\frac{1}{2}{\eta }_n^{e}L}{1-\frac{1}{2}{\eta }_n^{e}L}|E_n|)+c{h}^3\le \frac{1+\frac{1}{2}(h+c_n\mathcal{O}(h^2))L}{1-\frac{1}{2}(h-c_n\mathcal{O}(h^2))L}+c{h}^3{c}_n>0$$

(39)

Let $R=max\{c_0\mathcal{O}(h^2),c_1\mathcal{O}(h^2),\cdots ,c_n\mathcal{O}(h^2)\}$ and $r=min\{c_0\mathcal{O}(h^2),c_1\mathcal{O}(h^2),\cdots ,$$c_n\mathcal{O}(h^2)\}$; then,

$$|E_{n+1}|\le \frac{1+\frac{1}{2}(h+R)L}{1-\frac{1}{2}(h-r)L}+c{h}^3$$

(40)

Now, the relations of the form

$$|E_{n+1}|\le A|E_n|+B,n\ge 0,E_0=0$$

(41)

provide the estimate

$$|E_{n+1}|\le B(1+A+A^2+\cdots +A^n)=B\frac{A^{n+1}-1}{A-1}$$

(42)

Since $\frac{1}{2}\le 1-\frac{1}{2}(h-r)L\le 1$ for small h, we have

$$A-1=\frac{(h+\frac{1}{2}(R+r))L}{1-\frac{1}{2}(h-r)L}\ge (h+\frac{1}{2}(R+r))L$$

(43)

$$A^{n+1}={(1+\frac{(h+\frac{1}{2}(R+r))L}{1-\frac{1}{2}(h-r)L})}^{n+1}\le {(1+2(h+\frac{1}{2}(R+r))L)}^{n+1}\le e^{2(h+\frac{1}{2}(R+r))L(n+1)}$$

(44)

whence

$$|E_{n+1}|\le \frac{1}{(h+\frac{1}{2}(R+r))L}{e}^{2(h+\frac{1}{2}(R+r))L(n+1)}c{h}^3\le c_1{h}^3$$

(45)

with a constant $c_1$ that only depends on L and t, i.e., we have uniform convergence in any finite interval [0, t].

2.2.2. Stability Analysis

Just as previously discussed in Section 2, according to Equation (17), we obtain

$$\begin{array}{cc}\hfill {\mathbf{x}}_{n+1}& ={\mathbf{x}}_n+{\eta }_n^e\lambda ({\mathbf{x}}_n+\frac{{\eta }_n^e}{2}\lambda {\mathbf{x}}_n)\hfill \end{array}$$

(46)

$$\begin{array}{cc}\hfill & ={\mathbf{x}}_n(1+{\eta }_n^e\lambda +\frac{1}{2}{({\eta }_n^e\lambda )}^2)\hfill \end{array}$$

(47)

Similarly, if $\lambda \in \mathbb{R}$ is adopted and $\lambda <0$; then, substituting Equation (27) into Equation (30)

$${\mathbf{x}}_{n+1}={\mathbf{x}}_n(\frac{e^{h\lambda }+1}{2})$$

(48)

For any $h>0$, $0<\frac{e^{h\lambda }+1}{2}<1$; thus, it also has intervals of absolute stability containing $(-\infty ,0)$.

2.3. Scalability

As a matter of fact, along this line of thought, more classes of midpoint-like group preserving schemes can be extended.

$$\left\{\begin{array}{cc}\hfill {\mathbf{x}}_{n+1}& ={\mathbf{x}}_n+{\eta }_{n,\frac{1}{2}}^e{\mathbf{f}}_n\hfill \\ \hfill {\overline{\mathbf{x}}}_{n+1}& ={\mathbf{x}}_n+{\eta }_{n,\frac{1}{2}}^e{\mathbf{f}}_n({\mathbf{x}}_{n+1})\hfill \\ \hfill {\mathbf{x}}_{n+1}& ={\mathbf{x}}_n+{\eta }_n^e\mathbf{f}(t_{n+\frac{1}{2}},\frac{{\mathbf{x}}_n+{\overline{\mathbf{x}}}_{n+1}}{2})\hfill \end{array}\right.$$

(49)

Comparing Equations (30) and (49), it is intuitively clear that the prediction function has been further optimized, while the correction function remains in the same format. This optimization is worthwhile as evidenced by the following examples 3 and 4. This scheme is shorthanded as OMGPS hereafter.

3. Experiments

In this section, the experiments will be conducted in two types of examples: one for regular continuous dynamical systems (Examples 1 and 2) and the other for two special types of systems, Kolmogorov systems (Example 3) and Hamiltonian systems (Example 4). The proposed method and other second-order methods are adopted and compared. All tests are performed in MATLAB R2020a.

3.1. Example 1

The following nonlinear differential equation is use to describe population growth in a nonlinear hybrid system [11]

$$\left\{\begin{array}{c}\dot{x_1}=(2000-x_2-5x_1)x_1\hfill \\ \dot{x_2}=(4x_1-2600-4x_2)x_2\hfill \end{array}\right.$$

(50)

where initial value ${\mathbf{x}}_0=\{900,150\}$ and time $t=0.008$. Its true solutions are simulated using explicit Euler scheme with a step size of $h=0.00001$, and other numerical methods have the same step size $h=0.001$. Figure 2 shows their computation results. It is evident that GPS deviates from the real curve by a large margin. The simulation outcomes of midpoint, Heun, enhanced-GPS, and MGPS are quite close to each other. By comparison clearly, OMGPS is on top of them. This is also in accordance with the optimization point of Section 2.3.

3.2. Example 2

The second example is taken from the literature [26], which is a two-dimensional continuous system. Starting from initial set ${\mathbf{x}}_0=\{1,0\}$, the system evolves from time 0 to 25.

$$\left\{\begin{array}{c}\dot{x_1}=x_2\hfill \\ \dot{x_2}=-x_1+\frac{1}{3}{x}_1^3-x_2\hfill \end{array}\right.$$

(51)

Similarly, applying explicit Euler method with $h=0.0001$ simulates the true curve. $h=0.5$ is adopted for the remaining methods. The simulation results are exhibited in Figure 3. Apparently, The conclusion reached is identical to the first example.

3.3. Example 3

The Lotka–Volterra system is an important mathematical model for depicting predator–prey interactions in the dynamics of biological systems. The Lotka–Volterra equations can be written simply as follows:

$$\left\{\begin{array}{c}\dot{x_1}=a{x}_1-b{x}_1{x}_2\hfill \\ \dot{x_2}=-c{x}_2+d{x}_1{x}_2\hfill \end{array}\right.$$

(52)

It was mentioned earlier that [12] proposes a nonstandard finite-difference scheme (abbreviated here as NSFD) for this system, and to facilitate comparison, the same parameters as in [12] are employed here for the test, namely, $a=b=c=d=1$. As can be seen in Figure 4, the numerical integration results of MGPS form a closed curve, which is consistent with the results of NSFD and is very close to the true value, which is simulated using the Euler method with a step size of 0.0001. The Figure 5 exhibits more clearly the error accuracy of several methods, and the superiority of MGPS is noticeable.

3.4. Example 4

The last example is oriented towards Hamiltonian systems and comes from the literature [19], which demonstrates the GR family of methods($\mathrm{GR},\mathrm{MOD}\_\mathrm{GR},\mathrm{GR}\_\mathrm{LEX},\mathrm{GR}\_\mathrm{SLEX}$). From this, we chose $\mathrm{GR}\_\mathrm{SLEX}$ for comparison. The Hamiltonian system is the simple pendulum $H=\frac{p^2}{2}-cosx$. The corresponding ODE is

$$\left\{\begin{array}{c}\dot{x}=H_p=p\hfill \\ \dot{p}=-H_x=-sinx\hfill \end{array}\right.$$

(53)

Here, two sets of experiments were conducted adopting the same parameters as those listed above in the literature—that is, $h=0.25$ and $h=0.01$. The rest of parameters are $x_0=0,p_0=1.8,t=20$. In Figure 6, GPS performs very poorly and yields erroneous results. In Figure 7, however, the true value is better simulated, but the accuracy in comparison with the other two methods is inferior. This illustrates the fact that GPS is less stable in the Hamilton system. When dealing with this type of problem with GPS, a relatively small h is recommended. In both sets of experiments, MGPS and $\mathrm{GR}\_\mathrm{SLEX}$ are reliably capable of modeling the true values. Figure 8 and Figure 9 express that $\mathrm{GR}\_\mathrm{SLEX}$ has better accuracy. With the step size becoming smaller, the precision gap of MGPS and $\mathrm{GR}\_\mathrm{SLEX}$ is narrowing, and it can be seen from Figure 9 that MGPS nearly has the same level of error accuracy as $\mathrm{GR}\_\mathrm{SLEX}$. Regarding the running efficiency, we executed two sets of trials—one with 20,000 iterations and the other with 200,000 iterations, to minimize the effect of other factors—and averaged the running time after three runs for each group. The outcome of the experimental data is placed in Appendix A and shows that when the number of iterations is not large, the running times are comparable, and vice versa, $\mathrm{GR}\_\mathrm{SLEX}$ is superior. Although $MGPS$ is not computed serially, this strength cannot be fully exploited at low dimensions.

4. Conclusions

Computing flow pipes is a highly significant problem in the verification of hybrid systems. This paper constructs an easy-to-implement and extensible midpoint-series group preserving scheme from a numerical analysis point of view, making use of the classical midpoint method and the symmetry of the group structure. Yet, the verification of hybrid systems faces several more problems than this one. How to effectively detect the termination of continuous evolution is another very important issue. To the best of our knowledge, there is no viable numerical approach. How to better solve this problem and the solution of high-dimensional nonlinear dynamic equations will be the focus of our future work.