Discretized Fast-Slow Systems with Canards in Two Dimensions

TL;DR

This paper investigates how a specific discretization scheme, Kahan's method, preserves the delicate canard phenomena in two-dimensional fast-slow systems, ensuring the occurrence of maximal canards similar to continuous systems.

Contribution

It demonstrates that Kahan's discretization preserves the structure of maximal canards in quadratic fast-slow systems, extending continuous-time results to discrete schemes.

Findings

Kahan discretization preserves maximal canards in quadratic systems.

The analysis uses blow-up and Melnikov methods for discrete maps.

Maximal canards occur between attracting and repelling manifolds in the discretized setting.

Abstract

We study the problem of preservation of maximal canards for time discretized fast-slow systems with canard fold points. In order to ensure such preservation, certain favorable structure preserving properties of the discretization scheme are required. Conventional schemes do not possess such properties. We perform a detailed analysis for an unconventional discretization scheme due to Kahan. The analysis uses the blow-up method to deal with the loss of normal hyperbolicity at the canard point. We show that the structure preserving properties of the Kahan discretization for quadratic vector fields imply a similar result as in continuous time, guaranteeing the occurrence of maximal canards between attracting and repelling slow manifolds upon variation of a bifurcation parameter. The proof is based on a Melnikov computation along an invariant separating curve, which organizes the dynamics of the map similarly to the ODE problem.

Equations

$$\displaystyle\begin{array}[]{r@{\;\,=\;\,}l}x^{\prime}&f(x,y,\lambda,\varepsilon),\\ y^{\prime}&\varepsilon g(x,y,\lambda,\varepsilon),\end{array}$$

$$\displaystyle\begin{array}[]{rcrcl}\varepsilon\dfrac{\mathrm{d}x}{\mathrm{d}\tau}&=&\varepsilon\dot{x}&=&f(x,y,\varepsilon),\\ \dfrac{\mathrm{d}y}{\mathrm{d}\tau}&=&\dot{y}&=&g(x,y,\varepsilon),\quad\ x\in\mathbb{R}^{m},\quad y\in\mathbb{R}^{n},\quad 0<\varepsilon\ll 1\,,\end{array}$$

$$\displaystyle\begin{array}[]{r@{\;\,=\;\,}r}x^{\prime}&f(x,y,\varepsilon),\\ y^{\prime}&\varepsilon g(x,y,\varepsilon).\end{array}$$

$$\displaystyle\begin{array}[]{r@{\;\,=\;\,}l}0&f(x,y,0),\\ \dot{y}&g(x,y,0),\end{array}$$

$$\displaystyle\begin{array}[]{r@{\;\,=\;\,}l}x^{\prime}&f(x,y,0),\\ y^{\prime}&0.\end{array}$$

$$S_0=\{(x,y)\in {\mathbb{R}}^{m+n}:f(x,y,0)=0\}.$$

$$\dot{x}=\frac{g(x,\phi (x),0)}{{\phi }^{\prime }(x)}.$$

$$\displaystyle\begin{array}[]{r@{\;\,=\;\,}l}\varepsilon\dot{x}&-y+x^{2},\\ \dot{y}&x,\end{array}$$

$$(x_{0}(\tau),y_{0}(\tau))=\Big{(}\frac{\tau}{2},\frac{\tau^{2}}{4}\Big{)}.$$

$$(x_{0,\varepsilon}(\tau),y_{0,\varepsilon}(\tau))=\Big{(}\frac{\tau}{2},\frac{\tau^{2}}{4}-\frac{\varepsilon}{2}\Big{)}$$

$$S_{\epsilon }=\left\{(x,y)\in {\mathbb{R}}^2:y=x^2-\frac{\epsilon }{2}\right\},$$

$$\displaystyle\begin{array}[]{r@{\;\,=\;\,}l}x^{\prime}&f(x,y,\lambda,\varepsilon),\\ y^{\prime}&\varepsilon g(x,y,\lambda,\varepsilon).\end{array}$$

$$\displaystyle\begin{array}[]{r@{\;\,=\;\,}l}x^{\prime}&-yk_{1}(x,y,\lambda,\varepsilon)+x^{2}k_{2}(x,y,\lambda,\varepsilon)+\varepsilon k_{3}(x,y,\lambda,\varepsilon),\\ y^{\prime}&\varepsilon(xk_{4}(x,y,\lambda,\varepsilon)-\lambda k_{5}(x,y,\lambda,\varepsilon)+yk_{6}(x,y,\lambda,\varepsilon)),\end{array}$$

$$\displaystyle\begin{array}[]{r@{\;\,=\;\,}l}k_{i}(x,y,\lambda,\varepsilon)&1+\mathcal{O}(x,y,\lambda,\varepsilon)\,,\quad i=1,2,4,5,\\ k_{i}(x,y,\lambda,\varepsilon)&\mathcal{O}(x,y,\lambda,\varepsilon)\,,\quad i=3,6.\end{array}$$

$$\displaystyle\begin{array}[]{r@{\;\,=\;\,}l}a_{1}&\frac{\partial k_{3}}{\partial x}(0,0,0,0),\ a_{2}=\frac{\partial k_{1}}{\partial x}(0,0,0,0),\ a_{3}=\frac{\partial k_{2}}{\partial x}(0,0,0,0),\\ a_{4}&\frac{\partial k_{4}}{\partial x}(0,0,0,0),\ a_{5}=k_{6}(0,0,0,0),\end{array}$$

$$C=\frac{1}{8}(4a_1-a_2+3a_3-2a_4+2a_5).$$

$${\lambda }_c(\sqrt{\epsilon })=-C\epsilon +\mathcal{O}({\epsilon }^{3/2}),$$

$$z^{\prime }=f(z)=Q(z)+Bz+c,$$

$$\frac{\tilde{z}-z}{h}=\bar{Q}(z,\tilde{z})+\frac{1}{2}B(z+\tilde{z})+c,$$

$$\bar{Q}(z,\tilde{z})=\frac{1}{2}(Q(z+\tilde{z})-Q(z)-Q(\tilde{z}))$$

$$\tilde{z}=F_{f}(z,h)=z+h\Big{(}{\textnormal{Id}}-\frac{h}{2}\mathrm{D}f(z)\Big{)}^{-1}f(z).$$

$$\displaystyle\begin{array}[]{r@{\;\,=\;\,}l}\varepsilon\dot{x}&-y+x^{2}+\varepsilon a_{1}x-a_{2}xy,\\ \dot{y}&x-\lambda+a_{5}y+a_{4}x^{2},\end{array}$$

$$\displaystyle\begin{array}[]{r@{\;\,=\;\,}l}x^{\prime}&-y+x^{2}+\varepsilon a_{1}x-a_{2}xy,\\ y^{\prime}&\varepsilon(x-\lambda)+\varepsilon a_{5}y+\varepsilon a_{4}x^{2},\end{array}$$

$$\frac{\tilde{z}-z}{h}=-\frac{1}{2}f(z)+2f\left(\frac{z+\tilde{z}}{2}\right)-\frac{1}{2}f(\tilde{z}).$$

$$\begin{array}[]{r@{\;\,=\;\,}l}\dfrac{\varepsilon}{h}(\tilde{x}-x)&-\dfrac{1}{2}(y+\tilde{y})+x\tilde{x}+\dfrac{\varepsilon a_{1}}{2}(x+\tilde{x})-\dfrac{a_{2}}{2}(\tilde{x}y+x\tilde{y}),\\ \dfrac{1}{h}(\tilde{y}-y)&\dfrac{1}{2}(x+\tilde{x})-\lambda+\dfrac{a_{5}}{2}(y+\tilde{y})+a_{4}x\tilde{x}.\end{array}$$

$$S_{0,h}=\big{\{}(x,y)\in\mathbb{R}^{2}:y=\varphi_{0,h}(x)\big{\}}$$

$$\begin{array}[]{r@{\;\,=\;\,}l}\dfrac{\varepsilon}{h}(\tilde{x}-x)&-\dfrac{1}{2}(y+\tilde{y})+x\tilde{x},\\ \dfrac{1}{h}(\tilde{y}-y)&\dfrac{1}{2}(x+\tilde{x}).\end{array}$$

$$\begin{array}[]{r@{\;\,=\;\,}l}0&-\dfrac{1}{2}(\tilde{y}+y)+\tilde{x}x,\\ \dfrac{1}{h}(\tilde{y}-y)&\dfrac{1}{2}(\tilde{x}+x)\end{array}$$

$$S_{0,h}=\Big{\{}(x,y)\in\mathbb{R}^{2}:y=x^{2}-\frac{h^{2}}{8}\Big{\}}.$$

$$S_{\varepsilon,h}=\Big{\{}(x,y)\in\mathbb{R}^{2}:y=x^{2}-\frac{\varepsilon}{2}-\frac{h^{2}}{8}\Big{\}},$$

Full text

Keywords: slow manifolds, invariant manifolds, blow-up method, loss of normal hyperbolicity, discretization, maps, canards.

Mathematics Subject Classification (2010): 34E15, 34E20, 37M99, 37G10, 34C45, 39A99.

1 Introduction

In this paper, we study the effect of the time discretization upon systems of ordinary differential equations (ODEs) which exhibit the phenomenon called “canards”. It takes place, under certain conditions, in singularly perturbed (slow-fast) systems exhibiting fold points. The simplest form of such a system is

[TABLE]

where we interpret $\varepsilon>0$ as a small time scale parameter, separating between the fast variable $x$ and the slow variable $y$. For $\lambda=0$, the origin is assumed to be a non-hyperbolic fold point, possessing an attracting slow manifold and a repelling slow manifold. One says that the system admits a maximal canard if there are trajectories connecting the attracting and the repelling slow manifolds [4, 15, 32]. This is a non-generic phenomenon which only becomes generic upon including an additional parameter $\lambda$, for the region of $\lambda$’s which is exponentially narrow as $\varepsilon\to 0$. This makes the study of maximal canards especially challenging.

Krupa and Szmolyan [30] have analyzed maximal canards for equation (1.3) by using the blow-up method which allows to effectively handle the non-hyperbolic singularity at the origin. The key idea to use the blow-up method [13, 14] for fast-slow systems goes back to Dumortier and Roussarie [15]. They observed that non-hyperbolic singularities can be converted into partially hyperbolic one by means of an insertion of a suitable manifold, e.g. a sphere, at such a singularity. The dynamics on this inserted manifold are partially hyperbolic, and truly hyperbolic in its neighborhood. The dynamics on the manifold are usually analyzed in different charts. See e.g. [35, Chapter 7] for an introduction into this technique. A non-exhaustive list of different applications to planar fast-slow systems includes [7, 6, 8, 25, 30, 34, 36].

The main ingredient for a proof of maximal canards in [30] is the existence of a constant of motion for the dynamics in the rescaling chart in the blown-up space. This constant of motion can be used for a Melnikov method to compute the separation of the attracting and repelling manifold under perturbations, in particular to find relations between parameters $\varepsilon$ and $\lambda$ under which the manifolds intersect, leading to a maximal canard. The role of this constant of motion suggests that, in order to retain the existence of maximal canards, the right choice of the time discretization scheme becomes of a crucial importance. Indeed, one can show that conventional discretization schemes like the Euler method do not preserve maximal canards. The concept of a structure preserving discretization method is necessary. We investigate time discretization of the ODE (1.3) via the Kahan method which has been shown to preserve various integrability attributes in many examples (and known also as Hirota-Kimura method in the context of integrable systems, see e.g. [29, 42]). We apply the blow-up method, which so far has been mainly used for flows, to the discrete time fast-slow dynamical systems induced by the Kahan discretization procedure. We show that these dynamical systems exhibit maximal canards for $\lambda$ and $\varepsilon$ related by a certain a functional relation existing in a region which exponentially narrow with $\varepsilon\to 0$. Thus, we extend to the discrete time context the previously known feature of the continuous time systems, provided an intelligent choice of the discretization scheme. We would like to stress that, despite the similarity of results to the continuous time case, the techniques of the proofs for the discrete time had to be substantially modified. In particular, the arguments based on the conserved quantity cannot be directly transferred into the discrete time context, since the conserved quantities there are only formal (divergent asymptotic series). Thus, it turned out to be necessary to use more general arguments based on the existence of an invariant measure and an invariant separating curve characterized as a singular curve of an invariant measure. We use also a more general version of the Melnikov method, similar to the one presented in [43].

Note that the application of the blow-up method to the discrete-time problem of folded canards is a considerable extension compared to the Euler discretizations for transcritical singularities, as studied in [18]. The folded canard case has specific dynamic structure, as explained above, such that a structure-preserving discretization method is needed, now performing the blow-up for the rational Kahan mapping. Compared to [18], the kind of map is different, the structure in the singular limit is richer, there is an additional parameter $\lambda$, also rescaled in the blow-up, and the type of result, namely a continuation of the critical object along a two-parameter curve, is new.

Based on observations of this paper, the employment of Kahan’s method for a treatment of canards can also be found in [17]. There, the simplest canonical form for folded, pitchfork and transcritical canards is studied and the focus lies on the linearization along trajectories. While it is demonstrated that explicit Runge-Kutta methods cannot provide symmetry of entry-exit relations, the linearization along the Kahan scheme and similar symmetric, A-stable methods are shown to preserve the typical continuous-time behaviour. Hence, the discussion of symmetry and linear stability in [17] supplements the paper at hand; here, we establish the existence and extension of maximal canards along parameter combinations for the nonlinear problem of folded canards with additional quadratic perturbation terms, in particular using the blow-up technique.

The paper is organized as follows. Section 2 recalls the setting of fast-slow systems in continuous time and summarizes the main result on maximal canards, Theorem 2.2, with a short sketch of the proof, as given in [30]. In Section 3, we study the problem of a maximal canard for systems with folds in discrete time. We establish the Kahan discretization of the canard problem in Section 3.1 and discuss the reduced subsystem of the slow time scale in Section 3.2. In Section 3.3, we introduce the blow-up transformation for the discretized problem. We discuss the dynamics for the entering and exiting chart in Section 3.4, and for the rescaling chart in Section 3.5. In Section 3.6, we explore the dynamical properties of the Kahan map in the rescaling chart, including a formal conserved quantity, an invariant measure and an invariant separating curve. Following this, we conduct the Melnikov computation along the invariant curve in Section 3.7, leading to the proof of the main Theorem 3.11, which is the discrete-time analogue to Theorem 2.2. Finally, we provide various numerical illustrations in Section 3.8 and conclude with an outlook in Section 4.

Thus, we succeeded in adding the problem of maximal canards to the recent list of results, where a geometric analysis shows that certain features of fast-slow systems with non-hyperbolic singularities can be preserved via a suitable discretization, including the cases of the fold singularity [39], the transcritical singularity [18] and the pitchfork singularity [1]. More broadly viewed, our results also provide a continuation of a line of research on discrete-time fast-slow dynamical systems, which includes the study of canard/delay behavior in iterated maps via normal form transformations [38], non-standard analysis [20, 21], renormalization [2], Gevrey series [3], complex-analytic methods [22], and phase plane partitioning [37].

Acknowledgments: The authors gratefully acknowledge support by DFG (the Deutsche Forschungsgemeinschaft) via the SFB/TR 109 “Discretization in Geometry and Dynamics”. ME acknowledges support by Germany’s Excellence Strategy – The Berlin Mathematics Research Center MATH+ (EXC-2046/1, project ID: 390685689), and CK acknowledges support by a Lichtenberg Professorship of the VolkswagenFoundation.

2 Maximal canard through a fold in continuous time

2.1 Fast-slow systems

We start with a brief review and notation for continuous-time fast-slow systems. Consider a system of singularly perturbed ordinary differential equations (ODEs) of the form

[TABLE]

where $f,g,$ are $C^{k}$-functions with $k\geq 3$. Since $\varepsilon$ is a small parameter, the variables $x$ and $y$ are often called the fast and the slow variables, respectively. The time variable $\tau$ in (2.3) is termed the slow time scale. The change of variables to the fast time scale $t:=\tau/\varepsilon$ transforms the system (2.3) into ODEs

[TABLE]

To both systems (2.3) and (2.6) there correspond respective limiting problems for $\varepsilon=0$: the reduced problem (or slow subsystem) is given by

[TABLE]

and the layer problem (or fast subsystem) is

[TABLE]

The reduced problem (2.9) can be understood as a dynamical system on the critical manifold

[TABLE]

Observe that the manifold $S_{0}$ consists of equilibria of the layer problem (2.12). $S_{0}$ is called normally hyperbolic if for all $p\in S_{0}$ the matrix $\textnormal{D}_{x}f(p)\in\mathbb{R}^{m\times m}$ has no eigenvalues on the imaginary axis. For a normally hyperbolic $S_{0}$, Fenichel theory [19, 28, 35, 44] implies that, for sufficiently small $\varepsilon$, there is a locally invariant slow manifold $S_{\varepsilon}$ such that the restriction of (2.3) to $S_{\varepsilon}$ is a regular perturbation of the reduced problem (2.9). Furthermore, it follows from Fenichel’s perturbation results that $S_{\varepsilon}$ possesses an invariant stable and unstable foliation, where the dynamics behave as a small perturbation of the layer problem (2.12).

2.2 Main result on maximal canards in slow-fast systems with a fold

A challenging phenomenon is the breakdown of normal hyperbolicity of $S_{0}$ such that Fenichel theory cannot be applied. Typical examples of such a breakdown are found at bifurcation points $p\in S_{0}$, where the Jacobi matrix $\mathrm{D}_{x}f(p)$ has at least one eigenvalue with zero real part. The simplest examples are folds in planar systems ($m=n=1$), i.e., points $p=(x_{0},y_{0})\in\mathbb{R}^{2}$ (without loss of generality $p=(x_{0},y_{0})=(0,0)$) where $\partial f/\partial x$ vanishes and in whose neighbourhood $S_{0}$ looks like a parabola. The left part of $S_{0}$ (with $x<0$) is denoted by $S_{a}$ ($a$ for “attractive”), while its right part (with $x>0$) is denoted by $S_{r}$ ($r$ for “repelling”). These notations refer to the properties of dynamics of the layer problem in the region $y>0$ (see e.g. [35, Figure 8.1]). By standard Fenichel theory, for sufficiently small $\varepsilon>0$, outside of an arbitrarily small neighborhood of $p$, the manifolds $S_{a}$ and $S_{r}$ perturb smoothly to invariant manifolds $S_{a,\varepsilon}$ and $S_{r,\varepsilon}$.

In the following we focus on the particularly challenging problem of fold points admitting maximal canards. In this case, the critical curve $S_{0}=\{f(x,y,0)=0\}$ can be locally parametrized as $y=\varphi(x)$ such that the reduced dynamics on $S_{0}$ are given by

[TABLE]

In our setting, the function at the right-hand side is smooth at the origin, so that the reduced flow goes through the origin via a maximal solution $x_{0}(t)$ of (2.13) with $x_{0}(0)=0$. The solution $(x_{0}(t),y_{0}(t))$ with $y_{0}(t)=\varphi(x_{0}(t))$ connects both parts $S_{a}$ and $S_{r}$ of $S_{0}$. However, there is no reason to expect that for $\varepsilon>0$, the (extension of the) solution parametrizing $S_{a,\varepsilon}$ will coincide with the (extension of the) solution parametrizing $S_{r,\varepsilon}$, unless there are some special reasons, like symmetry, forcing such a coincidence.

Definition 2.1**.**

We say that a planar slow-fast system admits a maximal canard, if the extension of the attracting slow manifold $S_{a,\varepsilon}$ coincides with the extension of a repelling slow manifold $S_{r,\varepsilon}$.

Example. Consider the system

[TABLE]

corresponding to $f(x,y,\varepsilon)=x^{2}-y$ and $g(x,y,\varepsilon)=x$. For the reduced system ($\varepsilon=0$) we obtain $y=\varphi(x)=x^{2}$ and $2x\dot{x}=x$, hence $\dot{x}=1/2$ (regular at $x=0$). The solution $x_{0}(t)$ is given by $x_{0}(t)=\tau/2$ so that

[TABLE]

Observe that the system is symmetric with respect to the reversion of time $\tau\mapsto-\tau$ simultaneously with $x\mapsto-x$. This ensures the existence of the maximal canard also for any $\varepsilon>0$. In this particular example, one can easily find the maximal canard explicitly. Indeed, one can easily check that, for any $\varepsilon>0$,

[TABLE]

is a solution of (2.16) which parametrizes the invariant set

[TABLE]

which consists precisely of the attracting branch $S_{a,\varepsilon}=\left\{(x,y)\in S_{\varepsilon}\,:\,x<0\right\}$ and the repelling branch $S_{r,\varepsilon}=\left\{(x,y)\in S_{\varepsilon}\,:\,x>0\right\}$, such that trajectories on $S_{\varepsilon}$ go through $x=0$ with the speed $\dot{x}=\varepsilon/2$. However, any generic perturbation of this example, e.g. with $g(x,y,\varepsilon)=x+x^{2}$, will destroy its peculiarity and will not display a maximal canard.

Thus, maximal canards are not a generic phenomenon in the above setting. In order to find a context where they become generic, we have to consider families depending on an additional parameter $\lambda$:

[TABLE]

We assume that at $\lambda=\varepsilon=0$, the vector fields $f$ and $g$ satisfy the above conditions. By a local change of coordinates, the problem can be brought into the canonical form

[TABLE]

where

[TABLE]

The main result on existence of maximal canards, as given in [30, Theorem 3.1], can be summarized as follows. Set

[TABLE]

and

[TABLE]

Theorem 2.2**.**

Consider system (2.23) such that the solution $(x_{0}(t),y_{0}(t))$ of the reduced problem for $\varepsilon=0$, $\lambda=0$ connects $S_{a}$ and $S_{r}$. Assume that $C\neq 0$. Then there exist $\varepsilon_{0}>0$ and a smooth function

[TABLE]

defined on $[0,\varepsilon_{0}]$ such that for $\varepsilon\in[0,\varepsilon_{0}]$ there is a maximal canard, that is, the extended attracting slow manifold $S_{a,\varepsilon}$ coincides with the extended repelling slow manifold $S_{r,\varepsilon}$, if and only if $\lambda=\lambda_{c}(\sqrt{\varepsilon})$.

The main result of this paper will be a discretized version of Theorem 2.2 restricted to quadratic vector fields, proving for Kahan maps the extension of canards along a parameter curve, as opposed to [17] where only Example (2.16) and its linearization are studied.

The proof of Theorem 2.2 is based on the blow-up technique, transforming the singular problem to a manifold where the dynamics can be desingularized and studied in two different charts. The crucial step in the second chart $K_{2}$ is the continuation of center manifold connections via a Melnikov method based on an integral of motion $H$. In the Appendix A, we summarize this procedure from [30], adding several observations on the dynamics; its separatrix, its invariant measure and an alternative non-Hamiltonian expression that relates to the discrete-time proof we will provide in the following.

3 Maximal canard for a system with a fold in discrete time

3.1 Kahan discretization of canard problem

We discretize system (2.23) with the Kahan method. It was introduced in [29] as an unconventional discretization scheme applicable to arbitrary ODEs with quadratic vector fields. It was demonstrated in [40, 41, 42] and in [5] that this scheme tends to preserve integrals of motion and invariant volume forms. There are few general results available to support this claim, in particular, two general cases of preservation of invariant volume forms in [41, Section 2] and a similar result for Hamiltonian systems with a cubic Hamilton function in [5]. However, the number of particular results not covered by any general theory and reviewed in the above references, is quite impressive. Our study here will contribute an additional evidence, as the result of Section 3.6.2 also belongs to this category, i.e., is not covered by known general statements.

Consider an ODE with a quadratic vector field:

[TABLE]

where each component of $Q:\mathbb{R}^{n}\to\mathbb{R}^{n}$ is a quadratic form, $B\in\mathbb{R}^{n\times n}$ and $c\in\mathbb{R}^{n}$. The Kahan discretization of this system reads as

[TABLE]

where

[TABLE]

is the symmetric bilinear form such that $\bar{Q}(z,z)=Q(z)$. Note that equation (3.2) is linear with respect to $\tilde{z}$ and therefore defines a rational map $\tilde{z}=F_{f}(z,h)$, which approximates the time $h$ shift along the solutions of the ODE (3.1). Further note that $F_{f}^{-1}(z,h)=F_{f}(z,-h)$ and, hence, the map is birational. An explicit form of the map $F_{f}$ defined by equation (3.2) is given by

[TABLE]

In order to be able to apply the Kahan discretization scheme, we restrict ourselves to systems (2.3), (2.6) which are quadratic, that is, to

[TABLE]

resp.

[TABLE]

which corresponds to normal forms (2.23) with $k_{1}=1+a_{2}x$, $k_{2}=1$, $k_{3}=a_{1}x$, $k_{4}=1+a_{4}x$, $k_{5}=1$, and $k_{6}=a_{5}$.

Remark 3.1*.*

It was demonstrated in [5, Proposition 1] that Kahan map (3.3) coincides with the map produced by the following implicit Runge-Kutta scheme, when the latter is applied to a quadratic vector field $f$:

[TABLE]

This opens the way of extending our present results for more general (not necessarily quadratic) systems (2.23). In the present paper we restrict ourselves to the case (3.9), since the algebraic structure keeps the calculations clear and explicit and demonstrates the central methodological aspects of our proofs. However, we additionally apply the scheme (3.10) to the folded canard problem with cubic nonlinearity in Section 3.8, illustrating its numerical capacity beyond the quadratic case. A proof of maximal canards for the non-quadratic case remains an open problem for future work.

3.2 Reduced subsystem of the slow flow

Kahan discretization of (3.6) reads:

[TABLE]

Proposition 3.2**.**

The reduced system (3.11) with $\varepsilon=0$ defines an evolution on a curve

[TABLE]

which supports a one-parameter family of solutions $x_{h}(n;x_{0})$ with $x_{h}(0;x_{0})=x_{0}$. For small $\varepsilon>0$, this curve is perturbed to normally hyperbolic invariant curves $S_{a,h,\varepsilon}$ resp. $S_{r,h,\varepsilon}$ of the slow flow (3.12) for $x<0$, resp. for $x>0$.

For the simplest case $a_{1}=a_{2}=a_{4}=a_{5}=0$ and $\lambda=0$,

[TABLE]

everything can be done explicitly. Straightforward computations lead to the following results.

The reduced system

[TABLE]

has an invariant critical curve

[TABLE]

The evolution on this curve is given by $\tilde{x}=x+\frac{h}{2}$, so that $x_{h}(n;x_{0})=x_{0}+\frac{nh}{2}$.

For the full system (3.12), the symmetry $x\mapsto-x$, $h\to-h$ ensures the existence of an invariant curve

[TABLE]

whose parts with $x<0$, resp $x>0$ are the invariant curves $S_{a,h,\varepsilon}$ resp. $S_{r,h,\varepsilon}$. This curve supports solutions with $x(n)=x_{0}+\frac{nh}{2}$. Thus, system (3.12) exhibits a maximal canard. Our goal is to establish the existence of a maximal canard for system (3.11).

3.3 Blow-up of the fast flow

Kahan discretization of the fast flow (3.9) is the system (3.11) with $h\mapsto h\varepsilon$:

[TABLE]

We introduce a quasi-homogeneous blow-up transformation for the discrete time system, interpreting the step size $h$ as a variable in the full system. Similarly to the continuous time situation, the transformation reads

[TABLE]

where $(\bar{x},\bar{y},\bar{\varepsilon},\bar{\lambda},r,\bar{h})\in B:=S^{2}\times[-\kappa,\kappa]\times[0,\rho]\times[0,h_{0}]$ for some $h_{0},\rho,\kappa>0$. The change of variables in $h$ is chosen such that the map is desingularized in the relevant charts.

This transformation is a map $\Phi:B\to\mathbb{R}^{5}$. If $F$ denotes the map obtained from the time-discretization, the map $\Phi$ induces a map $\overline{F}$ on $B$ by $\Phi\circ\overline{F}\circ\Phi^{-1}=F$. Analogously to the continuous time case, we are using the charts $K_{i}$, $i=1,2$, to describe the dynamics. The chart $K_{1}$ (setting $\bar{y}=1$) focuses on the entry and exit of trajectories, and is given by

[TABLE]

In the scaling chart $K_{2}$ (setting $\bar{\varepsilon}=1$) the dynamics arbitrarily close to the origin are analyzed. It is given via the mapping

[TABLE]

The change of coordinates from $K_{1}$ to $K_{2}$ is denoted by $\kappa_{12}$ and, for $\varepsilon_{1}>0$, is given by

[TABLE]

Similarly, for $y>0$, the map $\kappa_{21}=\kappa_{12}^{-1}$ is given by

[TABLE]

3.4 Dynamics in the entering and exiting chart $K_{1}$

Here we extend the dynamical equations (3.16) by

[TABLE]

and then introduce the coordinate chart $K_{1}$ by (3.17):

[TABLE]

defined on the domain

[TABLE]

where $\rho,\delta,\nu>0$ are sufficiently small.

To transform the map (3.16) into the coordinates of $K_{1}$, we start with the particular case $a_{1}=a_{2}=a_{4}=a_{5}=0$, generated by difference equations

[TABLE]

supplied, as usual, by (3.21). Written explicitly, this is the map

[TABLE]

where

[TABLE]

Upon substitution $K_{1}$, we have:

[TABLE]

where

[TABLE]

Setting

[TABLE]

[TABLE]

we come to the following expression for the map (3.25) in the chart $K_{1}$:

[TABLE]

Now it is straightforward to extend these results to the general case of the map (3.16) with arbitrary constants $a_{i}$. For this, we observe:

• –

in the first equation, the terms $y$ and $x^{2}$ on the right-hand side scale as $r_{1}^{2}$ and $r_{1}^{2}x_{1}^{2}$, while the terms $\varepsilon x$ and $xy$ scale as $r_{1}^{3}\varepsilon_{1}x_{1}$ and $r_{1}^{3}x_{1}$, respectively;

• –

in the second equation, the terms $\varepsilon x$ and $\varepsilon\lambda$ on the right-hand side scale as $r_{1}^{3}\varepsilon_{1}x_{1}$ and $r_{1}^{3}\varepsilon_{1}\lambda_{1}$, while the terms $\varepsilon y$ and $\varepsilon x^{2}$ scale as $r_{1}^{4}\varepsilon_{1}$ and $r_{1}^{4}\varepsilon_{1}x_{1}^{2}$, respectively.

Therefore, we can treat all terms involving $a_{1},a_{2},a_{4},a_{5}$ as $\mathcal{O}(r_{1})$. The resulting map is given by formulas analogous to (3.42), with $X_{1}(x_{1},\varepsilon_{1},\lambda_{1},h_{1})$, $Y_{1}(x_{1},\varepsilon_{1},\lambda_{1},h_{1})$ replaced by certain functions

[TABLE]

We now analyze the dynamics of this map.

• •

The subset $\{r_{1}=0,\;\varepsilon_{1}=0,\;\lambda_{1}=0\}\cap D_{1}$ is invariant, and on this subset we have $Y_{1}(x_{1},r_{1},\varepsilon_{1},\lambda_{1},h_{1})=1$, so that

[TABLE]

Hence, it contains two curves of fixed points

[TABLE]

We have:

[TABLE]

for $h_{1}\leq\nu<1$, hence the point $p_{a,1}(h_{1})$ is attracting in the $x_{1}$-direction and the point $p_{r,1}(h_{1})$ is repelling in the $x_{1}$-direction. In all other directions, the multipliers of these fixed points are equal to 1.

• •

Similarly, we have on $\{\varepsilon_{1}=0,\lambda_{1}=0\}\cap D_{1}$ for small $r_{1}>0$:

[TABLE]

By the implicit function theorem, we can conclude that on $\{\varepsilon_{1}=0,\;\lambda_{1}=0\}\cap D_{1}$, there exist two families of normally hyperbolic (for $h_{1}>0$) curves of fixed points denoted as $S_{a,1}(h_{1})$ and $S_{r,1}(h_{1})$, parametrized by $r_{1}\in[0,\rho]$ and ending for $r_{1}=0$ at $p_{a,1}(h_{1})$ and $p_{r,1}(h_{1})$, respectively. For the map (3.25), corresponding to difference equation (3.24) (that is, to (3.16) with all $a_{i}=0$), the $\mathcal{O}(r_{1})$-term vanishes, and the above families are simply given by

[TABLE]

• •

On the invariant set $\{r_{1}=0,\lambda_{1}=0\}\cap D_{1}$, the dynamics of $x_{1}$, $\varepsilon_{1}$ and $h_{1}$ are given by

[TABLE]

We compute the Jacobi matrices of the map (3.46) at $p_{a,1}(h_{1})$ and $p_{r,1}(h_{1})$, restricting to the invariant set $\{r_{1}=0,\lambda_{1}=0\}\subset D_{1}$,

[TABLE]

The matrix $A_{a}$ has a two-dimensional invariant space corresponding to the eigenvalue 1, spanned by the vectors $v_{a}^{(1)}=(0,0,1)^{\top}$ and $v_{a}^{(2)}=(-1,4,0)^{\top}$, such that

[TABLE]

Similarly, the matrix $A_{r}$ has a two-dimensional invariant space corresponding to the eigenvalue 1, spanned by the vectors $v_{r}^{(1)}=(0,0,1)^{\top}$ and $v_{r}^{(2)}=(1,4,0)^{\top}$, such that

[TABLE]

It is instructive to compare this with the continuous-time case $h_{1}\to 0$ (see, e.g., [30, Lemma 2.5]), where both vectors $v_{a}^{(1)}$ and $v_{a}^{(2)}$ are eigenvectors of the corresponding linearized system, with $v_{a}^{(1)}$ being tangent to $S_{a,1}$ and $v_{a}^{(2)}$ corresponding to the center direction in the invariant plane $r_{1}=0$ (and similarly for $v_{r}^{(1)}$ and $v_{r}^{(2)}$).

We summarize these observations into the following statement.

Proposition 3.3**.**

For system (3.42), there exist a center-stable manifold $\widehat{M}_{a,1}$ and a center-unstable manifold $\widehat{M}_{r,1}$, with the following properties:

For $i=a,r$, the manifold $\widehat{M}_{i,1}$ contains the curve of fixed points $S_{i,1}(h_{1})$ on $\{\varepsilon_{1}=0,\ \lambda_{1}=0\}\subset D_{1}$, parametrized by $r_{1}$, and the center manifold $N_{i,1}$ whose branch for $\varepsilon_{1},h_{1}>0$ is unique (see Figure 3 (b)). In $D_{1}$, the manifold $\widehat{M}_{i,1}$ is given as a graph $x_{1}=\hat{g}_{i}(r_{1},\varepsilon_{1},\lambda_{1},h_{1})$. 2. 2.

For $i=a,r$, there exist two-dimensional invariant manifolds $M_{i,1}$ which are given as graphs $x_{1}=g_{i}(r_{1},\varepsilon_{1})$.

Proof.

The first part follows by standard center manifold theory (see, e.g., [26]). There exist two-dimensional center manifolds $N_{a,1}$ and $N_{r,1}$, parametrized by $h_{1},\varepsilon_{1}$, which at $\varepsilon_{1}=0$ coincide with the sets of fixed points

[TABLE]

respectively (see Figure 3 (b)). Note that, by (3.46), on $\{r_{1}=0,\ \lambda_{1}=0,\ h_{1}>0\}\cap D_{1}$ we have $\tilde{\varepsilon}_{1}>\varepsilon_{1}$ and $\tilde{h}_{1}<h_{1}$ for $x_{1}\leq 0$. Hence, for $\delta$ small enough, the branch of the manifold $N_{a,1}$ on $\{r_{1}=0,\varepsilon_{1}>0,\lambda_{1}=0,h_{1}>0\}\cap D_{1}$ is unique. On the other hand, we observe that for $x_{1}\geq\frac{1}{K}$ with a constant $K>1$, we have $\tilde{\varepsilon}_{1}<\varepsilon_{1}$ and $\tilde{h}_{1}>h_{1}$, if and only if $h_{1}<\frac{2K}{1+K^{2}}$. Thus, for $x_{1}$ from a neighborhood of $1$, we see that $\nu<\frac{2K}{1+K^{2}}<1$ guarantees that, for $\delta$ small enough depending on $K$, the branch of the manifold $N_{r,1}$ on $\{r_{1}=0,\ \varepsilon_{1}>0,\ \lambda_{1}=0,\ h_{1}>0\}\cap D_{1}$ is unique.

The second part follows from the invariances $\tilde{r}_{1}\tilde{\lambda}_{1}=r_{1}\lambda_{1}$ and $\tilde{h}_{1}/\tilde{r}_{1}=h_{1}/r_{1}$, compare [18, Proposition 3.3 and Figure 2] for details. ∎

3.5 Dynamics in the scaling chart $K_{2}$

Next, we investigate the dynamics in the scaling chart $K_{2}$, in order to find a trajectory connecting $\widehat{M}_{a,1}$ with $\widehat{M}_{r,1}$, or $M_{a,1}$ with $M_{r,1}$ respectively. Recall from (3.18) that in chart $K_{2}$ we have

[TABLE]

In this chart and upon the time rescaling $t=t_{2}/r_{2}$, equation (3.9) takes the form

[TABLE]

where the prime now denotes the derivative with respect to $t_{2}$, compare (A.12). Since in this chart $r_{2}=\sqrt{\varepsilon}$ is not a dynamical variable (remains fixed in time), we will not write down explicitly differential, resp. difference evolution equations for $\lambda_{2}=\lambda/\sqrt{\varepsilon}$ and for $h_{2}=h\sqrt{\varepsilon}$. We will restore these variables as we come to the matching with the chart $K_{1}$. The Kahan discretization of equation (3.51) with the time step $h_{2}$ can be written as

[TABLE]

On the blow-up manifold $r_{2}=0$, we are dealing with the simple model system

[TABLE]

This yields the birational map

[TABLE]

This gives the following expressions for the map $F=(F_{1},F_{2})$ and $\hat{J}=(\hat{J}_{1},\hat{J}_{2})$ in (3.54):

[TABLE]

and

[TABLE]

Explicit expressions for the functions $\hat{G}_{1}$ and $\hat{G}_{2}$ can be easily obtained, as well, but are omitted here due to their length.

3.6 Dynamical properties of the model map in the scaling chart

For a better readability, we omit index “2” referring to the chart $K_{2}$ starting from here. In particular, we write $x$, $y$, $r$, $\lambda$, $h$ for $x_{2}$, $y_{2}$, $r_{2}$, $\lambda_{2}$, $h_{2}$ rather than for the original variables (before rescaling). Similarly to the continuous-time case, we start the analysis in $K_{2}$ with the case $\lambda=0$, $r=0$ for $h>0$ fixed. This means that we study the dynamics of the map given by $F$ (3.61),

[TABLE]

which comes as the solution of the difference equation

[TABLE]

We discuss in detail the most important properties of the model map (3.65).

3.6.1 Formal integral of motion

Recall that, for $r=\lambda=0$, the ODE system (A.14) in the chart $K_{2}$ has a conserved quantity (A.15). Its level set $H(x,y)=0$ supports the special canard solution (A.19),

[TABLE]

In general, Kahan discretization has a distinguished property of possessing a conserved quantity for unusually numerous instances of quadratic vector fields. For (A.14), it turns out to possess a formal conserved quantity in the form of an asymptotic power series in $h$. However, there are indications that this power series is divergent, so that map $F$ (3.65) does not possess a true integral of motion. Nevertheless, it possesses all nice properties of symplectic or Poisson integrators, in particular, a truncated formal integral is very well preserved on very long intervals of time. Moreover, as we will now demonstrate, the zero level set of the formal conserved quantity supports the special family of solutions of the discrete time system crucial for our main results.

We recall a method for constructing a formal conserved quantity

[TABLE]

for the the Kahan discretization $F_{f}$ (3.3) for an ODE of the form (3.1) admitting a smooth conserved quantity $H:\mathbb{R}^{n}\to\mathbb{R}$. The latter means that

[TABLE]

The ansatz (3.67) containing only even powers of $h$ is justified by the fact that the Kahan method is a symmetric linear discretization scheme. Writing $\tilde{z}=F_{f}(z,h)$, we formulate our requirement of $\bar{H}$ being an integral of motion for $F_{f}$ as $\bar{H}(z,h)=\bar{H}(\tilde{z},h)$ on $\mathbb{R}^{n}\times[0,h_{0}]$, i.e., up to terms $\mathcal{O}(h^{4})$,

[TABLE]

To compute the Taylor expansion of the left hand side, we observe:

[TABLE]

Here, the $h$ and the $h^{2}$ terms vanish, as follows from (3.68) and its Lie derivative:

[TABLE]

Thus, we find: $H(\tilde{z})=H(z)+\mathcal{O}(h^{3})$, or, more precisely,

[TABLE]

Plugging this, as well as a Taylor expansion of $H_{2}(\tilde{z})$ similar to $H(\tilde{z})$, into (3.69), we see that vanishing of the $h^{3}$ terms is equivalent to

[TABLE]

This is a linear PDE defining $H_{2}$ up to an additive term which is an arbitrary function of $H$.

Following terms $H_{4},H_{6},\dots$ can be determined in a similar manner, from linear PDEs like (3.72) with recursively determined functions on the right hand side.

We now apply this scheme to obtain (the first terms of) the formal conserved quantity $\bar{H}(x,y,h)$ for (3.61). It turns out to be possible to find it in the form

[TABLE]

where

[TABLE]

with $\bar{H}_{2k}(x,y)$ being polynomials of degree $2k+2$. The symbol $\approx$ reminds that this is only a formal asymptotic series which does not converge to a smooth conserved quantity. A Taylor expansion of $H(\tilde{x},\tilde{y})$ as in (3.71) gives

[TABLE]

with

[TABLE]

The differential equation (3.72) reads in the present case:

[TABLE]

A solution for $\bar{H}_{2}$ which is a polynomial of degree 4 reads:

[TABLE]

Hence, we obtain the approximation

[TABLE]

A straightforward computation shows that on the curve $y-x^{2}+\frac{1}{2}=0$ (the level set $H(x,y)=0$), the function $\bar{H}_{2}(x,y)$ takes a constant value $\frac{1}{8}$. Therefore, the level set $\bar{H}(x,y,h)=0$ is given, up to $\mathcal{O}(h^{4})$, by

[TABLE]

Remarkably, we have the following statement.

Proposition 3.4**.**

The curve (3.78) represents a zero level set of the (divergent) formal integral $\bar{H}(x,y,h)$. More precisely, on this curve

[TABLE]

We will not prove this statement, but rather derive a different dynamical characterization of the curve (3.78).

3.6.2 Invariant measure

Proposition 3.5**.**

The map $F$ given by (3.65) admits an invariant measure

[TABLE]

with $\varphi_{h}(x,y)$ given in (3.78). This measure $\mu_{h}$ is singular on the curve $\varphi_{h}(x,y)=0$.

Proof.

Difference equations (3.66) can be written as a linear system for $(\tilde{x},\tilde{y})$:

[TABLE]

Differentiating with respect to $x,y$, we obtain:

[TABLE]

Computing determinants, we find:

[TABLE]

Next, we derive from the first equation in (3.65):

[TABLE]

Since the system (3.66) is symmetric with respect to interchanging $(x,y)\leftrightarrow(\tilde{x},\tilde{y})$ with the simultaneous change $h\mapsto-h$, we can perform this operation in the latter equation, resulting in

[TABLE]

Comparing the last two formulas, we obtain:

[TABLE]

or, equivalently,

[TABLE]

Together with (3.80), this results in

[TABLE]

which is equivalent to the statement of proposition. ∎

3.6.3 Invariant separating curve

It turns out that the singular curve of the invariant measure $\mu_{h}$ is an invariant curve under the map (3.65).

Proposition 3.6**.**

The parabola

[TABLE]

is invariant under the map $F$ given by (3.65). Solutions on $S_{h}$ are given by

[TABLE]

For $(x,y)\in S_{h}$, we have:

[TABLE]

Proof.

Plugging $y=x^{2}-\frac{1}{2}-\frac{h^{2}}{8}$ into formulas (3.65), we obtain upon a straightforward computation:

[TABLE]

This proves the first two claims.

As for the last claim, we compute by differentiating the first equation in (3.65):

[TABLE]

For $(x,y)\in S_{h}$, this gives:

[TABLE]

which implies inequalities (3.85). (We remark that the right hand side tends to infinity as $x\to(1+\frac{h^{2}}{4})/h$.) ∎

The invariant set $S_{h}$ (3.83) plays the role of a separatrix for $F$ (3.61): bounded orbits of $F$ lie above $S_{h}$, while unbounded orbits of $F$ lie below $S_{h}$, as illustrated in Figures 1, 2.

We can show the following connection to the chart $K_{1}$:

Lemma 3.7**.**

The trajectory $\gamma_{h}(n)$, transformed into the chart $K_{1}$ via

[TABLE]

for large $\left|n\right|$, lies in $\widehat{M}_{a,1}$ as well as in $\widehat{M}_{r,1}$.

Proof.

From (3.20) there follows that for sufficiently large $\left|n\right|$, the component $\varepsilon_{1}(n)$ of $\gamma_{h}^{1}(n)$ is sufficiently small such that $\gamma_{h}^{1}$, which lies on the invariant manifold $\kappa_{21}(S_{h},h)$, has to be in $N_{a,1}$ for $n<0$, and in $N_{r,1}$ for $n>0$ respectively, due to the uniqueness of the invariant center manifolds (see Proposition 3.3). In particular, observe that, if $h$ is small enough, $\gamma_{h}^{1}$ reaches an arbitrarily close vicinity of some $p_{a,1}(h_{1}^{*})$ for sufficiently large $n<0$ and of some $p_{r,1}(h_{1}^{*})$ for sufficiently large $n>0$, within $N_{a,1}\subset\widehat{M}_{a,1}$ and $N_{r,1}\subset\widehat{M}_{r,1}$ respectively (see also Figure 3 (b)). This finishes the proof. ∎

The trajectory $\gamma_{h}$ is shown in global blow-up coordinates as $\gamma_{\bar{h}}$ in Figure 3 (a), in comparison to the ODE trajectory $\bar{\gamma}_{0}$ corresponding to $\gamma_{0,2}$ in $K_{2}$.

3.7 Melnikov computation along the invariant curve

We consider a Melnikov-type computation for the distance between invariant manifolds, which is a discrete time analogue of continuous time results in [31] and, for a more general framework, in [43].

Consider an invertible map depending on a parameter $\mu$:

[TABLE]

where $(x,y)\in\mathbb{R}^{2}$, and $F=(F_{1},F_{2})^{\top}$ , $G=(G_{1},G_{2})^{\top}$ are $C^{k}$, vector-valued maps, $k\geq 1$. The following theory can be easily extended to $\mu\in\mathbb{R}^{m}$, like in [43], but for reasons of clarity we formulate it for $\mu\in\mathbb{R}$.

We formulate the following assumptions:

1. (A1)

There exist invariant center manifolds $M_{\pm}$ of the dynamical system (3.90), given as graphs of $C^{k}$-functions $y=g_{\pm}(x,\mu)$ and intersecting at $\mu=0$ along the smooth curve

[TABLE]

where $g_{\pm}(x,0)=g(x,0)$. 2. (A2)

Orbits of the map (3.90) with $\mu=0$ passing through a point $(x_{0},g(x_{0},0))$ on the invariant curve are given by a one-parameter family of solutions $(\gamma_{x_{0}}(n),0)^{\top}$ of dynamical system (3.90) with $\mu=0$, such that $\gamma_{x_{0}}(n)$ and $G(\gamma_{x_{0}}(n),0)$ are of a moderate growth when $n\to\pm\infty$ (to be specified later). 3. (A3)

There exist solutions $\phi_{\pm}(n)=(w_{\pm}(n),1)^{\top}$ of the linearization of (3.90) along $(\gamma_{x_{0}}(n),0)^{\top}$,

[TABLE]

such that

[TABLE]

and $w_{\pm}(n)$ are of a moderate growth (to be specified later) when $n\to\pm\infty$, respectively. 4. (A4)

The solutions $\psi_{x_{0}}(n)$ of the adjoint difference equation

[TABLE]

with initial vectors $\psi_{x_{0}}$ satisfying $\langle\psi_{x_{0}}(0),\partial_{x_{0}}\gamma_{x_{0}}(0)\rangle=0$, rapidly decay at $\pm\infty$ (the rate of decay to be specified later).

For a given $x_{0}$, we define $\psi_{x_{0}}(0)$ to be a unit vector in $\mathbb{R}^{2}$ orthogonal to $\partial_{x_{0}}\gamma_{x_{0}}(0)$, and set

[TABLE]

the intersections $M_{\pm}\cap\Sigma$ are then given by $(\Delta_{\pm}(\mu)\psi_{x_{0}}(0),\mu)$, where $\Delta_{\pm}$ are $C^{k}$-functions.

The following proposition is a discrete time analogue of [31, Proposition 3.1].

Proposition 3.8**.**

The first order separation between $M_{+}$ and $M_{-}$ at the section $\Sigma$ is given by

[TABLE]

Proof.

Equations (3.91) and (3.92) read:

[TABLE]

There follows:

[TABLE]

Choose initial data $w_{\pm}(0)=\frac{\mathrm{d}\Delta_{\pm}}{\mathrm{d}\mu}(0)\psi_{x_{0}}(0)$. Assuming that the growth of $w_{\pm}(n)$ and the decay of $\psi_{x_{0}}(n)$ at $n\to\pm\infty$, mentioned in (A3) and (A4), are such that

[TABLE]

we derive:

[TABLE]

and

[TABLE]

From this formula (3.93) follows immediately. ∎

We now apply Proposition 3.8 (or, better to say, its generalization for the case of two parameters $\mu=(r_{2},\lambda_{2})$) to the Kahan map (3.54) in the rescaling chart $K_{2}$. First of all, we have to justify Assumptions (A1)–(A4) for this case. Assumption (A1) follows from the fact that for $\mu=(r,\lambda)=0$, the center manifolds $\widehat{M}_{a,2}$ and $\widehat{M}_{r,2}$ intersect along the curve $S_{h}$ given in (3.83). Assumption (A2) follows from the explicit formula (3.84) for the solution $\gamma_{h,x_{0}}$, as well as from formulas (3.64) for the functions $\hat{J}$ and similar formulas for the functions $\hat{G}$. Assumption (A3) follows from the existence of the center manifolds away from $\mu=(r,\lambda)=0$, established in Proposition 3.3. Turning to the assumption (A4), we have the following results.

Proposition 3.9**.**

For problem (3.54), the adjoint linear system (3.92),

[TABLE]

has the decaying solution

[TABLE]

where

[TABLE]

and

[TABLE]

We have:

[TABLE]

Here the symbol $\approx$ relates quantities whose quotient has a limit as $n\to\pm\infty$.

Proof.

Fix $x_{0}\in\mathbb{R}$ and set

[TABLE]

Let

[TABLE]

be a fundamental matrix solution of the linear difference equation

[TABLE]

with $\det\Phi(0)=1$. The first column of the fundamental matrix solution $\Phi(n)$ can be found as $\partial_{x_{0}}\gamma_{h,x_{0}}$. Using formula (3.84) for $\gamma_{h,x_{0}}$, we have:

[TABLE]

A fundamental solution of the adjoint difference equation

[TABLE]

is given by

[TABLE]

Its second column is a solution of the adjoint system as given in (3.95), with $X(n)=\det\Phi(n)$. To compute $X(n)$, we observe that from

[TABLE]

and from $\det\Phi(0)=1$, there follows a discrete analogue of Liouville’s formula: for $n>0$,

[TABLE]

which coincides with (3.96) with $a(k)=\det A(k)$. Expression (3.97) for these quantities follows from (3.80).

To prove the estimate (3.98), we observe:

[TABLE]

Therefore, for $n>0$,

[TABLE]

Using the formula $\Gamma(n+c)\sim n^{c}\Gamma(n)$ by $n\to+\infty$ (in the sense that the quotient of the both expressions tends to 1), we obtain for $n\to+\infty$:

[TABLE]

This completes the proof. ∎

With the help of estimates of Proposition 3.9, we derive from Proposition 3.8 the following statement:

Proposition 3.10**.**

For the separation of the center manifolds $\widehat{M}_{a,2}$ and $\widehat{M}_{r,2}$, and for sufficiently small $h$, we have the first order expansion

[TABLE]

where $\mathcal{O}(2)$ denotes terms of order $\geq 2$ with respect to $\lambda,r$, and

[TABLE]

[TABLE]

In particular, convergence of the series in equation (3.101) is obtained for any $h>0$ and convergence of the series in equation (3.102) is obtained for $0<h<\sqrt{4/3}$.

Proof.

The form of the first order separation follows from Proposition 3.8. Furthermore, recall from equation (3.64) that

[TABLE]

Using Proposition 3.9, this yields (3.101) for any $h>0$. Note from equation (3.9) that the highest order $n^{\kappa}$ we can obtain in the terms $\hat{G}(\gamma_{h}(n),h)$ is $\kappa=3$ (coming from the term with factor $a_{2}$) such that for large $\left|n\right|$ we have

[TABLE]

This means that the convergence in (3.102) is given for $-4/h^{2}+2<-1$ such that the claim follows. ∎

We are now prepared to show our main result.

Theorem 3.11**.**

Consider the Kahan discretization for system (3.9). Then there exist $\varepsilon_{0},h_{0}>0$ and a smooth function $\lambda_{c}^{h}(\sqrt{\varepsilon})$ defined on $[0,\varepsilon_{0}]$ such that for $\varepsilon\in[0,\varepsilon_{0}]$ and $h\in(0,h_{0}]$ the following holds:

The attracting slow manifold $S_{a,\varepsilon,h}$ and the repelling slow manifold $S_{r,\varepsilon,h}$ intersect, i.e. exhibit a maximal canard, if and only if $\lambda=\lambda_{c}^{h}(\sqrt{\varepsilon})$. 2. 2.

The function $\lambda_{c}^{h}$ has the expansion

[TABLE]

where $C$ is given as in (2.30) (for $a_{3}=0$).

Proof.

First, we will work in chart $K_{2}$ and show that the quantities $d_{h,x_{0},\lambda}$, $d_{h,x_{0},r}$ in (3.101), (3.102) with $x_{0}=0$ approximate the quantities $d_{\lambda}$, $d_{r}$ in (A.21), (A.22)(up to change of sign). We prove:

[TABLE]

where, recall,

[TABLE]

[TABLE]

the function $\hat{J}$ is defined as in (3.64), and similar formulas hold true also for the function $\hat{G}$. Further recall that the Melnikov integrals can be solved explicitly, yielding

[TABLE]

where $a_{i}$ and $C$ are as introduced in Section 2.2 (for $a_{3}=0$, see (3.6) and (3.9)).

We show (3.103) — the simpler case (3.104) then follows similarly. We observe:

The remainder of the integral satisfies

[TABLE]

for $T>0$ and some $M\in\mathbb{N}$. Hence, we can keep $S(T)=\mathcal{O}(h^{2-c})$ for any $c>0$ with the choice $T\geq(4\ln\frac{1}{h})^{1/2}$. 2. 2.

For $N=T/h$, we turn to estimate

[TABLE]

We denote by $n^{*}$ the closest integer to $\alpha=2/h^{2}+1/2$, and recall that $\beta=2/h^{2}+3/2$. Since

[TABLE]

we can write, for all $n\geq 2/h^{2}+3/2$,

[TABLE]

Since, with Proposition 3.9 the summands of $\hat{S}(N)$ converge to zero even faster for smaller $h$, we obtain by choosing $N\geq\left\lceil{2/h^{2}+3/2}\right\rceil$, and hence $T\geq 2/h+5h/2$, that

[TABLE] 3. 3.

For $T=3/h$, we get by the standard methods the estimate

[TABLE]

Hence, we can conclude that equations (3.103) and (3.104) hold, and, in particular, that $d_{h,0,\lambda}$ and $d_{h,0,r}$ are bounded away from zero for sufficiently small $h$. Recall from (3.100) that

[TABLE]

where $D_{h,0}(0,0)=0$. Hence, the fact that $d_{h,0,\lambda}$ and $d_{h,r}$ are not zero implies, by the implicit function theorem, that there is a smooth function $\lambda^{h}(r)$ such that

[TABLE]

in a small neighborhood of $(0,0)$. Transforming back from $K_{2}$ into original coordinates then proves the first claim.

Furthermore, we obtain

[TABLE]

Transformation into original coordinates gives

[TABLE]

Hence, the second claim follows. ∎

Numerical computations show that $h_{0}$ in Theorem 3.11 does not have to be extremely small but that our results are quite robust for different step sizes. In Figure 4, we display such computations for the case $a_{1}=1$, $a_{2}=a_{4}=a_{5}=0$. In this case, the rescaled Kahan discretization in chart $K_{2}$ is given by

[TABLE]

Hence, we obtain

[TABLE]

For different values of $h$ and $N$ we calculate

[TABLE]

and, for the situation of (3.109) with $\hat{G}$ as in (3.110),

[TABLE]

We compare these quantities with the values of the respective continuous-time integrals $d_{\lambda}=-\sqrt{2\pi}$ and $d_{r}=-\sqrt{2\pi}/2$ (we have $C=1/2$ in this case).

We observe in Figure 4 that the sums converge very fast for relatively small $hN$ in both cases. Additionally, we see that $\left|d_{h,\lambda}(N)-d_{\lambda}\right|$ is significantly smaller than $\left|d_{h,r}(N)-d_{r}\right|$ for the same values of $h$. Note that the computations indicate that Theorem 3.11 holds for the chosen values of $h$ since $d_{h,0,\lambda}\approx\sqrt{2\pi}+\left(d_{\lambda}-d_{h,\lambda}(N)\right)$ is clearly distant from [math].

3.8 Numerical illustrations for $\varepsilon>0$

We illustrate the results by some additional numerics for $\varepsilon>0$, supplemenenting the illustrations of the dynamics in the rescaling chart $K_{2}$, as given by Figures 1 and 2. Firstly, we consider the simplest case where $a_{i}=0$ for all $i$, i.e., situation (3.12) with invariant curve $S_{\varepsilon,h}$ (3.15). Figure 5 shows different trajectories of the map (3.12) for $\varepsilon=0.1$ and $h=0.02$, illustrating the organization of dynamics around $S_{\varepsilon,h}$ analogously to the dynamics of (3.66) around $S_{h}$ (3.83) (see Figure 1).

Secondly, we consider the map (3.16) with $a_{1}=1$, i.e., a small additional perturbation of the canonical form, similarly to the end of the previous section. We take $\varepsilon=0.1$, $h=0.02$ and $\lambda=-(a_{1}/2)\varepsilon$, as a leading order approximation of $\lambda_{c}^{h}(\sqrt{\varepsilon})$ (see Theorem 3.11). In Figure 6, we observe that the numerics given by the Kahan discretization approximate very well the maximal canard, which slightly deviates from $S_{\varepsilon,h}$, again illustrating the organization of dynamics into bounded and unbounded trajectories sperated by the maximal canard. Note that we have chosen $\varepsilon=0.1$ to demonstrate the extension up to a relatively large $\varepsilon$.

In addition, we consider a model with cubic nonlinearity in order to demonstrate the application of the Kahan method beyond the purely quadratic case. Consider the equation

[TABLE]

as an example of equation (2.23), i.e., $a_{3}=1/3$ and $a_{i}=0,i=1,2,4,5$. Equation (3.113) is the van der Pol equation with constant forcing after transformation around one of the fold points (see [35, Example 8.1.6]). The Kahan discretization (3.10) of this equation yields

[TABLE]

such that the cubic nonlinearity does not vanish and we do not directly obtain an explicit form. However, we can use (3.114) as a numerical scheme by always taking the unique real solution $\tilde{x}$, closest to $x$ in absolut value, of the cubic polynomial.

In Figure 7, we illustrate the results of the Kahan discretization (3.114) of the van der Pol equation (3.113), again for $\varepsilon=0.1$ and $h=0.02$, taking $\lambda=-(3a_{3}/8)\varepsilon$, as a leading order approximation of $\lambda_{c}(\sqrt{\varepsilon})$ (see Theorem 2.2). Observe that the numerics indicate the existence of a maximal canard, also in this situation, separating bounded, now spiralling, orbits and unbounded orbits. Note that the implementation is based on the fact that the cubic polynomial in $\tilde{x}$ always has exactly one real solution, which we take as the next value, plus a complex conjugate pair with non-trivial imaginary part. A more general, algebraic analysis extends beyond the scope of this work and is left for additional research.

Note that the results on maximal canards for perturbations of the canonical form are local and do not make statements on the global stability. The preservation of canards for the van der Pol equation as depicted in Figure 7 is apparently also of predominantly local nature. Hence, we take a closer look in Figure 8, zooming into a neighbourhood of the inward-spiralling orbits from Figure 7. Here, we observe that the Kahan discretization even seems to capture the occurrence of a Hopf bifurcation in a neighbourhood of the maximal canard, as we slightly vary the parameter $\lambda$. Furthermore, the scheme seems to avoid crossing trajectories near the fold, which do occur as spurious solutions for some forward numerical methods near maximal canards. Indeed, there are also robust methods from boundary value problems (BVPs) [10, 23] and control theory [16, 27] to track canards for the van der Pol equation. However, these approaches do not take direct advantage of the polynomial structure, nor of the particular locally approximately integrable or symmetry structures of the van der Pol equation. Hence, building on the presented insights for the Kahan method, we consider an analytical treatment of the discretized cubic canard problem an intriguing direction for future work.

4 Conclusion

Our results show the importance of combining geometric invariants or integrable structures hidden in blow-up coordinates with suitable discretization schemes. Although we have just treated a very low-dimensional fast-slow fold case, one anticipates similar results also to be relevant for various other higher-dimensional singularities and bifurcation points, where blow-up is a standard tool. For example, it is well-known that in the Bogdanov-Takens unfolding one obtains small homoclinic orbits via a hidden integrable structure visible only after re-scaling. A thorough discretization analysis of higher-dimensional canards, similarly to the one at hand, would also deserve further investigation.

From a numerical perspective, forward integration schemes often provide an exploratory perspective to actually detect interesting dynamics or find a suitable invariant solution for fixed parameter values. In several cases these particular forward solutions are then used as starting conditions in numerical continuation techniques [11, 12] to study parametric dependence in a setting of BVPs. BVPs have also been successfully adapted to parametrically continue canard-type solutions [9, 10, 24, 33]. In particular, BVPs for canards turn out to be well-posed with a small numerical error, yet to set up the problem purely by continuation one already needs a very good understanding of phase space for the initial canard orbits. Therefore, a direct numerical integration scheme can be very helpful to automatically yield suitable starting solutions close to a maximal canard.

The Kahan method has mainly turned out to be favorable, since explicit, for quadratic vector fields; hence, in our analysis we have focused on this situation. However, we have seen in the numerical investigations in Section 3.8 that, by using its implicit form, also non-quadratic problems can be tackled, at least numerically. A further investigation into the dynamical and algebraic properties of the scheme, in particular for cubic nonlinearities, is a highly intriguing research question for the future, in general, and also in particular with respect to geometric multiscale problems as the one presented in this work.

Appendix A Existence of maximal canards for ODEs

In order to use specific geometric methods in singular perturbation theory, we consider $\varepsilon$ and $\lambda$ as variables, writing equation (2.20) as

[TABLE]

Note that in equation (2.23) the Jacobi matrix of the vector field in $(x,y,\lambda,\varepsilon)$ has a quadruple zero eigenvalue at the origin. A well established way to gain (partial) hyperbolicity at such a singularity is the blow-up technique which replaces the singularity by a manifold on which the dynamics can be desingularized. An important technical assumption for this technique is quasi-homogeneity of the vector field $f:\mathbb{R}^{n}\to\mathbb{R}^{n}$ of the ODE (cf. [35, Definition 7.3.2]), which means that there are $(a_{1},\dots,a_{n})\in\mathbb{N}^{n}$ and $k\in\mathbb{N}$ such that for every $r\in\mathbb{R}$ and each component $f_{j}:\mathbb{R}^{n}\to\mathbb{R}$ of $f$ we have

[TABLE]

The proof of Theorem 2.2 in [30] uses the quasi-homogeneous blow-up transformation $\Phi:B\to\mathbb{R}^{4}$,

[TABLE]

where $(\bar{x},\bar{y},\bar{\varepsilon},\bar{\lambda},r)\in B=S^{2}\times[-\kappa,\kappa]\times[0,\rho]$, where $S^{2}=\{(\bar{x},\bar{y},\bar{\varepsilon})\,:\,\bar{x}^{2}+\bar{y}^{2}+\bar{\varepsilon}^{2}=1\}$, with some $\kappa,\rho>0$. We assume that $\rho$ and $\kappa$ sufficiently small, so that the dynamics on $\Phi(B)$ can be described by the normal form approximation. Let $\overline{X}=\Phi^{*}(X)$ be the pull-back of the vector field $X$ to $B$. The dynamics of $\overline{X}$ on $B$ are analyzed in two charts $K_{1}$, $K_{2}$:

• •

the entering and exiting chart $K_{1}$ projecting the neighborhood of $(0,1,0)$ on $S^{2}$ to the plane $\bar{y}=1$:

[TABLE]

• •

and the scaling chart $K_{2}$ projecting the neighborhood of $(0,0,1)$ on $S^{2}$ to the plane $\bar{\varepsilon}=1$:

[TABLE]

The dynamics in the chart $K_{2}$ is of a primary interest. Here, the transformed equations admit a time rescaling allowing to divide out a factor $r_{2}$, which is possible due to the quasi-homogeneity of the leading part of the vector field $X$ (the new time being denoted by $t_{2}=r_{2}t$). Upon this operation, equations of motion take the form

[TABLE]

where $G=(G_{1},G_{2})$ can be written explicitly as

[TABLE]

On the invariant set $\{r_{2}=0,\lambda_{2}=0\}$, we have

[TABLE]

Let us list some crucially important qualitative features of system (A.14).

• •

As pointed out in [30, Lemma 3.3], system (A.14) possesses an integral of motion

[TABLE]

• •

Moreover, one can put (A.14) as a generalized Hamiltonian system

[TABLE]

• •

As a generalized Hamiltonian system, (A.14) preserves the measure $\textnormal{e}^{-2y_{2}}\mathrm{d}x_{2}\wedge\mathrm{d}y_{2}$. Since the density of an invariant measure is defined up to a multiplication by an integral of motion, the following is an alternative invariant measure:

[TABLE]

• •

System (A.14) has an equilibrium of center type at $(0,0)$, surrounded by a family of periodic orbits coinciding with the level curves $\{H(x_{2},y_{2})=c\}$ for $0<c<\frac{1}{2}$. The level curves for $c<0$ correspond to unbounded solutions. These two regions of the phase plane are separated by the invariant curve $\{H(x_{2},y_{2})=0\}$, or

[TABLE]

Thus, we have two alternative characterizations of the separatrix (A.18): on one hand, it is the level set $\{H(x_{2},y_{2})=0\}$, and on the other hand, it is the singular curve of the invariant measure (A.17).

• •

Separatrix (A.18) supports a special solution of (A.14):

[TABLE]

Pulled back to the manifold $B$, the special solution $\bar{\gamma}_{0}$ connects the endpoint $p_{a}$ of the critical attracting manifold $S_{a}$ across the sphere $S^{2}$ to the endpoint $p_{r}$ of the critical repelling manifold $S_{r}$ (see e.g. [35, Figure 8.2]). In other words, the center manifolds $\overline{M}_{a}$ and $\overline{M}_{r}$, corresponding to $p_{a}$ and $p_{r}$ respectively, and written in chart $K_{2}$ as $M_{a,2}$ and $M_{r,2}$, intersect along $\gamma_{0,2}$ for $r_{2}=\lambda_{2}=0$.

The difference between $M_{a,2}$ and $M_{r,2}$ for $(r_{2},\lambda_{2})\neq(0,0)$ is measured by the difference $y_{a,2}(0)-y_{r,2}(0)$, where $\gamma_{a,2}(t)=(x_{a,2}(t),y_{a,2}(t))$ and $\gamma_{r,2}(t)=(x_{r,2}(t),y_{r,2}(t))$ are the trajectories in $M_{a,2}$ and $M_{r,2}$ respectively, for given $r_{2},\lambda_{2}$ with the initial data $x_{a,2}(0)=x_{r,2}(0)=0$. This distance can be expressed as [30, Proposition 3.5]

[TABLE]

where

[TABLE]

are the respective Melnikov integrals. Since $d_{\lambda}\neq 0$, one concludes by the implicit function theorem that for sufficiently small $r_{2}$ there exists $\lambda_{2}$ such that the manifolds $M_{a,2}$ and $M_{r,2}$ intersect. Transforming back into the original variables yields Theorem 2.2.

It will be important for us that formulas (A.21), (A.22) admit also a non-Hamiltonian expression given in [43], where $\operatorname{grad}H(\gamma_{0,2}(t_{2}))$ is replaced by

[TABLE]

The function $\psi(t_{2})$ admits a more intrinsic interpretation as the only exponentially decaying solution of the adjoint system for the system (A.14) linearized along the solution $\gamma_{0,2}(t_{2})$,

[TABLE]

while the expression $2y_{0,2}(t_{2})=t_{2}^{2}/2$ in the exponent is interpreted as

[TABLE]

the matrix of the system (A.14) linearized along the solution $\gamma_{0,2}(t_{2})$ being given by

[TABLE]