Preservation of the invariants of Lotka-Volterra equations by iterated deferred correction methods

TL;DR

This paper demonstrates that applying iterated deferred correction methods to Kahan's discretization of 3D Lotka-Volterra equations preserves invariants and periodicity more accurately, while significantly speeding up computations.

Contribution

It introduces a novel combination of Kahan's discretization with deferred correction methods to enhance invariant preservation and computational efficiency in Lotka-Volterra systems.

Findings

Invariants are preserved with increasing accuracy using deferred correction.

The method achieves substantial speedups over Kahan's original method.

Periodic solutions remain well-preserved over long-term integration.

Abstract

In this paper we apply Kahan's nonstandard discretization to three dimensional Lotka-Volterra equations in bi-Hamiltonian form. The periodicity of the solutions and all polynomial and non-polynomial invariants are well preserved in long-term integration. Applying classical deferred correction method, we show that the invariants are preserved with increasing accuracy as a results of more accurate numerical solutions. Substantial speedups over the Kahan's method are achieved at each run with deferred correction method.

Tables

Table 1. Table 1: Example  5.1 : Convergence rates and speed-ups over the system without CDC

$S$/$n$	5	6	7	8	9	10	11	$\Delta t$	Speed-Up
1	4.48	4.53	4.65	4.47	4.75	4.72	4.83	0.01	3.0
2	4.51	6.99	6.78	7.07	6.82	6.85	6.79	0.05	7.1
3	4.84	6.95	6.76	8.89	8.53	10.73	10.39	0.15	10.7
4	4.84	7.02	6.76	8.89	9.29	11.02	11.03	0.25	11.4

Equations

$${\dot{u}}_i=u_i(r_i+\sum _j{\alpha }_{ij}{u}_j),i=1,\dots ,m,$$

$$\dot{\mathbf{u}}=J_1\nabla H_2=J_2\nabla H_1,$$

$$\dot{\mathbf{u}}=\nabla H_1\times \nabla H_2.$$

$${\dot{u}}_1$$

$${\dot{u}}_2$$

$${\dot{u}}_3$$

$$J_{1}=\left(\begin{array}[]{ccc}0&cu_{1}u_{2}&bcu_{1}u_{3}\\ -cu_{1}u_{2}&0&-u_{2}u_{3}\\ -bcu_{1}u_{3}&u_{2}u_{3}&0\end{array}\right),$$

$$J_{2}=\left(\begin{array}[]{ccc}0&cu_{1}u_{2}(au_{3}+\mu)&cu_{1}u_{3}(u_{2}+\nu)\\ -cu_{1}u_{2}(au_{3}+\mu)&0&u_{1}u_{2}u_{3}\\ -cu_{1}u_{3}(u_{2}+\nu)&-u_{1}u_{2}u_{3}&0\end{array}\right),$$

$$H_1=ab\ln u_1-b\ln u_2+\ln u_3,H_2=ab{u}_1+u_2-a{u}_3+\nu \ln u_2-\mu \ln u_3.$$

$${\dot{u}}_1$$

$${\dot{u}}_2$$

$${\dot{u}}_3$$

$$J_{1}=\left(\begin{array}[]{ccc}0&-1&1\\ 1&0&-1\\ -1&1&0\end{array}\right),\quad J_{2}=\left(\begin{array}[]{ccc}0&u_{1}u_{2}&-u_{1}u_{3}\\ -u_{1}u_{2}&0&u_{2}u_{3}\\ u_{1}u_{3}&-u_{2}u_{3}&0\end{array}\right),$$

$$\dot{\mathbf{u}}=f(\mathbf{u}):=Q(\mathbf{u})+B\mathbf{u},$$

$$\frac{{\mathbf{u}}_{n+1}-{\mathbf{u}}_n}{\Delta t}=Q({\mathbf{u}}_n,{\mathbf{u}}_{n+1})+\frac{1}{2}B({\mathbf{u}}_n+{\mathbf{u}}_{n+1}),$$

$$Q({\mathbf{u}}_n,{\mathbf{u}}_{n+1})=\frac{1}{2}\left(Q({\mathbf{u}}_n+{\mathbf{u}}_{n+1})-Q({\mathbf{u}}_n)-Q({\mathbf{u}}_{n+1})\right).$$

$$\frac{{\mathbf{u}}_{n+1}-{\mathbf{u}}_n}{\Delta t}$$

$$\left(I-\frac{\Delta t}{2}{f}^{\prime }({\mathbf{u}}_n)\right)\tilde{\mathbf{u}}=\Delta tf({\mathbf{u}}_n),{\mathbf{u}}_{n+1}={\mathbf{u}}_n+\tilde{\mathbf{u}}.$$

$$\frac{{\mathbf{u}}_{n+1}-{\mathbf{u}}_n}{\Delta t}=-\frac{1}{2}f({\mathbf{u}}_n)+2f\left(\frac{{\mathbf{u}}_{n+1}+{\mathbf{u}}_n}{2}\right)-\frac{1}{2}f({\mathbf{u}}_{n+1}).$$

$$0=t_0<t_1<\cdots <t_j<\cdots <t_J=T,t_{j+1}-t_j=\Delta t,$$

$$t_j=t_{j,1}<t_{j,2}<\cdots <t_{j,i}<\cdots <t_{j,n}=t_{j+1},$$

$${\mathbf{u}}_{j,i}^{[0]}=\mathbf{u}(t_{j,i})+\mathcal{O}((\Delta t{)}^{p_0}),i=1,\dots ,n.$$

$$\frac{d}{dt}{\mathbf{e}}_j^{[s-1]}(t,{\mathbf{U}}_j^{[s-1]}(t))$$

$${\mathbf{e}}_j^{[s-1]}(t_j,{\mathbf{U}}_j^{[s-1]}(t_j))$$

$${\mathbf{e}}_j^{[s]}(t,{\mathbf{U}}_j^{[s]}(t))=\mathbf{u}(t)-{\mathbf{U}}_j^{[s]}(t),t\in [t_j,t_{j+1}].$$

$${\mathbf{U}}_j^{[s]}(t)=\sum _{k=1}^n{l}_k(t)\cdot {\mathbf{u}}_{j,k}^{[s]},l_k(t)=\prod _{i\ne k}\frac{t-t_{j,i}}{t_{j,k}-t_{j,i}},$$

$${\mathbf{U}}_j^{[s]}={\mathbf{U}}_j^{[s-1]}+{\mathbf{E}}_j^{[s-1]}.$$

$$L^2(H)={\left(\Delta t\sum _{i=1}^J[H(t_i)-H(0){]}^2\right)}^{1/2},L^2(\mathbf{u})={\left(\Delta t\sum _{i=1}^J[{\mathbf{u}}_i-{\mathbf{u}}_{exact}(t_i){]}^2\right)}^{1/2},$$

$$\text{order}=\frac{1}{\log 2}\log \left(\frac{{\text{Err}}_{\Delta t_i}}{{\text{Err}}_{\Delta t_{i+1}}}\right),$$

Full text

1 Introduction

In the last two decades, many structure preserving geometric integrators are developed to preserve symplectic structure, energy and other invariants, phase space volume, reversing symmetries, dissipation approximately or exactly (up to the round-off errors) [1, 2] of dynamical systems. These are symplectic and variational integrators for Hamiltonian systems [1, 3], integral preserving methods [2] and discrete gradient methods [4]. For special classes of ordinary differential equations (ODEs), there exist non-standard discretization methods [5, 6] which preserve the conserved quantities and other features approximately or exactly. Among them Kahan’s method, also known as Hirota-Kimura method, applied to ODEs with quadratic vector fields, preserves the integrals or conserved quantities of many Hamiltonian and integrable systems [7, 8, 9, 10, 11, 12]. It was introduced by W. Kahan as ”unconventional” discretization method [5] for quadratic vector fields and applied to a scalar Riccati equation and a two-dimensional Lotka-Volterra system [13].

The Lotka-Volterra systems (LVSs) are first order ODEs, initially designed as an ecological predator-prey model. They occur in epidemiology, in laser physics [14], in evolutionary game theory [15] and as spatial discretizations of the Korteweg de Vries equation [13, 16]. Most of the two and three dimensional LVSs have periodic solutions and posses polynomial and non-polynomial integrals. They can be written as Poisson systems in bi-Hamiltonian form [17] and Nambu systems [18]. Many numerical methods are applied to LVSs which preserve the integrals, periodic solution, attractors and son on [6, 19, 20, 21, 22, 23].

For Hamiltonian systems, higher order accuracy for integrals can be achieved by composing symplectic integrators in time [24, 25, 26, 27]. Starting with a basic method, arbitrary orders of accuracy can be obtained by applying the composition to a lower order symplectic method recursively. Another class of numerical methods designed for the construction of high-order approximations to the solution of differential equations are the deferred correction methods. A numerical solution of an initial-value problem (IVP) for ODEs is computed by a low order method and then subsequently refined by solving the IVP constructed by the error between the numerical and continuous solutions. Under suitable assumptions, this process can be repeated to produce solutions with an arbitrarily high order of accuracy. Deferred correction methods have been extensively applied to IVPs such as, classical deferred correction (CDC) methods [28, 29], spectral deferred correction methods [28, 30] and integral deferred correction methods [31].

The LVSs are ODEs with quadratic vector fields. In this paper, two three dimensional (3D) LVSs in bi-Hamiltonian form are solved by the CDC method based on Kahan’s method. We show that the periodicity of the solutions and integrals are preserved in long term integration. At each correction step, more accurate solutions are obtained and the integrals are preserved more accurately. Iterated deferred correction methods are more efficient than the composition methods, because at the correction step the same grid is used. Therefore substantial speedups can be obtained by the CDC methods over the basic method, i.e. Kahan’s method. To the best of our knowledge, the deferred correction methods are used first time to preserve the conserved quantities of dynamical systems with higher accuracy.

The paper is organized as follows. In Section 2, we present two 3D LVSs in bi-Hamiltonian form. In Section 3, we give a short description of Kahan’s method applied to ODEs with a quadratic vector field. Algorithm for the CDC methods is discussed and given in Section 4. Numerical results for Kahan’s method with CDC methods and composition methods are compared in Section 5. The paper ends with some conclusions in Section 6.

2 Lotka-Volterra systems

The LVSs are systems of first order ODEs in the following form [32, 33]

[TABLE]

where ${\mathbf{u}}:=(u_{1},\ldots u_{m})^{T}$ is the $m$-dimensional state vector and $\dot{u}_{i}=du_{i}/dt$ denotes the derivative with respect to time. In ecology, $u_{i}$ describe the densities of each species and $r_{i}$ are the intrinsic growth or decay rates. The interaction between the species is specified by the coefficient matrix $A=(\alpha_{ij}),\;i,j=1,\ldots,m$. All variables in (1) are real and the densities $u_{i}$ are positive. There are no closed solutions of LVSs when $m\geq 2$, they have to be integrated numerically.

2.1 Bi-Hamiltonian 3D Lotka-Volterra systems

Many 3D LVSs can be written in the following bi-Hamiltonian form

[TABLE]

where $J_{1}$ and $J_{2}$ are the skew-symmetric Poisson matrices satisfying the Jacobi identity. There exists two independent integrals $H_{1}$ and $H_{2}$, associated with $J_{1}$ and $J_{2}$ such that $H_{2}$ is the Casimir for one Poisson structure while $H_{1}$ is the Casimir for the other [34]. Bi-Hamiltonian systems are completely integrable [17]. 3D LVSs can also be written as Nambu systems [18], as generalization of Hamiltonian systems with multiple Hamiltonians. Nambu form of (2) is given as

[TABLE]

Vector fields of Nambu systems are divergence free and the flow is volume preserving.

A well known 3D LVS possessing bi-Hamiltonian structure is [35, 21, 36]

[TABLE]

where $\lambda,\mu,\nu>0$, and with $abc=-1$ and $\nu=\mu b-\lambda ab$. The skew-symmetric Poisson matrices and Hamiltonians then are given by

[TABLE]

[TABLE]

[TABLE]

$H_{1}$ and $H_{2}$ are Casimirs of $J_{1}$ and $J_{2}$, respectively, i.e. $J_{1}\nabla H_{1}=0$ and $J_{2}\nabla H_{2}=0$.

Another example of 3D LVS is the reversible 3D LVS with the circulant coefficient matrix $A$ [37]

[TABLE]

It has a game-theoretical interpretation [15] and possesses bi-Hamiltonian form with the Poisson matrices

[TABLE]

and with the linear Hamiltonian $H_{1}=u_{1}+u_{2}+u_{3}$ and the cubic Hamiltonian $H_{2}=u_{1}u_{2}u_{3}$. It is reversible with respect to $\rho={\mathrm{d}iag}(-1,-1,-1)$, $\rho f(\mathbf{x})=-f(\rho{\mathbf{x}})$. It can also be written as Nambu system. The flow generated by (4) is source free, i.e. the volume is preserved. The linear integral $H_{1}$ represents the volume. The $n$-dimensional extension of (4) as integrable discretization of the Korteweg de Vries equation was integrated with a Poisson structure preserving integrator in [16]. Necessary and sufficient conditions for conservation laws of $n$-dimensional LVSs (1) including the two Poisson systems are derived in [33].

3 Kahan’s method

The LVS (1) is an autonomous ODE system in the following form

[TABLE]

with the quadratic vector field $(Q({\mathbf{u}}))_{i}=u_{i}\sum_{j}\alpha_{ij}u_{j}$ and the diagonal matrix $B=diag(r_{1},\ldots,r_{m})$. In the system (5), the unknown solution vector is ${\mathbf{u}}=(u_{1},\ldots,u_{m})^{T}$, and it is prescribed the vector of initial conditions ${\mathbf{u}}(t_{0})=(u_{1}(t_{0}),\ldots,u_{m}(t_{0}))^{T}$.

For the ODE system (5), Kahan introduced in 1993 the ”unconventional” discretization as [5]

[TABLE]

where $\Delta t$ is the step size of the integration, ${\mathbf{u}}_{n+1}$ and ${\mathbf{u}}_{n}$ are the approximations at the time instances $t_{n+1}$ and $t_{n}$, respectively. The symmetric bilinear form $Q(\cdot,\cdot)$ is obtained by the polarization of the quadratic vector field $Q(\cdot)$ [38]

[TABLE]

The Kahan’s method is second order and time-reversal [13]:

[TABLE]

where $I\in\mathbb{R}^{m\times m}$ is the identity matrix and $f^{\prime}$ denotes the Jacobian of $f$. Moreover, Kahan’s method is linearly implicit and it coincides with a certain Rosenbrock method on quadratic vector fields, i.e. ${\mathbf{u}}_{n+1}$ can be computed by solving a single linear system

[TABLE]

Symplectic integrators like the implicit mid-point rule [1], energy preserving average vector field method [22] and conservative methods [39] require at each time step more than one Newton iteration to solve nonlinear implicit equations to preserve the integrals accurately. Due to the linearly implicit nature, Kahan’s method is a very efficient structure preserving integrator for ODEs with quadratic vector fields.

Kahan’s method is also a Runge-Kutta method, with negative weights, restricted to quadratic vector fields [40]:

[TABLE]

Kahan’s method was independently rediscovered by Hirota and Kimura [41, 9], which preserves the integrability for a large number of integrable quadratic vector fields like Euler top, Lagrange top [9, 10, 12] Suslov and Ishii systems, Nambu systems, Riccati systems, and the first Painlevé equation [40, 38]. Kahan’s method is generalized in [38] to cubic and higher degree polynomial vector fields.

4 Iterative deferred correction method

In this section, we apply the CDC method [28, 29] to the ODE system (5) with quadratic vector field, related to the LVSs (1). For a given time interval $[0,T]$, we subdivide it into $J$ equidistant intervals $[t_{j},t_{j+1}]$, $j=0,1,\ldots,J-1$:

[TABLE]

on which we define the approximate solutions by ${\mathbf{u}}_{j}\approx{\mathbf{u}}(t_{j})$, $j=1,\ldots,J$, and for $j=0$ we use the initial condition, ${\mathbf{u}}_{0}={\mathbf{u}}(0)$. Further, each interval $[t_{j},t_{j+1}]$ is subdivided into $n-1$ equidistant intervals forming $n$ nodes including the end points $t_{j}$ and $t_{j+1}$ as

[TABLE]

and we define the approximate solutions on these nodes by ${\mathbf{u}}_{j,i}\approx{\mathbf{u}}(t_{j,i})$, $i=1,\ldots,n$, given that ${\mathbf{u}}_{j,1}={\mathbf{u}}_{j}$. The CDC method, on each subinterval $[t_{j},t_{j+1}]$, starts by solving the ODE system (5) for the solutions at the nodes $\{t_{j,i}\}_{i=1}^{n}$, with a method of order $p_{0}$. Then, the approximate solutions of the ODE system (5) on the interval $[t_{j},t_{j+1}]$ are defined by ${\mathbf{U}}_{j}^{[0]}:=({\mathbf{u}}^{[0]}_{j,1},\ldots,{\mathbf{u}}^{[0]}_{j,n})$, and satisfy that

[TABLE]

We apply here the second order Kahan’s method described in Section 3 as the basic method, so $p_{0}=2$. After, it follows the correction procedure. At the $s$-th correction step, the CDC method computes an improved (corrected) solution ${\mathbf{U}}_{j}^{[s]}:=({\mathbf{u}}^{[s]}_{j,1},\ldots,{\mathbf{u}}^{[s]}_{j,n})$ of the following error system

[TABLE]

by a method of order $p_{s}$. In (6), $\mathbf{e}_{j}^{[s]}(t,\mathbf{U}_{j}^{[s]}(t))$ denotes the error function on the $s$-th iteration step given by

[TABLE]

The differences between the deferred correction methods are based on the formation of an error system; on the continuous level they are equivalent. In the CDC method, the error system (6) is used, which is obtained by the differentiation of the error equation (7) with respect to the time variable $t$. The function ${\mathbf{U}}_{j}^{[s]}(t)$ stands for the continuous approximation of the discrete solutions ${\mathbf{U}}_{j}^{[s]}=({\mathbf{u}}^{[s]}_{j,1},\ldots,{\mathbf{u}}^{[s]}_{j,n})$. Here, we construct the continuous approximation ${\mathbf{U}}_{j}^{[s]}(t)$ based on the Lagrange interpolation as

[TABLE]

where $l_{k}(t)$ are the Lagrange basis functions.

The error system (6) is non-autonomous due the occurrence of the time dependent terms $\mathbf{U}_{j}^{[s-1]}(t)$ and their derivatives. Because Kahan’s method is designed for autonomous systems, in the correction steps we use the second order mid-point rule.

After, defining the vector of error approximations ${\mathbf{E}}_{j}^{[s]}:=({\mathbf{e}}^{[s]}_{j,1},\ldots,{\mathbf{e}}^{[s]}_{j,n})$ where ${\mathbf{e}}^{[s]}_{j,i}$ are the discrete solutions of the error system (6) on the nodes $\{t_{j,i}\}$, we obtain the corrected numerical approximations through the update formula

[TABLE]

An outline of the CDC method can be found in Algorithm 1.

Expected order of accuracy of the CDC methods for uniformly spaced nodes is given by $\min\{P_{S},n-1\}$, where $P_{S}=\sum_{s=0}^{S}p_{s}$, $S$ is the number of corrections and $n$ is the number of nodes used in each interval $[t_{j},t_{j+1}]$ [28, 29]. Since we use Kahan’s method and mid-point method, both of which are second order methods, we have $p_{s}=2$ for all $s=0,1,\ldots,S$, and then the expected order of accuracy becomes $\min\{2S+2,n-1\}$. According to this fact, we set $n=2S+3$ in the simulations to obtain the expected order of accuracy as $n-1$. When non-uniform nodes like Gauss–Lobatto, Gauss–Legendre, and Chebyshev nodes are used, for CDC methods the accuracy improves with more corrections although the order of accuracy stagnates at two [29]. When a low order Lagrange interpolation is used on small intervals $[t_{j},t_{j+1}]$, the CDC method can produce accurate results on uniform nodes, as it will be shown in the next Section.

5 Numerical results

In this section, we present numerical results for the Lotka-Volterra systems described in Section 2 solving by Kahan’s method, and demonstrate the performance of the CDC method. In all examples, we give the results of the runs using Kahan’s method with a small time step-size without CDC method, and the ones with the CDC method using Kahan’s method for the ODE system (5) and the mid-point rule for the error system (6), with a larger time step-size. Hamiltonian errors $H(0)-H(t)$ are plotted over $t$. We set $S=1$ and accordingly $n=2S+3=5$ in the CDC procedure.

The $L^{2}$-error for a Hamiltonian $H(t)$, and the $L^{2}$-error between the exact solution $\mathbf{u}_{exact}(t)$ and the numerical solution $\mathbf{u}$ are measured using the following norms

[TABLE]

where the exact solution $\mathbf{u}_{exact}(t)$ is obtained by MatLab’s ode45 solver in which we set the relative and the absolute tolerances as $10^{-13}$. The order of accuracy is calculated as

[TABLE]

where $\text{Err}_{\Delta t_{i}}$ and $\text{Err}_{\Delta t_{i+1}}$ stand for the $L^{2}$-error of an Hamiltonian or the solution, obtained by the consecutive step sizes $\Delta t_{i}$ and $\Delta t_{i+1}=\Delta t_{i}/2$, respectively.

5.1 Bi-Hamiltonian LVS

We consider the 3D LVS (3) on the interval $[0,100]$, with the parameter values $(a,b,c,\lambda,\mu,\nu)=(-1,-1,-1,0,1,-1)$ [21]. The initial condition is taken as $(u_{1}(0),u_{2}(0),u_{3}(0))^{T}=(1,1.9,0.5)^{T}$.

We show that Kahan’s method preserves the periodicity of the solutions and the Hamiltonians in Figure 1. It was proved in [19] that Kahan’s method preserves the periodicity of LVSs (1). The average vector field method, which preserves the Poisson structure, was also applied to LVSs (1) in [22]. It was shown there that the first Hamiltonian $H_{1}$ of (3) is preserved, but the Casimir $H_{2}$ shows a drift in long term integration.

When the CDC method is applied with $S=1$, $n=2S+3=5$ and with the use of larger time step-size $\Delta t=0.01$, the periodicity of the solutions and Hamiltonians $H_{1}$ and $H_{2}$ are preserved in Figure 2. Compared with the numerical results in Figure 1, it turns out that the use of CDC method is more efficient in terms of the preservation of the Hamiltonians. In Figure 2, a slow drift in the preservation of the Hamiltonians is observed. When composition methods are applied to Kahan’s method, a comparatively more rapid Hamiltonian drift is observed [40].

In Figure 3, we give the $L^{2}$-errors and convergence orders of the solutions, Hamiltonians $H_{1}$ and $H_{2}$ for different correction number $S$ related with the number of nodes $n=2S+3$. When the errors reach about $10^{-10}$ level, the computations are stopped. With increasing number of correction step $S$, larger time steps are used to attain a prescribed order, which demonstrates the computational efficiency of CDC methods.

For varying correction number $S$ and node number $n$, the convergence orders are presented in Table 1 for the preservation of the Hamiltonian $H_{1}$. The results for the Hamiltonian $H_{2}$ and solutions are similar. The red labeled orders in Table 1 correspond to the setting $n=2S+3$, and they agree with the expected orders of accuracy. The time-step size $\Delta t$ in the table are chosen in order to reach the target error level $10^{-10}$ related to the runs with the red labeled orders. The efficiency with respect to the step-size can be seen by the speedup factors in the last column, which are calculated as the ratio of the Wall Clock time required for the run without CDC method over the Wall Clock time required for the one with the CDC method, on the same level of accuracy $10^{-10}$.

5.2 Reversible LVS

We consider the reversible LVS (4) on the interval $[0,100]$, and with the initial conditions $(u_{1}(0),u_{2}(0),u_{3}(0))^{T}=(0.3,0.3,0.4)^{T}$ [21]. Kahan’s method preserves again the periodicity of the reversible LVS (4), Figure 4, top. The reversible LVS (4) was solved in [39] using a conservative multiplier method. It was shown that the linear Hamiltonian $H_{1}$ is preserved with an accuracy $10^{-15}$ and the cubic Hamiltonian $H_{2}$ with an accuracy $10^{-14}$. Kahan’s method also preserves the linear Hamiltonian $H_{1}$ and cubic Hamiltonian $H_{2}$ accurately in Figure 4, bottom.

Preservation of the periodicity of the solutions and the Hamiltonians in the case of CDC method is similar to the LVS (3) in the previous example. Again, similar convergence orders are attained in Figure 5.

6 Conclusions

We have shown that the Hamiltonians of 3D LVSs can be preserved with a high accuracy, when we use CDC methods based on the Kahan’s discretization for quadratic vector fields. In a future work, the integral and spectral correction methods on non-uniform grids will be applied.