Dynamics of Numerics of Nonautonomous Equations with Periodic Solutions: Introducing the Numerical Floquet Theory

Abstract

Nonautonomous systems with periodic solutions are encountered frequently in applications. In this paper, we will consider simple systems whose solutions are periodic with a known period. Their transformation under linearized collocation methods is investigated, using a technique called stroboscopic sampling, a discrete version of the well-known Poincaré map. It is shown that there is an inextricable relationship between AN stability (or BN stability) of the numerical methods and the correct qualitative behaviour of solutions.

1. Introduction

Let 𝑥^′=𝑓(𝑥,𝑡),𝑥(0)=𝑥₀, (1) where 𝑓∶𝐼⊂ℝ×ℝ→ℝ, be a scalar ordinary differential equation in which 𝑓(𝑥,𝑡) is a periodic function of 𝑡 with prime period 𝑇.

The detailed dynamics of numerics for nonautonomous ODEs has notably been lacking. Although any nonautonomous ODE can be transformed to an autonomous one, thereby increasing the dimension by one, the familiar dynamics of autonomous equations which is centered around the notion of equilibrium points [1] is lost. In certain special cases, this notion is replaced by that of periodicity. It is on these special cases that we will focus our attention. Stuart [2] proved using the bifurcation theory that for reaction-diffusion-convection equations, linearized instability implies the existence of spurious periodic solutions. Our approach differs from that of Stuart, who considered partial differential equations. Nonautonomous ODEs where 𝑓 is periodic in 𝑡 are very common in applications such as population dynamics with seasonal parameters or periodically forced systems.

Under certain conditions on 𝑓, (1) has a unique 𝑇-periodic solution [3]. We will assume that the solution is approximated by a linearized one-point collocation method as in Foster and Khumalo [4]. Our objective is to determine, for each 𝑓 under consideration, whether the numerical scheme has the same dynamical behaviour as the differential equation. In particular, we will consider cases in which the ODE has a unique, asymptotically stable periodic solution and establish conditions under which the numerical methods have the same dynamics. These special cases will take the following form:(i)𝑓linear: 𝑓(𝑥,𝑡)=𝑎(𝑡)𝑥+𝑏(𝑡), where 𝑎(𝑡) and 𝑏(𝑡) are 𝐶¹𝑇-periodic functions of 𝑡,(ii)𝑓 nonlinear: 𝑓=𝑎(𝑡)𝑔(𝑥)+𝑏(𝑡), where 𝑔(𝑥) is a 𝐶² nonlinear function of 𝑥.

The linearized one-point collocation methods for the scalar nonautonomous equation (1) are given by 𝑥_𝑛+1=𝑥_𝑛+ℎ𝑓𝑡_𝑛,𝑥_𝑛+𝑐₁ℎ²𝜕𝑓/𝜕𝑡𝑡_𝑛,𝑥_𝑛−𝑐₁ℎ𝑥_𝑛𝜕𝑓/𝜕𝑥𝑡_𝑛,𝑥_𝑛×1−𝑐₁ℎ𝜕𝑓/𝜕𝑥𝑡_𝑛,𝑥_𝑛⁻¹, (2) where 0≤𝑐₁≤1.

For a discussion of collocation methods in general, see Hairer et al. [5].

We begin with a description of the dynamical systems theory approach, which will be used in determining the conditions under which the methods have the same dynamical behaviour as the differential equations. Upon establishing these conditions, we compare them with those imposed by nonautonomous stability theory.

1.1. Dynamical Systems Approach

In what follows, we will use a technique known as stroboscopic sampling to reduce the problem of determining existence and stability of periodic solutions to existence and stability of fixed points.

Let 𝑥_𝑛+1=𝑝(𝑥_𝑛;𝑛;ℎ) be the discrete system representing the numerical method, applied with fixed stepsize, to the nonautonomous differential equation.

Step 1. Using inductive arguments, write the method in the form 𝑥_𝑛+1=𝜙(𝑥₀;𝑛;ℎ).

Step 2. Choose ℎ such that the period 𝑇=ℎ𝑘. Then, 𝑥_𝑘=𝜙(𝑥₀;𝑘;ℎ), and then we establish the new discrete system 𝑋_𝑛+1=𝜙(𝑋_𝑛;𝑘;ℎ). This is known as stroboscopic sampling.

Step 3. The fixed points of the last system correspond to 𝑇-periodic solutions of the method. These are determined with their stability types.

The above procedure is analogous to the Poincaré map of (1). Let Φ(𝑡,𝑥₀) be the 𝑇-periodic solution of (1), with the starting value 𝑥(0)=𝑥₀. Then, the Poincaré map of (1) is the scalar mapping Π∶ℝ⟶ℝ;𝑥⟼Φ(𝑇,𝑥). (3)

2. Linear Case

Suppose that 𝑓=𝑎(𝑡)𝑥+𝑏(𝑡), where 𝑎 and 𝑏 are 𝑇-periodic functions of 𝑡. Then, the linear nonautonomous differential equation (1) becomes 𝑥^′=𝑎(𝑡)𝑥+𝑏(𝑡),for𝑡≥0,𝑥(0)=𝑥₀. (4) If 𝜈₁=∫^𝑇₀𝑎(𝑡)𝑑𝑡<0, then (4) has a unique 𝑇-periodic asymptotically stable solution. If 𝑎(𝑡)=0, then (4) has a unique 𝑇-periodic solution that is asymptotically stable if 𝜈₂=∫^𝑇₀𝑏(𝑠)𝑑𝑠=0 (Hale and Koçak, [3]).

Now, the linearized one-point collocation methods are given by 𝑥_𝑛+1=𝑥_𝑛+ℎ1−ℎ𝑐₁𝑎()𝑎()𝑥𝑛ℎ𝑛ℎ_𝑛+𝑏(𝑛ℎ)+ℎ𝑐₁𝑎^′(𝑛ℎ)𝑥_𝑛+𝑏^′(𝑛ℎ). (5)

2.1. Linear Case with 𝑎(𝑡)=0

Theorem 1. Suppose that a linearized one-point collocation method is used to solve the linear nonautonomous differential equation (4) with 𝑎(𝑡)=0. The method tends (as 𝑛→∞, ℎ>0 fixed) to a periodic solution for any starting value if and only if 𝑏̃^𝑘−1_𝑟=0𝑟=0, (6) where 𝑏_𝑟=𝑏(𝑟ℎ)+ℎ𝑐₁𝑏^′(𝑟ℎ) for each 𝑟.

Proof. Assume that 𝜈₂=∫^𝑇₀𝑏(𝑠)𝑑𝑠=0. Taking 𝑎(𝑡)=0 in (5), we obtain 𝑥_𝑛+1=𝑥_𝑛+ℎ𝑏(𝑛ℎ)+ℎ𝑐₁𝑏^′(𝑛ℎ). (7) We define 𝑏̃_𝑖=𝑏(𝑖ℎ)+ℎ𝑐₁𝑏^′(𝑖ℎ)for𝑖=0,1,2,… (8) and denote Π𝑥_𝑛=𝑥_𝑛+ℎ𝑏(𝑛ℎ)+ℎ𝑐₁𝑏^′(𝑛ℎ), (9) so that for a given value of 𝑥₀, 𝑥_𝑛+1=𝑥₀+ℎ^𝑛_𝑟=0𝑏̃_𝑟. (10)
Fix 𝑘∈ℕ such that 𝑇=ℎ𝑘. Consider the 𝑘th iterate of 𝑥₀ under Π: 𝑥_𝑘=Π^𝑘𝑥₀=𝑥₀+𝑇𝑘^𝑘−1_𝑟=0𝑏̃_𝑟 (11) and the related iteration, which corresponds to stroboscopic sampling 𝑋_𝑛+1=𝑋_𝑛+𝑇𝑘^𝑘−1_𝑟=0𝑏̃_𝑟, (12) where 𝑋₀=𝑥₀. If the summation on the right-hand side of (12) is zero, then the discrete system is fixed at 𝑋₀ for all 𝑛, which corresponds to a periodic solution. If it is nonzero, then the stroboscopic iteration has no fixed point and diverges.

Remark 2. The second term on the right-hand side of (12) can be viewed as an application of the left rectangular quadrature rule to the integral ^𝑇₀𝑏(𝑠)+𝑐₁𝑇𝑘𝑏^′(𝑠)𝑑𝑠, (13) and (12) can be written as the simple map 𝑋_𝑛+1=𝑋_𝑛+𝑐, (14) where |𝑐|=𝑇𝑘^𝑘−1_𝑟=0𝑏̃_𝑟. (15)
The above theorem can then be restated as follows. Suppose that a linearized one-point collocation method is used to solve the linear nonautonomous differential equation (4) with 𝑎(𝑡)=0. The method tends (as 𝑛→∞, ℎ>0 fixed) to a periodic solution for any starting value if and only if the rectangular quadrature rule, used to approximate the integral in (13), gives a zero.

Illustration. We examine the stroboscopic sampling of the numerical solution of the differential equation 𝑥^′=cos𝑡𝑒^sin𝑡, 𝑥(0)=1. Figure 1 shows the numerical results for 𝑘=4, 5, and 10 with 𝑐₁=0.5. For 𝑘=4, the method diverges quicker from the periodic solution than for 𝑘=5. For 𝑘=10, the divergence is negligible.

Figure 1 

Numerical results for (12) with 𝑐1=1/2. 𝑘=4 (××), and 𝑘=5 (—), 𝑘=10 (------).

2.2. Linear Case with 𝑎(𝑡)=−1

If 𝑎(𝑡) is a negative constant, (4) could be scaled in such a way that 𝑎(𝑡)=−1. Then, clearly, 𝜈₁=∫^𝑇₀𝑎(𝑠)𝑑𝑠<0, and the differential equation has a unique, asymptotically stable periodic solution.

Theorem 3. Suppose that a linearized one-point collocation method is used to solve the linear nonautonomous differential equation (4) with 𝑎(𝑡)=−1. Then, for fixed 𝑐₁ and 𝑘, the method admits a unique periodic solution that is asymptotically stable, provided 𝑇𝑘>21−2𝑐₁. (16)

Proof. If 𝑎(𝑡)=−1, (5) simplifies to 𝑥_𝑛+1=𝑥_𝑛+𝑇𝑘+𝑇𝑐₁−𝑥_𝑛+𝑏𝑛𝑇𝑘+𝑇𝑐₁𝑘𝑏^′𝑛𝑇𝑘, (17) or, 𝑥_𝑛+1=𝑥_𝑛𝑘+𝑇𝑐₁−1𝑘+𝑇𝑐₁+𝑇𝑘+𝑇𝑐₁𝑏𝑛𝑇𝑘+𝑇𝑐₁𝑘𝑏^′𝑛𝑇𝑘. (18) Proceeding in a manner analogous to the above, we can show by induction that 𝑥_𝑛+1=𝑥₀𝑘+𝑇𝑐₁−1𝑘+𝑇𝑐₁^𝑛+1+𝑇𝑘+𝑇𝑐₁^𝑛_𝑟=0𝑘+𝑇𝑐₁−1𝑘+𝑇𝑐₁^𝑟𝑏̃_𝑛−𝑟, (19) where each 𝑏̃_𝑖 is defined by (8).
Denoting the right-hand side of (19) by Π(𝑥_𝑛), we can perform stroboscopic sampling and consider 𝑥_𝑘=Π^𝑘𝑥₀=𝑥₀𝑘+𝑇𝑐₁−1𝑘+𝑇𝑐₁^𝑘+𝑇𝑘+𝑇𝑐₁^𝑘−1_𝑟=0𝑘+𝑇𝑐₁−1𝑘+𝑇𝑐₁^𝑟𝑏̃_{𝑘−𝑟−1} (20) and the associated map 𝑋_𝑛+1=𝑋_𝑛𝑘+𝑇𝑐₁−1𝑘+𝑇𝑐₁^𝑘+𝑇𝑘+𝑇𝑐₁^𝑘−1_𝑟=0𝑘+𝑇𝑐₁−1𝑘+𝑇𝑐₁^𝑟𝑏̃_{𝑘−𝑟−1}. (21) Equation (21) is just the linear map 𝑋_𝑛+1=𝑐𝑋_𝑛+𝑑 (22) with 𝑐𝑐=𝑘+𝑇₁−1𝑘+𝑇𝑐₁^𝑘,𝑇𝑑=𝑘+𝑇𝑐₁^𝑘−1_𝑟=0𝑘+𝑇𝑐₁−1𝑘+𝑇𝑐₁^𝑟𝑏̃_{𝑘−𝑟−1}. (23) The map has a single fixed point, 𝑋^∗=𝑑,1−𝑐 (24) which is asymptotically stable if and only if |𝑐|<1; that is, 𝑇𝑘>1−2𝑐₁2, (25) and the result is established.

The following results are simple consequences of the above theorem.

Corollary 4. The linearized one-point collocation method with any 𝑐₁∈[1/2,1] applied to the linear nonautonomous differential equation with 𝑎(𝑡)=−1 admits a unique periodic solution that is asymptotically stable for all 𝑘>0.

Corollary 5. The explicit Euler method for the linear nonautonomous differential equation with 𝑎(𝑡)=−1 admits a unique, asymptotically stable periodic solution if and only if 𝑘>𝑇/2.

We now attempt to bound |𝑋^∗|. The following lemma gives a bound on the solution Φ(𝑡) of (4).

Lemma 6. Let 𝑏(𝑡) be a 𝐶¹𝑇-periodic function of 𝑡, and assume that (4) with 𝑎(𝑡)=−1 has a unique 𝑇-periodic solution. There exist a number 𝑀>0 such that |𝑏(𝑡)|, |𝑏^′(𝑡)|≤𝑀 and the 𝑇-periodic solution, Φ(𝑡), satisfies |Φ(𝑡)|≤𝑀 (26) for 𝑡→∞.

Proof. The boundedness of 𝑏(𝑡) and 𝑏^′(𝑡) follows from the periodicity and continuity of both functions.
Let 𝜙(𝑡) be a periodic solution of (4) with 𝑎(𝑡)=−1. Then, 𝜙(𝑡) satisfies the inequality −𝜙−𝑀≤𝜙^′≤−𝜙+𝑀∀𝑡. (27) This implies that 𝑒^−𝑡𝑥₀+𝑀−𝑀≤𝜙≤𝑒^−𝑡𝑥₀−𝑀+𝑀; (28) see [3]. Thus, 𝜙(𝑡) is bounded for 𝑡≥0; therefore, it approaches a 𝑇-periodic solution Φ(𝑡). Taking 𝑡→∞ in (28) establishes the lemma.

From the above lemma, we have |𝑥|≤𝑀. Then, we have the following theorem.

Theorem 7. Let 𝑋^∗, given by (24), denote the fixed point generated by the stroboscopic sampling of the numerical solution of the linear nonautonomous differential equation with 𝑎(𝑡)=−1 by a linearized one-point collocation method. Then, the inequality 𝑋^∗≤𝑀1+𝑇𝑐₁𝑘 (29) holds.

Proof. Since 𝑏(𝑡) is continuous, there is a number 𝑀 such that |𝑏(𝑡)|, |𝑏^′(𝑡)|≤𝑀 for all 𝑡∈ℝ. Then, for large 𝑡, the solution |Φ(𝑡)|<𝑀, and hence |𝑥|<𝑀.
Observe that 𝑏̃_𝑟≤𝑀₁=𝑀1+𝑐₁𝑇𝑘,|𝑑|≤𝑇𝑘+𝑇𝑐₁⋅𝑀₁⋅^𝑘−1_𝑟=0𝑘+𝑇𝑐₁−1𝑘+𝑇𝑐₁^𝑟=𝑇𝑘+𝑇𝑐₁⋅𝑀₁⋅𝑐1+𝑘+𝑇₁−1𝑇𝑐1−𝑘+𝑇₁−1𝑘+𝑇𝑐₁^𝑘−1=𝑀₁⋅𝑐1−𝑘+𝑇₁−1𝑘+𝑇𝑐₁^𝑘=𝑀₁(1−𝑐). (30)
Therefore, 𝑋^∗≤𝑀₁(1−𝑐)1−𝑐=𝑀₁=𝑀1+𝑇𝑐₁𝑘. (31)

Hence, 𝑋^∗ has essentially the same bound as the periodic solution.

Example 8. Consider the linear nonautonomous equation 𝑥^′=−𝑥+sin𝑡,𝑥(0)=0, (32) which has solution 𝑥(𝑡)=(1/2)[sin𝑡−cos𝑡+𝑒^−𝑡]. Figure 2 shows the numerical results (stroboscopic sampling) of the linearized implicit midpoint method (𝑐₁=1/2) with 𝑘=3 and 𝑘=10. Figure 3 shows the results of the explicit Euler method (𝑐₁=0) with 𝑘=3 and 𝑘=10. In these experiments, 𝑇=2𝜋; hence, the convergence to a unique periodic solution is expected for 𝑘>𝜋(1−2𝑐₁).

Figure 2 

Numerical results of (32) using linearized implicit midpoint method: 𝑘=3 (—) and 𝑘=10 (------).

Figure 3 

Numerical results of (32) using the explicit Euler method: 𝑘=3 (—) and 𝑘=10 (------).

2.3. Linear Case with 𝑎(𝑡)=−1+𝜖𝜌(𝑡)

Let 𝑎(𝑡)=−1+𝜖𝜌(𝑡), (33) where 𝜖≥0 is a constant and 𝜌(𝑡) is a 𝑇-periodic function of 𝑡. We assume that −𝑇+𝜖∫^𝑇₀𝜌(𝑠)𝑑𝑠<0 so that (4) has a unique asymptotically stable periodic solution.

2.3.1. Linearized One-Point Collocation Methods

We establish conditions under which a one-point collocation method with fixed 𝑐₁ and step-size ℎ>0 exhibits the same dynamical behaviour as the nonautonomous linear differential equation with 𝑎(𝑡)=−1+𝜖𝜌(𝑡).

Notation. In what follows, we will denote 𝑏_𝑛=𝑏(𝑛ℎ), 𝑏^′_𝑛=𝑏^′(𝑛ℎ), 𝜌_𝑛=𝜌(𝑛ℎ), and 𝜌^′_𝑛=𝜌(𝑛ℎ). Here, as before, 𝑇=ℎ𝑘 (𝑘∈ℕ).

Theorem 9. Suppose that a linearized one-point collocation method is used to solve a linear nonautonomous differential equation with 𝑎(𝑡)=−1+𝜖𝜌(𝑡). Then, for fixed 𝑐₁ and 𝑘, the method will admit a unique periodic solution that is asymptotically stable, provided ^𝑘−1_𝑖=01+𝑘𝑇−1+𝜖𝜌_𝑖+𝑐₁𝑇²𝜖𝜌^′_𝑖𝑘²−𝑘𝑇𝑐₁−1+𝜖𝜌_𝑖<1. (34)

Proof. From (5), we deduce that the linearized collocation methods, applied to the nonautonomous linear ODE with 𝑎(𝑡) given by (33), are 𝑥_𝑛+1=ℎ1+−1+𝜖𝜌_𝑛+𝑐₁ℎ²𝜖𝜌^′_𝑛1−ℎ𝑐₁−1+𝜖𝜌_𝑛𝑥_𝑛+ℎ1−ℎ𝑐₁−1+𝜖𝜌_𝑛⋅𝑏_𝑛+𝑐₁ℎ𝑏^′_𝑛. (35) We can write the above as 𝑥_𝑛+1=𝑇_𝑛𝑥_𝑛+𝐻_𝑛𝑏̃_𝑛∶=Π𝑥_𝑛, (36) where 𝑇_𝑛ℎ=1+−1+𝜖𝜌_𝑛+𝑐₁ℎ²𝜖𝜌^′_𝑛1−ℎ𝑐₁−1+𝜖𝜌_𝑛,𝐻_𝑛=ℎ1−ℎ𝑐₁−1+𝜖𝜌_𝑛,𝑏̃_𝑛=𝑏_𝑛+𝑐₁ℎ𝑏^′_𝑛. (37)
Proceeding by induction, we establish that 𝑥_𝑛+1=𝑥₀^𝑛_𝑖=0𝑇_𝑖+^𝑛_𝑖=0𝐻_𝑖𝑏̃_𝑖^𝑛_𝑗=𝑖+1𝑇_𝑗, (38) from which we deduce that 𝑥_𝑘=𝑥₀^𝑘−1_𝑖=0𝑇_𝑖+^𝑘−1_𝑖=0𝐻_𝑖𝑏̃_𝑖^𝑘−1_𝑗=𝑖+1𝑇_𝑗∶=Π^𝑘𝑥₀. (39)
The discrete system that corresponds to stroboscopic sampling is the linear system 𝑋_𝑛+1=𝑐(𝜖)𝑋_𝑛+𝑑(𝜖), (40) where 𝑐(𝜖)=^𝑘−1_𝑖=0𝑇_𝑖,𝑑(𝜖)=^𝑘−1_𝑖=0𝐻_𝑖𝑏̃_𝑖^𝑘−1_𝑗=𝑖+1𝑇_𝑗. (41) This system has a unique fixed point, 𝑋^∗, given by (24). It is asymptotically stable if and only if |𝑐(𝜖)|<1; that is, |∏^𝑘−1_𝑖=0𝑇_𝑖|<1, or ^𝑘−1_𝑖=0ℎ1+−1+𝜖𝜌_𝑖+𝑐₁ℎ²𝜖𝜌^′_𝑖1−ℎ𝑐₁−1+𝜖𝜌_𝑖<1. (42)
Substituting ℎ=𝑇/𝑘 gives the result.

2.3.2. Examples

For each of the three special values of 𝑐₁, 𝜌(𝑡)=sin𝑡, and the increasing values of 𝜖, we determined, using (42), the minimum value of 𝑘 such that each method has dynamical behaviour that is the same as that of the differential equation. The results are illustrated in Figure 4.

Figure 4 

Least 𝑘 for unique periodic solution. The explicit Euler (------), linearized implicit midpt (—), and the linearized implicit Euler (××).

For 𝜖<2, the explicit Euler method is the most restrictive of the three (i.e., comparatively larger minimum values of 𝑘 must be used to obtain dynamical behaviour that is the same as that of the differential equation). However, as 𝜖 is increased, the Explicit Euler method outperforms the linearized implicit Euler method by becoming less restrictive than that method for 𝜖≥3. For 𝜖≥3, the explicit Euler and linearized midpoint methods give comparable results, and for those values of 𝜖, the linearized implicit Euler method becomes more and more restrictive in comparison to the other two.

Finally, we develop a bound on |𝑋^∗|. For the comparison purposes, we present the following lemma, which can be proved in a manner analogous to Lemma 6.

Lemma 10. Let 𝑏(𝑡), 𝜌(𝑡) be 𝐶¹ 𝑇-periodic functions of 𝑡. There exist numbers 𝑀, 𝑀₂>0 such that |𝑏(𝑡)|, |𝑏^′(𝑡)|≤𝑀, |𝜌(𝑡)|, |𝜌^′(𝑡)|≤𝑀₂, and the solution Φ(𝑡) satisfies |Φ(𝑡≤𝑀)|1−𝜖𝑀₂ (43) as 𝑡→∞.

From the lemma, we have |𝑥|≤𝑀/|1−𝜖𝑀₂|.

Theorem 11. Let 𝑋^∗, given by (24), be the fixed point generated by the stroboscopic sampling of the numerical solution of the linear nonautonomous differential equation with 𝑎(𝑡)=−1+𝜖𝜌(𝑡) by a linearized collocation method. Then, the following inequality holds: 𝑋^∗≤𝑀1+𝑇𝑐₁/𝑘1−ℎ𝑐₁𝜖𝑀₂−𝜖𝑀₂. (44)

Proof. Since 𝑏(𝑡) and 𝜌(𝑡) are 𝐶¹ and periodic, there are numbers 𝑀, 𝑀₂ such that |𝑏(𝑡)|, |𝑏^′(𝑡)|≤𝑀 and |𝜌(𝑡)|, |𝜌^′(𝑡)|≤𝑀₂ for all 𝑡∈ℝ.
For each 𝑖, 𝑇_𝑖≤1+ℎ𝑐₁+ℎ𝜖𝑀₂−ℎ1+ℎ𝑐₁−ℎ𝑐₁𝜖𝑀₂,𝐻_𝑖≤ℎ1+ℎ𝑐₁−ℎ𝑐₁𝜖𝑀₂,𝑏̃_𝑖≤𝑀₁=𝑀1+𝑐₁ℎ, (45) where ℎ=𝑇/𝑘. If 𝑣(𝜖)=1+ℎ𝑐₁−ℎ𝑐₁𝜖𝑀₂+ℎ𝜖𝑀₂+ℎ²𝑐₁𝜖𝑀₂−ℎ1+ℎ𝑐₁−ℎ𝑐₁𝜖𝑀₂, (46) then |𝑑(𝜖≤ℎ)|1+ℎ𝑐₁−ℎ𝑐₁𝜖𝑀₂⋅𝑀₁⋅^𝑘−1_𝑖=0𝑣^{𝑘−𝑖−1}=ℎ𝑀₁1+ℎ𝑐₁−ℎ𝑐₁𝜖𝑀₂⋅1−𝑣^𝑘ℎ1−𝑐₁ℎ𝜖𝑀₂−𝜖𝑀₂⋅1+ℎ𝑐₁−ℎ𝑐₁𝜖𝑀₂=𝑀₁⋅1−𝑣^𝑘1−𝑐₁ℎ𝜖𝑀₂−𝜖𝑀₂. (47)
On the other hand, we deduce from the definition of 𝑐(𝜖) that |(𝜖≥1−𝑐)|1−1+ℎ𝑐₁−ℎ+ℎ𝜖𝑀₂−ℎ𝑐₁𝜖𝑀₂+ℎ²𝑐₁𝜖𝑀₂1+ℎ𝑐₁−ℎ𝑐₁𝜖𝑀₂^𝑘. (48) Hence, 𝑋^∗(𝜖≤𝑀)¹1−𝑐₁ℎ𝜖𝑀₂−𝜖𝑀₂=𝑀1+𝑐₁ℎ1−𝑐₁ℎ𝜖𝑀₂−𝜖𝑀₂, (49) which is identical to the inequality (44).

If we substitute 𝜖=0 in (49), we obtain (31) as expected. Here, as well, the bound for 𝑋^∗(𝜖) is the same as that of the periodic solution as ℎ→0.

We would like to obtain a relationship between the dynamical approach study and stability analysis. We introduce a natural stability criterion for the differential equation as well as any numerical method used to discretize it.

2.4. Conditional AN Stability and AN Stability

We consider the problem of determining a criterion for some sort of “controlled behaviour” of the solutions of the methods. We adopt a linear stability criterion that is based on the scalar test equation 𝑥^′=𝑎(𝑡)𝑥, (50) where 𝑎(𝑡)∈ℂ. If Re(𝑎(𝑡))<0 for all 𝑡∈[𝛽₁,𝛽₂], then 𝑥(𝑡+ℎ)=𝐾𝑥(𝑡),|𝐾|≤1 (51) for all 𝑥∈[𝛽₁,𝛽₂] and ℎ>0.

Definition 12. A numerical method is said to be conditionally AN stable for some ℎ>0 if, when applied to the test equation (50) with Re(𝑎(𝑡))<0, for all 𝑡, 𝑥_𝑛+1=𝐾(ℎ)𝑥_𝑛,𝐾(ℎ)≤1 (52) holds for all 𝑛∈ℕ.
If this condition is satisfied for all ℎ>0, then the method is AN stable.

For a discussion of AN stability, see Lambert [6] and Stuart and Humphries [7]

The following simple result gives a condition under which the linearized one-point collocation methods are conditionally AN stable.

Theorem 13. Assuming real 𝑎(𝑡), the linearized one-point collocation methods are conditionally AN stable if and only if 1−2𝑐₁𝑎(𝑡)+𝑐₁ℎ𝑎^′(𝑡)2≥−ℎ (53) and 𝑎(𝑡)+𝑐₁ℎ𝑎^′(𝑡)≤0.

Proof. Use the above test equation in (2). On the other hand, if (2) is not satisfied, the method fails the stability criterion and is not AN stable.
Examples:(1)if 𝑎(𝑡)=−1, then the linearized one-point collocation methods are conditionally AN stable if and only if ℎ21−2𝑐₁≤1. (54) It is easy to observe that the methods are AN stable if 𝑐₁≥1/2,(2)if 𝑎(𝑡)=−1+𝜖𝜌(𝑡), where 𝜖, 𝜌(𝑡)∈𝐂, and 𝜖(𝜌(𝑡)+𝑐₁ℎ𝜌^′(𝑡))≤1, then the linearized one-point collocation methods are conditionally AN stable if and only if −1+𝜖𝜌(𝑡)+𝑐₁ℎ𝜖𝜌^′(𝑡)(𝑡2+2−2𝜖𝜌)≥−ℎ. (55)
We have proved the existence of a relationship between the linear stability theory of the collocation methods and the existence and asymptotic stability of periodic solutions, identified via stroboscopic sampling. This relationship is stated in the following theorem.

Theorem 14. Suppose that a linearized one-point collocation method is used to solve a linear nonautonomous equation of the form discussed in Sections 2.2 or 2.3 which has periodic coefficients and possesses a unique, asymptotically stable periodic solution. Then, the following conditions are equivalent:(a)the method is AN stable,(b)the method yields the same dynamics as the differential equation.

The last theorem is very significant, since it gives us a bridge connecting standard stability theory with dynamical systems. Naturally, we would like to find out if there is a corresponding result for the nonlinear case, which we now consider.

3. Nonlinear Case

We consider the nonlinear equation, 𝑥^′=𝑎(𝑡)𝑔(𝑥)+𝑏(𝑡),for𝑡≥0,𝑥(0)=𝑥₀, (56) where 𝑏(𝑡) is a 𝑇-periodic function of 𝑡 and 𝑔(𝑥) is a 𝐶² nonlinear function of 𝑥.

Massera [8] proved that if a nonlinear equation of the form (56) has the uniqueness property with respect to the initial conditions, the existence of a bounded solution implies the existence of a 𝑇-periodic solution.

3.1. Linearized One-Point Collocation Methods

The linearized collocation methods, applied to (56), are given by 𝑥_𝑛+1=𝑥_𝑛+ℎ𝑎_𝑛𝑔𝑥_𝑛1−ℎ𝑎_𝑛𝑐₁𝑔^′𝑥_𝑛ℎ𝑏̃ℎ+^𝑛1−ℎ𝑎_𝑛𝑐₁𝑔^′𝑥_𝑛, (57) where 𝑏̃_𝑛=𝑏_𝑛+𝑐₁ℎ𝑏^′_𝑛 and 𝑎_𝑛=𝑎_𝑛+𝑐₁ℎ𝑎^′_𝑛.

We perform a simplification on the third term of (57) that takes the form of evaluating the derivative of 𝑔 at the starting value, instead, at each step. The resulting method, that will be referred to as a simplified linearized one-point collocation method, is 𝑥_𝑛+1=𝑥_𝑛+ℎ𝑎_𝑛𝐺𝑥_𝑛;𝑛+𝐻_𝑛𝑏̃_𝑛, (58) where 𝐻_𝑛=ℎ1−ℎ𝑎_𝑛𝑐₁𝑔^′𝑥₀,𝐺()=𝑔(𝑥)𝑥;𝑛1−ℎ𝑎_𝑛𝑐₁𝑔^′(𝑥). (59)

3.2. The Dynamical Systems Approach

We would like to take the dynamical systems approach and determine the conditions under which (58), applied to the differential equation (56), yields the same dynamics as the continuous system.

We rewrite (58) as ℎ𝑎_𝑛𝐺𝑥_𝑛;𝑛=𝑥_𝑛+1−𝑥_𝑛−𝐻_𝑛𝑏̃_𝑛. (60) Inductively, we can show that 𝑥_𝑛+1=𝑥₀+^𝑛_𝑟=0𝐻_𝑟𝑏̃_𝑟+ℎ^𝑛_𝑟=0𝑎_𝑟𝐺𝑥_𝑟;𝑟. (61)

We choose an integer 𝑘 such that 𝑇=ℎ𝑘. Sampling stroboscopically in the iteration above, we get 𝑥_𝑘=𝑥₀+^𝑘−1_𝑟=0𝐻_𝑟𝑏̃_𝑟+ℎ^𝑘−1_𝑟=0𝑎_𝑟𝐺𝑥_𝑟;𝑟 (62) and associate this with the discrete system 𝑋_𝑛+1=𝑋_𝑛+^𝑘−1_𝑟=0𝐻_𝑟𝑏̃_𝑟+ℎ^𝑘−1_𝑟=0𝑎_𝑟𝐺𝑋_𝑟;𝑟, (63) where 𝑋₀=𝑋_𝑛,𝑋_𝑟=𝑋_𝑛+^𝑟−1_𝑖=0𝐻_𝑟𝑏̃_𝑟+ℎ^𝑟−1_𝑖=0𝑎_𝑖𝐺𝑋_𝑖;𝑖for𝑟=1,2,…,𝑘−1. (64)

The fixed points of (63) correspond to periodic solutions of (56). Fixed points are points, 𝑋^∗, such that ^𝑘−1_𝑟=0𝐻_𝑟𝑏̃_𝑟+ℎ^𝑘−1_𝑟=0𝑎_𝑟𝐺𝑋_𝑟;𝑟=0, (65) where 𝑋_𝑟=𝑋^∗+^𝑟−1_𝑖=0𝐻_𝑖𝑏̃_𝑖+ℎ^𝑟−1_𝑖=1𝑎_𝑖𝐺𝑋_𝑖;𝑖for𝑟=1,2,…,𝑘−1, (66) and 𝑋₀=𝑋^∗. Define the sequence of functions 𝐹₁(𝑥), 𝐹₂(𝑥), … by 𝐹_𝑘(𝑥)=^𝑘−1_𝑟=0𝑎_𝑟ℎ𝐺𝑥_𝑟;𝑟+𝐻_𝑟𝑏̃_𝑟, (67) where 𝑥₀=𝑥 and 𝑥_𝑟=𝑥+^𝑟−1_𝑖=0𝑎_𝑖ℎ𝐺𝑥_𝑖;𝑖+𝐻_𝑖𝑏̃_𝑖, (68) for 𝑟=1,2,…,𝑘−1.

The theorem below gives conditions under which the simplified linearized one-point collocation methods, applied to (56), exhibit dynamical behaviour that is the same as the differential equation.

Theorem 15. Assume that the differential equation (56) has a unique solution. If(i)𝑔^′′(𝑥)𝑔(𝑥)≤0,(ii)𝑐₁𝑎^′(𝑡)<−𝑎(𝑡)/ℎ for all 𝑡,(iii)Max𝐺^′(𝑥;𝑡)≤1/ℎ|𝑎(𝑡)|,(iv)|𝐹^′_𝑘(𝑥)|<1 for all 𝑘∈ℕ, where 𝐹_𝑘(𝑥) is given by (67),
then a simplified one-point collocation method has a periodic solution for any 𝑘=𝑇/ℎ; this solution is unique and asymptotically stable.

Proof. Finding possible fixed points of (63), hence periodic solutions of the methods, is the same as finding zeros of (67). This is equivalent to solving the nonlinear system: −ℎ𝑎₀𝐺𝑥^∗;0−ℎ𝑎₁𝐺𝑥₁;1−ℎ𝑎₂𝐺𝑥₂𝑎;2−⋯−ℎ_𝑘−1𝐺𝑥_𝑘−1;𝑘−1−𝐻₀𝑏̃₀−𝐻₁𝑏̃₁−⋯−𝐻_𝑘−1𝑏̃_𝑘−1𝑥=0,₁−𝑥^∗−ℎ𝑎₀𝐺𝑥^∗;0−𝐻₀𝑏̃₀𝑥=0,₂−𝑥^∗−ℎ𝑎₀𝐺𝑥^∗;0−ℎ𝑎₁𝐺𝑥₁;1−𝐻₀𝑏̃₀−𝐻₁𝑏̃₁⋮𝑥=0,_𝑘−1−𝑥^∗−ℎ𝑎₀𝐺𝑥^∗;0−ℎ𝑎₁𝐺𝑥₁;1−ℎ𝑎₂𝐺𝑥₂𝑎;2−⋯−ℎ_𝑘−2𝐺𝑥_𝑘−2;𝑘−2−𝐻₀𝑏̃₀−𝐻₁𝑏̃₁−⋯−𝐻_𝑘−2𝑏̃_𝑘−2=0. (69) The above system is of the form 𝐅(𝐱)=𝟎, where 𝐅 is a nonlinear function of 𝑥𝐱=^∗=𝑥₀𝑥₁𝑥₂⋮𝑥_𝑘−1. (70)
To prove existence, it is sufficient to show that the Jacobian matrix, J(𝐅), of 𝐅 is nonsingular. Now, J(𝐅)=𝑑₀𝑑₁𝑑₂𝑑₃⋯𝑑_𝑘−2𝑑_𝑘−1−1+𝑑₀100⋯00−1+𝑑₀𝑑₁100⋯0−1+𝑑₀𝑑₁𝑑₂⋮10⋯0−1+𝑑₀𝑑₁𝑑₂⋯𝑑_𝑘−2𝑑_𝑘−11, (71) where 𝑑_𝑖=−ℎ𝑎_𝑖𝐺^′(𝑥_𝑖;𝑖) for each 𝑖. We perform one elementary row operation: row 1 → row 1–row 𝑘. The matrix becomes 1000⋯0−1+𝑑_𝑘−1−1+𝑑₀100⋯00−1+𝑑₀𝑑₁100⋯0−1+𝑑₀𝑑₁𝑑₂⋮10⋯0−1+𝑑₀𝑑₁𝑑₂⋯𝑑_𝑘−2𝑑_𝑘−11. (72) It is easy to see that the above matrix is nonsingular.
To prove uniqueness, let 𝑥₁ and 𝑥₂ be two fixed points of (63). Then, from (67), 𝐹_𝑘𝑥₁=𝐹_𝑘−1𝑥₁+ℎ𝑎_𝑘−1𝐺𝑥₁+𝐹_𝑘−1𝑥₁;𝑘−1+𝐻_𝑘−1𝑏̃_𝑘−1𝐹=0,_𝑘𝑥₂=𝐹_𝑘−1𝑥₂+ℎ𝑎_𝑘−1𝐺𝑥₂+𝐹_𝑘−1𝑥₂;𝑘−1+𝐻_𝑘−1𝑏̃_𝑘−1=0. (73)(74) Subtracting (74) from (73) and using the mean value theorem gives 𝑥₁−𝑥₂𝐹^′_𝑘−1𝛼₁+ℎ𝑎_𝑘−1𝐺^′𝛼₂;𝑘−11+𝐹^′_𝑘−1𝛼₁=0, (75) which may alternatively be written as 𝑥₁−𝑥₂𝐹^′_𝑘−1𝛼₁1+ℎ𝑎_𝑘−1𝐺^′𝛼₂𝑎;𝑘−1+ℎ_𝑘−1𝐺^′𝛼₂;𝑘−1=0. (76)
In the above equations, 𝛼₁ is between 𝑥₁ and 𝑥₂, and 𝛼₂ is between 𝑥₁+𝐹_𝑘−1(𝑥₁) and 𝑥₂+𝐹_𝑘−1(𝑥₂).
From the hypotheses of the theorem, we have that 𝐹^′<0, 𝐺^′>0, 𝑎_𝑘−1<0, and 1≥1+ℎ𝑎_𝑘−1𝐺^′(𝛼₂;𝑘−1)≥0. Thus, 𝑥₁=𝑥₂.
The periodic solution is asymptotically stable since |1+𝐹^′|<1. This completes the proof.

3.3. Numerical Experiments

Consider the nonlinear nonautonomous ODE 𝑥^′(𝑡)=−𝑥³+sin𝑡,𝑥(0)=1. (77)

The numerical results (stroboscopic sampling) of each of the three cases are given in Figures 5, 6, and 7.

Figure 5 

Numerical results for (77) with 𝑐1=0. 𝑘=20.

Figure 6 

Numerical results for (77) with 𝑐1=1/2. 𝑘=20.

Figure 7 

Numerical results for (77) with 𝑐1=1. 𝑘=20.

The methods give comparable results, but the implicit midpoint converges faster to the periodic solutions than the other two.

3.4. Nonlinear Stability Theory

We wish to establish conditions under which numerical methods for the solution of (56) behave in a “controlled” manner. We will view such controlled behaviour in the system meaning that neighbouring solution curves get closer and closer together as 𝑡 increases. This concept, called contractivity, is discussed, for instance, in Lambert [6]. As in the linear case, we will contrast the conditions for stability with the existence and uniqueness of a periodic solution in the numerical methods.

We briefly state the concepts of contractivity and conditional BN stability.

Definition 16. Let 𝑥(𝑡) and 𝑥(𝑡) be any two solutions of the differential equation 𝑥^′=𝑓(𝑥,𝑡), satisfying initial conditions 𝑥(0)=𝜂, 𝑥(0)=𝜂, 𝜂≠𝜂. Then, if 𝑥𝑡₂−𝑥𝑡₂≤𝑥𝑡₁−𝑥𝑡₁ (78) holds under the ℝ^𝑚 norm ‖⋅‖ for all 𝑡₁, 𝑡₂ such that 𝛽₁≤𝑡₁≤𝑡₂≤𝛽₂, the solutions of the system are said to be contractive in [𝛽₁,𝛽₂].

The discrete analog of the above definition is given below.

Definition 17. Let {𝑥_𝑛} and {𝑥_𝑛} be two numerical solutions generated by a numerical method with different starting values. Then, if 𝑥_𝑛+1−𝑥_𝑛+1≤𝑥_𝑛−𝑥_𝑛,0≤𝑛≤𝑁, (79) the numerical solutions are said to be contractive for 𝑛∈[0,𝑁].

Definition 18. The system 𝑥^′=𝑓(𝑥,𝑡) is dissipative in [𝛽₁,𝛽₂] if 𝑓(𝑥,𝑡)−𝑓(𝑥,𝑡),𝑥−𝑥≤0 (80) holds for all 𝑥, 𝑥∈ℝ and for all 𝑡∈[𝛽₁,𝛽₂].

It is easy to show that the solutions of a dissipative system are contractive under the norm induced by the inner product in (80). It is desirable that a numerical method, used with fixed stepsize ℎ>0 to solve a dissipative system, gives contractive solutions. This brings us to the concepts of conditional BN stability, BN stability and nonlinear nonautonomous stability criteria.

Definition 19. If a numerical method, applied with fixed steplength ℎ>0 to (56) satisfying (80), generates contractive solutions, the method is said to be conditionally BN stable. If the method generates contractive solutions when applied with any ℎ>0, then it it is BN stable (Lambert [6]).

The concepts of AN and BN stability are equivalent for the nonconfluent Runge-Kutta methods.

To determine the conditional BN stability of the methods, we use the scalar test system 𝑥^′=𝑎(𝑡)𝑔(𝑥), (81) where, as before, 𝑎(𝑡)∈ℂ and 𝑔(𝑥)∈𝐶¹. This system is dissipative if 𝑎(𝑡)𝑔(𝑥)−𝑔(𝑥),𝑥−𝑥=𝑎(𝑡)𝑔^′(𝜉)|𝑥−𝑥|²≤0, (82) where 𝜉 lies between 𝑥 and 𝑥 and ⟨⋅,⋅⟩ is an inner product in ℝ. The condition is satisfied if 𝑎(𝑡)≤0 for all 𝑡 and 𝑔^′(𝑥)≥0 for all 𝑥. Therefore, the existence and uniqueness of a periodic solution in (56) is a sufficient condition for the dissipativity of the system.

Now, we determine the conditions under which the linearized one-point collocation methods are conditionally BN stable. Applying the methods to the test system (81) for two different initial conditions gives 𝑥_𝑛+1=𝑥_𝑛+ℎ𝑎_𝑛𝐺𝑥_𝑛;𝑛,𝑥_𝑛+1=𝑥_𝑛+ℎ𝑎_𝑛𝐺𝑥_𝑛;𝑛. (83)(84)

Subtracting (84) from (83) and using the mean value theorem gives 𝑥_𝑛+1−𝑥_𝑛+1=𝑥_𝑛−𝑥_𝑛⋅1+ℎ𝑎_𝑛𝐺^′𝜉_𝑛;𝑛, (85) where 𝜉_𝑛 lies between 𝑥_𝑛 and 𝑥_𝑛.

The solutions generated by the methods are contractive if and only if |1+ℎ𝑎_𝑛𝐺^′(𝜉_𝑛;𝑛)|≤1. Assuming 𝑎(𝑡)<0 for all 𝑡 and 𝐺^′(𝑥;𝑡)≥0 for all 𝑥 and 𝑡, the methods are conditionally BN stable if and only if 𝐺^′()≤2𝑥;𝑡ℎ𝑎(𝑡). (86)

Remark 20. If we let 𝑔(𝑥)=𝑥 (the linear case), then condition (86) collapses to (53), which is the condition for conditional AN stability.

3.5. Discussion

We have, in Theorem 15, established sufficient conditions for the simplified linearized one-point collocation methods to exhibit the same dynamics as the nonlinear nonautonomous ODE. Conditions (i) and (ii) are required for conditional BN stability as well, but condition (iii) is not necessary for the existence of a unique, asymptotically stable solution. In fact, we can come to the same conclusion as in Theorem 15 if, in (75), 𝐹^′(𝑥)1+ℎ𝑎(𝑡)𝐺^′(𝑥;𝑡)+ℎ𝑎(𝑡)𝐺^′(𝑥;𝑡)≠0 (87) for all 𝑥 and 𝑡. The following theorem shows that condition (87) is satisfied if the method is conditionally BN stable.

Theorem 21. Suppose that a simplified linearized one-point collocation method is used to solve a nonlinear nonautonomous equation of the form (56) with periodic coefficients which has a unique, asymptotically stable periodic solution. Then, the method yields a unique periodic solution that is asymptotically stable if it is conditionally BN stable.

Proof. Assume the conditional BN stability. To prove the theorem, it is sufficient to establish condition (87). From conditional BN stability, we have −1≤1+ℎ𝑎𝐺^′≤1. Recall that −1≤𝐹^′<0 from condition (iv). Therefore, 𝐹^′+ℎ𝑎𝐺^′≤𝐹^′1+ℎ𝑎𝐺^′+ℎ𝑎𝐺^′≤−𝐹^′+ℎ𝑎𝐺^′. (88) Since 0<−𝐹^′≤1, if −1≤1+ℎ𝑎𝐺^′<0, 𝐹^′1+ℎ𝑎𝐺^′+ℎ𝑎𝐺^′≤1+ℎ𝑎𝐺^′<0. (89) If 0≤1+ℎ𝑎𝐺^′≤1, 𝐹^′+ℎ𝑎𝐺^′≤𝐹^′1+ℎ𝑎𝐺^′+ℎ𝑎𝐺^′≤ℎ𝑎𝐺^′<0. (90) This proves the theorem.

Remark 22. The above theorem is somewhat similar to Theorem 14 (in the linear case) if we replace the concept of conditional BN stability by conditional AN stability.

4. Conclusion

We already knew that numerical methods can introduce spurious behaviours into the solution for autonomous equations. By concentrating on linear and nonlinear nonautonomous equations with unique periodic solutions and discretizing them using one-point collocation methods, we were interested in the existence of periodic solutions in the numerical methods.

We found that the results obtained from the dynamical systems approach are closely linked to those that are imposed by standard stability analysis. It has been shown that for linear and nonlinear nonautonomous differential equations of the form considered in this chapter, there is a relationship between conditional AN or BN stability of a one-point collocation method and the method yielding the same dynamical behaviour as the differential equation under consideration.