Computational enhancement of discrete gradient method

1 Introduction

Dicrete gradient method for ODEs [23] was developed for carying numeric calculations in which essential part was played by preserving a first integral or Lyapunov-type behavior of some quantity. Typically it was used on Hamiltonian systems

\begin{matrix} {˙ q}^{i} = \frac{\partial H}{\partial p_{i}}, {˙ p}_{i} = - \frac{\partial H}{\partial q^{i}}, \end{matrix}

(1.1)

where function $H (q, p)$ was the one to be preserved by the algorithm.

Of course, Hamiltonian systemss, due to their analytic properties [5], are quaitatively well discretized by the symplectic schemes (e.g. textbook [7], classical papers[13, 24], and for more modern geometric approach [17]) even when they admit some perturbation [15]. However, it would be useful to have at our disposal gradient scheme that would be efficient in general, non-conservative case, what this paper aims to deliver.

To convey the main principle on which the method is based, we concentrate exclusively on one-dimensional systems. Then the scheme looks as follows

\begin{matrix} q_{i + 1} = q_{i} + h \frac{H (q_{i + 1}, p_{i + 1}) - H (q_{i + 1}, p_{i})}{p_{i + 1} - p_{i}}, p_{i + 1} = p_{i} - h \frac{H (q_{i + 1}, p_{i}) - H (q_{i}, p_{i})}{q_{i + 1} - q_{i}}, \end{matrix}

(1.2)

where we choose one of possible discrete gradient forms [20], this one called coordinate increment discrete gradient (also one of its versions - may even be symmetric!).

Slightly rearranging and cross-multiplying above equations shows

H (q_{i + 1}, p_{i + 1}) = H (q_{i}, p_{i})

(1.3)

that is, Hamiltonian will be conserved at every step up to machine accuracy.

This brief article is organized as follows: In sec. 2 we use the test-bench example of damped harmonic oscillator, as one with accessible exact solution to measure parameters (local and global error, energy decrementation rate etc,) in our new scheme, then, in sec. 3 we use modified gradient method on van der Pol oscillator - system well described in literature in many different contexts (e.g. [22, 9] and [14] for FitzHugh-Nagumo generalization), of course here we perform only some basic tests. Sec. 4 is the last in the paper and is devoted to concluding remarks and possibilities the new method offers.

2 Modified discrete gradient scheme for non-conservative systems

Now we consider a little more general system

\begin{matrix} ˙ q = \frac{\partial H}{\partial p}, ˙ p = - \frac{\partial H}{\partial q} - D (q, p) . \end{matrix}

(2.1)

An obvious remark is that energy is non-constant for this system, actually it obeys

\frac{d H}{d t} = - p D (q, p),

(2.2)

hence if we would find quantity with exactly opposite time behavior, sum of the energy and this new quantity would be conserved.

Let us introduce new artificial reservoir variable (which, from point of view of solving (2.1), is completely redundant (however numerically situation is completely defferent)

w (q, p, q_{0}, p_{0}) = \int_{q_{0}}^{q} D (q, p) d q = \int_{t_{0}}^{t} D (q, p) p d t,

(2.3)

from what we see their possess proper time derivative (here there matters only numerical values at each moment). The quantity

K = H + w

(2.4)

will be referred to as effectively conserved, as numerically it should stay constant, in other words

˙ K = ˙ H + ˙ w = q ˙ q + p ˙ p + p D (q, p) = q p - q p - p D (q, p) + p D (q, p) = 0

(2.5)

as follows from (2.1).

Strikingly simple physical interpretation of $K$ (initial energy), makes it pretty natural guess for conserved quantity associated with dissipative system. This simple effect is obtained by considering dissipative forces as contained in the system.

Coming to an example of damped oscillator

\begin{matrix} ˙ q = p, ˙ p = - q - b p, ˙ w = b p^{2} \end{matrix}

(2.6)

where $b$ is the damping constant, we put

\begin{matrix} \frac{q_{i + 1} - q_{i}}{h} = \frac{K (q_{i + 1}, p_{i + 1}, w_{i + 1}) - K (q_{i + 1}, p_{i}, w_{i + 1})}{p_{i + 1} - p_{i}}, \frac{p_{i + 1} - p_{i}}{h} = \frac{K (q_{i}, p_{i}, w_{i}) - K (q_{i + 1}, p_{i}, w_{i + 1})}{q_{i + 1} - q_{i}}, \frac{w_{i + 1} - w_{i}}{h} = \frac{1}{2} D (q_{i}, q_{i + 1}, p_{i}, p_{i + 1}) (p_{i} + p_{i + 1}) \end{matrix}

(2.7)

where we are able to express evolution of $w$ in quite arbitrary way - its flow should only obey the condition of becoming $D (q, p) p$ in the continuous limit. Note that we change $w$ every time we change $q$ .

Above scheme guarantees that

K (q_{i}, p_{i}, w_{i}) = K (q_{i + 1}, p_{i + 1}, w_{i + 1})

(2.8)

as declared before. Note that in this way we provided global effective invariant for numerical scheme, hence our reasoning is of huge practical importance.

In the following, we will solve implicit equations with tolerance $ε = 10^{- 16}$ , time-step will take the value $h = 0, 001$ . The system on which we will test the new approach will be damped harmonic oscillator. For reservoir variable we will use method

w_{i + 1} = w_{i} + b \frac{p_{i} + p_{i + 1}}{2} (q_{i + 1} - q_{i}) .

(2.9)

In order to check how our new approach works, we perform numerical experiment, consisting in executing few different algorithms on the same set of initial data. We compare our modification of discrete gradient method (modDG, blue) with symplectic leap-frog scheme (pqpLF, red) and the implicit midpoint rule (IMR, magenta). Note that all schemes used here are of order two.

We use initial conditions $q_{0} = 1.3, p_{0} = - 2.2, w_{0} = 0.0$ .

Continuous system is

\begin{matrix} ˙ q = p, ˙ p = - q - b p, ˙ w = b p^{2}, \end{matrix}

(2.10)

where $b = 0.1$ is the damping constant and we have already included reservoir in the description.

The leap-frog scheme as indicating by equations of motion, is

\begin{matrix} p_{i + \frac{1}{2}} = \frac{p_{i} - \frac{h}{2} q_{i}}{1 + \frac{h b}{2}}, q_{i + 1} = q_{i} + h p_{i + \frac{1}{2}}, p_{i + 1} = p_{i + \frac{1}{2}} - \frac{h}{2} (q_{i + 1} + b p_{i + \frac{1}{2}}) \end{matrix}

(2.11)

so it is an explicit scheme (because equations of motion are linear).

To compare with modGR scheme, we adjoin the reservoir obeying numerical evolution

w_{i + 1} = w_{i} + b p_{i + \frac{1}{2}} (q_{i + 1} - q_{i}) .

(2.12)

Despite the classical features (order, errors behavior, etc.) we will care about the energy decrement at each step

R_{i + 1} = \frac{H (q_{i}, p_{i})}{H (q_{i + 1}, p_{i + 1})} = \frac{K_{0} - w_{i}}{K_{0} - w_{i + 1}},

(2.13)

where we made use of definition of the effectively conserved quantity.

It occurs modGR scheme is of order two for $q, p$ variables, but fortunately of order three for $w$ . hence it is wiser to use this variable, whenever possible. We note that for value of $D$ we used midpoint approximation.

Combination of previously stated remarks leads to fig. 4.2. IMR and pqpLF perform equally poor when compared with our modGR scheme.

Figure 4.1.: Deviation from initial value of $K$ . We adjoined reservoir to pqpLF to see, how it conserves $K$ .

Figure 4.2.: Comparison of deviation from theoretical value of energy decrement. modGR is better through six orders of magnitude.

Figures below center around behavior of $d_{K}$ and $d_{R}$ modulus maxima for different vlues of $b$ . We are not able to perform tests for larger $b$ values since the damped system finally comes to rest, and the bigger $b$ is, the faster it comes. If numerical scheme admits this kind of behavior, when continuous system does, we call corresponding numerical method A-stable. In the next section we will be able to perform more tests since van der Pol oscillator is alternately damped and stimulated.

Below, continuous lines with circles denote $h = 0.001$ , dashed lines with circles are for $h = 0.01$ and continuous lines with asterisks is $h = 0.1$ .

We use also $b = {0.001, 0.01, 0.1, 0.2, 0.3, 0.5, 0.7, 0.9, 1.0, 1.3, 1.5, 1.8}$ .

Figure 4.3.: Deviation from initial value of $K$ . Surprising is the fact, that (in reasonable interval) gradient performs better with bigger time step, but LF goes worse.

Figure 4.4.: $d_{R}$ for different schemes with different time steps. Nearly unprecedented winner is the modGR scheme.

3 Practical test: van der Pol oscillator

Van der Pol oscillator is a widely used physical systems in many contexts (e.g. [8, 9]). Since it is non-integrable, we are forced to study only its numerical properties and could examine how our method work in this particular case. Equations of motion are

\begin{matrix} ˙ q = p, ˙ p = - q + a (1 - q^{2}) p, \end{matrix}

(3.1)

which means for $a << 1$ it is a perturbative system (see e.g. [15] , especially in numerical context).

Modified discrete gradient works as usual, but we choose midpoint approximation to define evolution of artificial variable $w$ :

w_{i + 1} = w_{i} - a (1 - (\frac{q_{i} + q_{i + 1}}{2})^{2}) \frac{p_{i} + p_{i + 1}}{2} (q_{i + 1} - q_{i}) .

(3.2)

Figure 4.4.: On the left $a = 0.05$ showing perturbative behavior and on the right picture $a = 50$ being sstrongly non-linear. All the phase-trajectories admit stable limit cycles.

Proper phase curves behavior for a certain range of values of $a$ gives basic confirmation of our idea being reasonable. Next trial is to check how initial value of $K$ is conserved, when the dynamics is generated by van der Pol term with all the sifferent $a$ ’s.

We use $a = {0.05, 0.1, 0.3, 0.73, 1.0, 2.0, 4.9, 9.8, 16.0, 30.0, 50.0, 60.0, 79.0, 100.0, 140.0}$ and for every value we search for maximum modulus of deviation from initial value of $K$ . Fig. 4.5 is qualitative in nature, since it illustrates only the search during $T = 100.0$ units of time evolution. It does not follow from it that for huge $a$ simulation does not break down for longer times - usually cumulation of errors has this effect. One obvious remark we get from this experiment is that the bigger $a$ we consider, the smaller time step $h$ is needed for safe usage of gradient scheme.

Figure 4.5.: Deviations $d_{K}$ for different values of $a$ . Blue line represents simulation with $h = 0.001$ , red- $h = 0.01$ , green- $h = 10^{- 4}$ and magenta- $h = 10^{- 5}$ .

Effectively conserved quantity increases stability of the scheme (as any conserved quantity always does). In the former case of damped harmonic oscillator it guaranteed conditional A-stability (A-stability for some specified range of time-steps). As shown on the right, the same effect occurs for van der Pol oscillator, e.g. line for $h = 0.01$ (red) completely fails at the values of $a > 50$ , where $K$ iss not being kepr constant properly.

This simple experiment ends our brief survey of finding dissipative analogue of discrete gradient scheme. Nonlinear systems, as stiff by nature, will certainly benefit from such improved numerical treatment.

4 Concluding remarks

Goal of this paper was to show that modified gradient scheme with computational invariant works, and it was achieved. Of course further investigation has to be proceeded concerning technical details of implementation. Also higher-dimensional systems need to be inspected under the treatment, and non-atonomous systems are to be given a closer look (not speaking obout other interesting autonomous cases).

Scheme we proposed is of order two, but rising accuracy is not so easy like in the case of conservative systems. Ideas given in [12] has led to rising order for the method without improving accuracy even a bit. This weird effect is closely related to the form of dissipative systems.

Existence of a conserved $K$ quantity opens possibilities for geometric treatment of non-conservative systems in general, along the lines of GNI (e.g. [18, 16]) approach. Especially when combined with higher-dimensions it may be useful to confront our modification of gradient scheme with ideas of local invariants [19] and local exactness [10].

References

[1]
[2]
[3]
[4]
[5] W.I.Arnold: “Mathematical methods of classical mechanics”, PWN 1981 (in Polish),
[6] R.Belousov, F.Berger, A.J.Huspeth: “Volterra-series approach to stochastic nonlinear dynamics: The Duffing oscillator driven by white noise
”, Phys.Rev.E, 99, 042204 (2019),
[7] S.Blanes, F.Casas: “A Concise Introduction to Geometric Numerical Integration”, CRC Press, 2016,
[8] T.Botari, E.D.Leonel: “One-dimensional Fermi accelerator model with moving wall described by a nonlinear van der Pol oscillator”, Phys.Rev.E, 87, 012904 (2013),
[9] Z.Chen, B.Raman, A.Stern: “Structure preserving numerical integrators for Hodgkin-Huxley-type systems”, arXiv: 181100173v1 [math.NA],
[10] J.L.Cieśliński: “Locally exact modifications of numerical schemes”, Computers & Mathematics with Applications, Volume 65, Issue 12, 2013, Pages 1920-1938,
[11] J.L.Cieśliński, B.Ratkiewicz: “Long-time behavior of discretizations of the simple pendulum equation”,J.Phys.A: Math.Theor. 42 (2009), 105204 (29pp), IOP Publishing,
[12] J.L.Cieśliński, B.Ratkiewicz: “Improving the accuracy of the discrete gradient method in the one-dimensional case”, Phys. Rev. E 81, 016704 (2010),
[13] P.J.Channel, C.Scovel: “Symplectic integration of Hamiltonian systems”, Nonlinearity 3(1990), 231-259, IOP Publishing Ltd,
[14] A.Gawlik, S.Skurativskyi, V.Vladimirov: “Solitary waves dynamics governed by the modified FitzHugh-Nagumo equation”, arXiv: 1906.01865v1 [nlin.PS], 2019,
[15] E.Hairer, C.Lubich: “Invariant tori of dissipatively perturbed Hamiltonian systems under symplectic discretization”, Applied Numerical Mathematics 29 (1999), 57-71, Elsevier,
[16] E.Hairer, C.Lubich, G.Wanner: “Geometric Numerical Integration: Structure-preserving algorithms for ordinary differential equations
”, Springer Verlag, 2006,
[17] E.Hairer, C.J.Zbinden: “Conjugate Symplectic B-Series”, AIP Conference Proceedings 2012, DOI: 10.1063/1.4756053
[18] A.Iserles: “A first course in the numerical analysis of differential equations”, Cambridge University Press, 2009,
[19] T.Matsuda, T.Matsuo: “A new geometric integration approach based on local invariants”, JSIAM Letters Vol.5 (2013), pp. 37-40,
[20] R.I.McLachlan, G.R.W.Quispel, N.Robidoux: “Geometric integration sing discrete gradients”, Phil.Trans.R.Soc.Land. A (1999) 357, 1021-1045, The Royal Society Publishing,
[21] E.Süli, D.Meyers: “An Introduction to numerical analysis”, Cambridge University Press, 10th printing, 2014,
[22] U.Parlitz, W.Lauterborn: “Period-doubling cascades and devil’s staircases of the driven van der Pol oscillator”, Phys.Rev.A, vol 36 (3), 1987,
[23] G.R.W.Quispel, H.W.Capel: “Solving ODEs numerically while preserving a first integral”, Phys.Let.A 218 (1996), 223-228, Elsevier,
[24] H.Yoshida: “Recent progress in the theory and applications of symplectic integrators”, Celestial Mechanics and Dynamical Astronomy 56 (1993), 27-43, Kluwer Academic Publishers,

Computational enhancement of discrete gradient method

Authors

Conservative Integrators for Many-body Problems

DG Approach to Large Bending Plate Deformations with Isometry Constraint

Forced extension of GNI techniques to dissipative systems

Combinatorial Gradient Fields for 2D Images with Empirically Convergent Separatrices

A New Cyclic Gradient Method Adapted to Large-Scale Linear Systems

Numerical analysis of the Landau-Lifshitz-Gilbert equation with inertial effects

Bayesian filtering for nonlinear stochastic systems using holonomic gradient method with integral transform

1 Introduction

2 Modified discrete gradient scheme for non-conservative systems

3 Practical test: van der Pol oscillator

4 Concluding remarks

References

Computational enhancement of discrete gradient method

Authors

Conservative Integrators for Many-body Problems

DG Approach to Large Bending Plate Deformations with Isometry Constraint

Forced extension of GNI techniques to dissipative systems

Combinatorial Gradient Fields for 2D Images with Empirically Convergent Separatrices

A New Cyclic Gradient Method Adapted to Large-Scale Linear Systems

Numerical analysis of the Landau-Lifshitz-Gilbert equation with inertial effects

Bayesian filtering for nonlinear stochastic systems using holonomic gradient method with integral transform

This week in AI

1 Introduction

2 Modified discrete gradient scheme for non-conservative systems

3 Practical test: van der Pol oscillator

4 Concluding remarks

References