On exact discretization of the cubic-quintic Duffing oscillator

# On exact discretization of the cubic-quintic Duffing oscillator

## Abstract

Application of intersection theory to construction of -point finite-difference equations associated with classical integrable systems is discussed. As an example, we present a few exact discretizations of one-dimensional cubic and quintic Duffing oscillators sharing form of Hamiltonian and canonical Poisson bracket up to the integer scaling factor.

## 1 Introduction

A completely integrable system on symplectic manifold with form of dimension is defined by smooth functions in involution

 {fi,fj}=0,i,j=1,…,n

with the independent differentials at each cotangent space , .

If is a regular value of , then the corresponding level is a smooth -dimensional Lagrangian submanifold of . Geometrically this means that locally around the regular value the map collecting the integrals of motion is a Lagrangian fibration, i.e. it is locally trivial and the fibers are Lagrangian submanifolds.

Let us consider a finite-difference equation

 Y(Pk−ℓ,…,Pk,…Pk+m;k)=0,k,ℓ,m∈Z. (1.1)

relating points of submanifold . Ordinary finite difference equations of this type can be viewed as a dynamical -point map, see [15].

Because all points in (1.1) belong to the given Lagrangian submanifold we may suppose that the corresponding map preserves functions and symplectic form up to a scaling factor. Finite-difference equation sharing integrals of motion with the continuous time system and the symplectic structure is the so-called exact discretization of integrable systems. Nowadays, refactorization in the Poisson-Lie groups is viewed as one of the most universal constructions of finite difference equations (1.1), see discussion in [3, 4, 15, 18, 19, 27] and references within.

The idea is to identify points in (1.1) with intersection points of with auxiliary curve . If and are algebraic, then we can consider the standard equation for their intersection divisor

 div(X⋅Y)=0

as the finite-difference equation (1.1) for the corresponding completely integrable system. Here is the intersection divisor of two algebraic varieties and is a suitable equivalence relation [6, 11, 13, 16]. Our main objective is to study properties of such -point finite-difference equations for different integrable systems [32, 33, 34, 35, 36]. In this paper we restrict ourselves by consideration of cubic and quintic nonlinear Duffing oscillators in order to clarify our view point on relations between the exact discretizations and the intersection divisors.

Thus, we consider integrable systems on two-dimensional plane with a pair coordinates and symplectic form . Because any smooth curve on the plane is a Lagrangian submanifold, we can directly apply classical intersection theory [2, 12] to exact discretization of one-dimensional Hamiltonian systems with the algebraic Hamilton function. Below we consider Hamiltonians associated with hyperelliptic curves on the projective plane defined by equation

 X:y2=a2g+2x2g+2+a2g+1x2g+1+⋯+a1x+a0,aj∈C,

at and various intersections of with the line, quadric and cubic on the plane defined by

 Y:y=P(x),P(x)=bmxm+⋯+b0,m=1,2,3. (1.2)

From now on and mean coordinates on the projective plane, whereas and are coordinates on the phase space .

## 2 Cubic oscillator

The Hamilton function

 H(q,p)=p2−a4q4−a3q3−a2q2−a1q (2.1)

and canonical Poisson bracket determine Hamiltonian equations

 ˙q=∂H∂p=2p,˙p=−∂H∂q=4a4q3+3a3q2+2a2q+a1 (2.2)

and equation of motion

 ¨q=8a4q3+6a3q2+4a2q+2a1, (2.3)

for the generalized oscillator with the cubic nonlinearity [22].

At this integrable system is called a cubic Duffing oscillator without forcing. Duffing oscillators have received remarkable attention in recent decades due to the variety of their engineering applications. For instance magneto-elastic mechanical systems, large amplitude oscillations of centrifugal governor systems, nonlinear vibration of beams, plates and fluid flow induced vibration, seismic waves before earthquake, ecology or cancer dynamics, financial fluctuations and so on are modeled by the nonlinear Duffing equations.

In the numerical integration of nonlinear differential equations, discretization of the nonlinear terms poses extra ambiguity in reducing the differential equation to a discrete difference equation. For instance, in the framework of the standard-like discretization differential equation

 ¨q+Aq+Bq3=0 (2.4)

can be transformed to the finite-difference equation

 qn+1−2qn+qn−1h2+Aqn+B(qn+1+qn−1)q2n=0,

where is a discrete time interval [23, 24]. This equation may be reduced to the expression with the mapping function

 qn+1−2qn+qn−1=F(qn),F(qn)=−Aqn−Bq3nh−2+1/2Bq2n

or to the area preserving map on the plane

 pn+1=pn+ϕ(qn),qn+1=qn+pn+1,

where is a rational control function [20, 21]. This integrable map admits the invariant integral

 ~H=p2n+1+Aqnqn+1+Bq2nq2n+1,

see details in [20, 21, 23, 24, 25, 26, 27], but it is not exact discretization of the Duffing oscillator, i.e. trajectories of the discrete flow do not coincide with the trajectories of the continuous flow.

### 2.1 Exact discretization and intersection divisors

In order to get exact discretization sharing integrals of motion with the continuous time system and the Poisson bracket we can start with the well-known analytical solutions of the Duffing equation (2.4), which are expressed via Jacobi elliptic functions.

Indeed, let us consider the equation (2.4) with initial condition

 q(0)=α,˙q(0)=0.

For and periodic solution is

 q(t)=αcn(2(A+2α2B)1/2t;m),m=α2BA+2α2B.

For and periodic solution reads as

 q(t)=αdn(2B1/2t;m),m=2(1+A2α2).

For and periodic solution has the form

 q(t)=αsn(2(A+α2B)1/2t;m),m=−α2BA+α2B.

Here and are the Jacobi elliptic functions. Discussion of the non-periodic solution can be also found in [23, 24].

Following [3] we can construct exact discretizations of the Duffing equation using these explicit solutions and well-known addition theorems for Jacobi elliptic functions, for instance

 sn(X+Y)=snXcnYdnY+snYcnXdnX1−m2sn2Xsn2Y.

However, it is more easy and convenient to apply standard algorithms of the intersection theory for this purpose.

In order to apply the intersection theory to the exact discretization of the cubic oscillator we put and consider the corresponding level curve on the projective plane defined by equation

 X:y2=f(x),f(x)=a4x4+a3x3+a2x2+a1x+a0. (2.5)

where . Any partial solution and of the Hamiltonian equations (2.2) at is a point of with abscissa and ordinate . It allows us to study relations between points of instead of relations between solutions of the differential equations (2.2).

Let be a smooth nonsingular algebraic curve on a projective plane. Prime divisors are rational points on denoted and is a point at infinity. Divisor

 D=∑miPi,mi∈Z

is a formal sum of prime divisors, and the degree of divisor is a sum deg of multiplicities of points in support of the divisor. Group of divisors is an additive Abelian group under the formal addition rule

 ∑miPi+∑niPi=∑(mi+ni)Pi.

Two divisors are linearly equivalent

 D≈D′

if their difference is principal divisor

 D−D′=div(ψ)≡0modPrinX,

i.e. divisor of rational function on .

Intersection divisor of with some auxiliary smooth nonsingular plane curve

 div(X⋅Y)=0modPrinX

is equal to zero with respect to the linear equivalence of divisors. It allows us to identify intersection divisor with some finite-difference equation (1.1)

 Y(Pk−ℓ,…,Pk,…,Pk+m;k)=k+m∑i=k−ℓmiPi+∑iniP(k)i=0modPrinX. (2.6)

Here we divide intersection divisor in two parts

 div(X⋅Y)=k+m∑i=k−ℓmiPi+∑iniP(k)i

where prime divisors are parameters of discretization implicitly depending on .

### 2.2 Examples of the intersections

Let us consider the intersection of plane curve (2.5) with a parabola

 Y:y=P(x),P(x)=b2x2+b1x+b0

and the corresponding intersection divisor of degree four, see [2], p.113 or [12], p.166. Following Abel we substitute into (2.5) and obtain the so-called Abel polynomial

 ψ(x)=P(x)2−f(x).

Divisor of this polynomial on coincides with , i.e. roots of this polynomial are abscissas of intersection points and forming support of the intersection divisor .

At one of the intersection points is , see examples in Figure 1.

In this case polynomial is equal to

 ψ(x)=(2b1b2−a3)x3+(2b0b2+b21−a2)x2+(2b0b1−a1)x+b20−a0=(2b1b2−a3)(x−x1)(x−x2)(x−x3).

Equating coefficients of one gets relation between abscissas of the remaining rational points and in support of the intersection divisor

 x1+x2+x3=−2b0b2+b21−a22b1b2−a3. (2.6)

If as in Figure 1a, we can define parabola using the Lagrange interpolation by any pair of points , or . For instance, taking the following pair of points one gets

 P(x)=b2x2+b1x+b0=√a4(x−x1)(x−x2)+(x−x2)y1x1−x2+(x−x1)y2x2−x1,

which allows us to determine as functions on and . Substituting coefficients of into the equation (2.6) we obtain an explicit expression for abscissa as a function of coordinates and

 x3=−x1−x2+ϕ(x1,y1,x2,y2),ϕ=−2b0b2+b21−a22b1b2−a3 (2.7)

If we have a double intersection point, for instance as in Figure 1b, then

 x2=−2x1+ϕ(x1,y1),ϕ=−2b0b2+b21−a22b1b2−a3, (2.8)

where function is defined by due to the Hermite interpolation

 P(x)=b2x2+b1x+b0=√a4(x−x1)2+(x−x1)(4a4x31+3a3x21+2a2x1+a1)2y1+y1.

In modern terms, we consider two partitions of the intersection divisor

 div(X⋅Y)=(P1+P2)+P3+P∞anddiv(X⋅Y)=(2P1)+P2+P∞.

Using brackets we separate a part of the intersection divisor which is necessary for polynomial interpolation of auxiliary curve . Because

 div(X⋅Y)=0,

these partitions can be rewritten as addition and doubling of prime divisors

 −P3=P1+P2,−P2=2P1,

where we use standard hyperelliptic inversion , see Figure 1.

At support of the intersection divisor consists of four rational points up to multiplicity. Let us consider the following partitions of this divisor

 div(X⋅Y)=(P1+P2+P3)+P4,div(X⋅Y)=(2P1+P2)+P3,div(X⋅Y)=(3P1)+P2,

see Figure 2. In the first case parabola is defined by the Lagrange interpolation using three ordinary points and . In the second and third cases parabola is defined by the Hermite interpolation using either double and ordinary points or one triple point , respectively.

In the first case abscissa of the fourth intersection point is

 x4=−x1−x2−x3+φ(x1,x2,x3,y1,y2,y3),φ=−a3−2b1b2a4−b22, (2.9)

where function is defined using coefficients of quadratic polynomial

 P(x)=(x−x2)(x−x3)y1(x1−x2)(x1−x3)+(x−x1)(x−x3)y2(x2−x1)(x2−x3)+(x−x1)(x−x2)y3(x3−x1)(x3−x2). (2.10)

In the second case expression for the abscissa looks like

 x3=−2x1−x2+φ(x1,x2,y1,y2),φ=−a3−2b1b2a4−b22. (2.11)

Here function is defined via coefficients of the same polynomial and Hermite interpolation formulae

 P(x)=(x−x1)2y2−(x−2x1+x2)(x−x2)y1(x1−x2)2+(x−x1)(x−x2)(4a4x31+3a3x21+2a2x1+a1)2y1(x1−x2)

In the third case, when we consider tripling the prime divisor on

 (x2,y2)=3(x1,y1),

second abscissa is equal to

 x2=−3x1+φ(x1,y1),φ=−a3−2b1b2a4−b22, (2.12)

where function is defined via coefficients of the polynomial

 P(x)=b2x2+b1x+b0=−(x−x1)2(4a4x31+3a3x21+2a2x1+a1)28y31+(x−x1)(x(6a4x21+3a3x1+a2)−2a4x31+a2x1+a1)2y1+y1. (2.13)

At we have the intersection divisor of with line , which can be represented in the following form

 div(X⋅Y)=(P1+P2)+P3+P4.

It means that line is interpolated by two points and

 P(x)=b1x+b0=x−x2x1−x2y1+x−x1x2−x1y2,

whereas abscissas of remaining two points and are the roots of polynomial

 ψ(x)(x−x1)(x−x2)=a4x2+(a4(x1+x2)+a3)x+a4(x21+x1x2+x22)+a3(x1+x2)+a2−b21.

Thus, are algebraic functions on coordinates and

 x3,4=−x12−x22−a32a4±√α(x1,x2,y1,y2)2a4 (2.14)

where

 α(x1,x2,y1,y2)=4a4((y1−y2x1−x2)2−a2)+a23−2a4a3(x1+x2)−a24(3x21+2x1x2+3x22). (2.15)

In the generic cases, using intersection divisors of plane curve with auxiliary curves

 Y:y=bNxN+bN−1xN−1+⋯+b0,N=1,2,3,…

we can describe multiplication of the prime divisor on any integer , which is a key ingredient of the modern elliptic curve cryptography, and other configurations of the prime divisors entering into the intersection divisor.

All the relations between abscissas (2.6-2.14) are well known, here we only repeat the fairly simple calculations based on Abel’s ideas and their geometric interpretation proposed by Clebsch, see the historical comments in [16]. The modern intersection theory gives a common language for the compact description of these partial cases of intersections [6, 11], whereas modern cryptography equips us with the effective algorithms for such computations [14].

### 2.3 Examples of finite-difference equations

Our aim is to interpret well-studied relations between prime divisors as the finite-difference equations (1.1) realizing various exact discretizations of the given Hamiltonian system and to study the properties of the corresponding discrete maps. For this purpose, we will identify partial solutions and of the Hamilton equations (2.2)

 ˙q=2p,˙p=4a4q3+3a3q2+2a2q+a1

at with a prime divisor , where and .

For instance, substituting

 x1=q1,y1=p1,x2=q2,y2=p2,x3=q3,y3=p3

in (2.6) and one gets finite-difference equations

 q1+q2+q3=ϕ(q1,p1,q2,p2),√a4(q1−q2)(q2−q3)(q3−q2)=q3(p1−p2)+p3(q1−q2)+q1p2−q2p1, (2.16)

where is the rational function on variables and

 ϕ1=a2−a4(q21+4q1q2+q22)+2√a4(q1p1−2(q1p2−q2p1)−q2p2)q1−q2−(p1−p2)2(q1−q2)2ϕ2=2√a4(p1−p2q1−q2−√a4(q1+q2))−a3.

We can directly verify the following properties of the corresponding discrete mapping.

###### Proposition 1

Relations (2.16 ) determine 3-point mapping

 (q1,q2p1,p2)→(q3p3),

preserving the form of Hamiltonain (2.1) and Poisson bracket, i.e. from , and (2.16) will follow that  .

In order to get an iterative system of finite-difference equations we identify a part of abscissas of the intersection points with arbitrary numbers

 xi=λik,yi=μik=±√f(λik),λik∈C.

In this case finite-difference equations (2.6) implicitly depend on the independent variable via parameters of discretization . For instance, addition of prime divisors

 P3=P1+P2,

at

 x1=qk,y1=pk,x3=qk+1,y3=pk+1andx2=λk,y2=μk

determines the following iterative system of 2-point invertible mappings

 qk+1=−qk−λk+ϕ(qk,λk),pk+1=−(b2q2k+1−b1qk+1−b0), (2.17)

where is given by (2.7) and

 b2=√a4,b1=−√a4(qk+λk)+qk−μkqk−λk,b0=√a4qkλk+qkμk−λkpkqk−λk.

Here are arbitrary numbers, whereas the corresponding ordinates

 μk=±√a4λ4k+a3λ3k+a2λ2k+a1λk+H,H=p2k−a4q4k−a3q3k−a2q2k−a1qk

are the functions on the phase space . We have to use this fact to calculate Poisson bracket between variables and (2.17), obtained from variables and .

###### Proposition 2

Relations (2.17) determine iterative system of 2-point invertible mappings

 ⋯−−−→λk−2(qk−1pk−1)−−−→λk−1(qkpk)−→λk(qk+1pk+1)−−−→λk+1⋯

preserving the form of Hamilton function

 H=p2k−1−a4q4k−1−a3q3k−2−a2q2k−1−a1qk−1=p2k+1−a4q4k+1−a3q3k+1−a2q2k+1−a1qk+1=p2k+1−a4q4k+1−a3q3k+2−a2q2k+1−a1qk+1⋯ (2.18)

and Poisson bracket, i.e. from and (2.17) will follow that .

The proof is a straightforward calculation.

Substituting

 x1=qk,y1=pk,x2=qk+1,y2=pk+1

in (2.8) and (2.12) one gets two other iterative systems of 2-point mappings

 qk+1=−qk+ϕ(qk),pk+1=−(b2q2k+1−b1qk+1−b0),

associated with multiplication of prime divisor on integer at . For the cubic Duffing oscillator at we present these mapping explicitly

 N=2,qk+1=p2k−2a4q4k−a2q2k2√a4qkpk,pk+1=q4k(2a4q2k+a2)2−(4a4q4k+p2k)p2k4√a4q2kp2k,N=3,qk+1=qk+4qkp2k(2a4q4k+a2q2k−p2k)4a24q8k+4a2a4q6k+a22q4k−8a4q4kp2k−2a2q2kp2k+p4k,pk+1=−pk+−(qk−qk+1)(a2(qk+1+qk)+2a4q2k(3qk+1−qk))2pk+q2k(qk−qk+1)2(2a4q2k+a2)22p3k (2.19)
###### Proposition 3

Relations (2.19) define two iterative systems of 2-point maps

 ⋯−→N(qk−1pk−1)−→N(qkpk)−→N(qk+1pk+1)−→N⋯

which can be considered as the counterparts of usual geometric progression. These 2-points maps are canonical transformations of valence preserving the form of Hamiltonian (2.18), i.e. from and (2.19) will follow that .

The proof is a straightforward calculation.

Let us now take intersection divisor

 div(X⋅Y)=(P1+P2+P3)+P4,

see Fig.2a. If we identify coordinates of all the intersection points with the partial solutions of the Hamilton equations, relations (2.9) and define 4-point mapping

 (q1,q2,q3p1,p2,p3)→(q4p4)

which has the standard properties.

###### Proposition 4

Discrete map

 q4=−q1−q2−q3+φ(q1,q2,q3,p1,p2,p3),p4=−(b2q24−b1q4−b0),

where and are given by (2.9,2.10), preserves the form of Hamiltonian and original Poisson bracket.

The proof is a straightforward calculation.

In order to get iterative systems of finite-difference equations we identify one of the intersection points with parameter of discretization. For instance, we can substitute

 x1=qk−1,y