Geometric properties of Kahan’s method

# Geometric properties of Kahan’s method

## Abstract

We show that Kahan’s discretization of quadratic vector fields is equivalent to a Runge–Kutta method. When the vector field is Hamiltonian on either a symplectic vector space or a Poisson vector space with constant Poisson structure, the map determined by this discretization has a conserved modified Hamiltonian and an invariant measure, a combination previously unknown amongst Runge–Kutta methods applied to nonlinear vector fields. This produces large classes of integrable rational mappings in two and three dimensions, explaining some of the integrable cases that were previously known.

## 1 Introduction: Kahan’s method for quadratic vector fields

Consider a system of differential equations arising from a quadratic vector field

 ˙x=f(x):=Q(x)+Bx+c,x∈Rn, (1)

where is an -valued quadratic form, , and . Consider the numerical integration method with step size given by

 x′−xh=Q(x,x′)+12B(x+x′)+c (2)

where

 Q(x,x′)=12(Q(x+x′)−Q(x)−Q(x′)) (3)

is the symmetric bilinear form obtained from the quadratic form by polarisation. We call (2) the Kahan method. It is symmetric (i.e., self-adjoint), and, crucially, it is only linearly implicit, that is, can be computed by solving a single linear system (because the right hand side of (2) is linear in ). The method (2) was introduced by Kahan in [8] for two examples, a scalar Riccati equation and a 2-dimensional Lotka–Volterra system ([8], p. 14) and written down in the general form (2) in [9] (see also references therein).

Because of the different treatment of each term and the unusual treatment of the quadratic term, Kahan called (2) an ‘unconventional’ method.

The map obtained from applying the Kahan method to various quadratic vector fields has been shown to be completely integrable in a number of cases (see [5], [6], [10], and references therein). In most cases, the conserved quantities depend on the step size . At present there is no single ‘integrability mechanism’ known which accounts for all integrable cases.

In this paper we show that the Kahan method is a Runge–Kutta method. As such it shares a number of features with all Runge–Kutta methods: it has a B-series, it is affine covariant, and it preserves all affine symmetries and all linear integrals of automatically. As a symmetric linear method it preserves all affine reversing symmetries of automatically, and the B-series of its modified vector field contains only even powers of .

We then consider the case that is a Hamiltonian vector field on either a symplectic vector space or a Poisson vector space with constant Poisson structure, in any dimension . We show that in this case the Kahan map has a conserved quantity that converges to the Hamiltonian of the vector field as . It also has a conserved measure which converges to the Euclidean measure as . These general properties explain some of the integrable cases considered in [10].

## 2 Kahan’s method as a Runge–Kutta method

###### Proposition 1.

The Kahan method coincides with the Runge–Kutta method

 x′−xh=−12f(x)+2f(x+x′2)−12f(x′) (4)

###### Proof.

We have

 (5)

(Many other Runge–Kutta methods also coincide with the Kahan method when restricted to quadratic vector fields. In this paper, we restrict our attention to (4).) As already noted by Kahan [9], the Kahan method also coincides with a certain Rosenbrock method on quadratic vector fields, for expanding in Taylor series about gives

 (6)

so

 x′−xh=(I−h2f′(x))−1f(x). (7)

From the symmetry of the method, or by expanding instead around , Kahan’s method can also be written

 x′−xh=(I+h2f′(x′))−1f(x′). (8)

The B-series of this method is

 x′=x+∞∑k=0hk+12kf′(x)kf(x),

that is, it contains only tall trees. (For nonquadratic vector fields, the methods (7) and (4) are not necessarily equivalent.)

The Runge–Kutta method (4) has 3 stages and Butcher tableau

0 12 1 0 0 0 −14 1 −14 −12 2 −12 −12 2 −12

The modified vector field of the Kahan method applied to quadratic vector fields can be calculated using standard methods [4]. Its first few terms are

A calculation using conjugation by B-series, considering only quadratic vector fields, now yields the following result. We omit the details.

###### Proposition 2.

Kahan’s method applied to general quadratic fields has order 2 and is conjugate to symplectic up to order 4. It is not conjugate by B-series to a method of order greater than 2 or conjugate-symplectic by B-series to order higher than 4.

## 3 Conservative properties of Kahan’s method

We now consider the conservative properties of the Kahan method in the case of canonical Hamiltonian systems where is the Hamiltonian or energy of the system. First, note that the method (4) is the member of the class of Runge–Kutta methods

 x′−xh=af(x)+(1−2a)f(x+x′2)+af(x′). (9)

These are all symmetric, A-stable, and second order. Some other members of this family are also known to have conservative properties:

1. When , we have the midpoint rule. It is symplectic for canonical Hamiltonian systems. Because it is symplectic, it conserves the Euclidean measure. When the Hamiltonian is analytic, the method has a formal invariant . When is quadratic (i.e. when is linear) this series converges to give a conserved quantity of the method.

2. When , we have the trapezoidal rule. It is conjugate to the midpoint rule (the conjugacy being half an Euler step), and so it is also conjugate to symplectic and hence conserves a measure close to the Euclidean measure, and it also has a formal invariant close to .

3. When , we have ‘Simpson’s method’ [2], so-called because the right hand side of (9) is Simpson’s quadrature of , appearing in the average vector field method , which conserves the Hamiltonian in canonical Hamiltonian systems. Simpson’s method preserves quartic Hamiltonians exactly because it coincides with the average vector field method in that case. It is not conjugate to symplectic in the sense of B-series [4].

###### Proposition 3.

Kahan’s method has a conserved quantity given by the modified Hamiltonian

 ˜H(x):=H(x)+13h∇H(x)T(I−12hf′(x))−1f(x) (10)

for all cubic Hamiltonian systems on symplectic vector spaces and on all Poisson vector spaces with constant Poisson structure. The modified Hamiltonian is (i) a rational function of ; (ii) an even function of ; and (iii) given by a convergent series of elementary Hamiltonians containing only even-order tall trees.

###### Proof.

We first consider the homogeneous case, i.e., we let where is an arbitrary (not necessarily invertible) constant antisymmetric matrix and where is a symmetric trilinear form. Note that for all , . For any of the methods (9), we have (writing )

The case is Simpson’s method, confirming that is conserved in that case.

For Kahan’s method, , and we have from Eqs. (7,8) that

Therefore

 (11)

because each matrix in square brackets is antisymmetric. (Each term is an elementary Hamiltonian corresponding to a superfluous tall tree.) The expression is a rational function of and so if its Taylor series in is zero, the function is zero. Therefore the Kahan method has a first integral . This can also be written in the symmetric form

 C(x,x,x′)=(C(x,x,x′)+C(x,x′,x′))/2=C(x,(x+x′)/2,x′)

or explicitly as a function of as

As in (11), for all , so we can also write in a form manifestly even in ,

 ˜H(x)=H(x)+16h2∇H(x)T(I−12h2f′(x)2)−1f′(x)f(x).

Now consider the case that is cubic but not homogeneous. We extend it to a homogeneous function so that , and extend to by adding a zero initial row and column, so that . The linear integral is conserved by Kahan’s method, and , so the modified Hamiltonian of Kahan’s method for reduces to a modified Hamilton of Kahan’s method for given by the same formula (10) as in the homogeneous case. ∎

The Kahan map and its conserved quantity are rational functions of whose degrees are described in the following proposition. When has full rank, for (resp. ), is degree 3 over degree 2 (resp. degree 5/2, degree 5/4) and the Kahan map is degree 2/2 (resp. degree 3/2, degree 4/4). In the planar case, Kahan’s method gives a rational map with cubic invariant curves. We conjecture that the dynamics of the Kahan map in two dimensions is related to the abelian group structure of elliptic curves (and in higher dimensions, to that of abelian varieties) as is the case for planar QRT maps [3].

###### Proposition 4.

Let be a cubic in and let be a rank antisymmetric matrix.

1. The degree of the denominator of is at most and the degree of the numerator of is at most . When the degree of the numerator of is at most .

2. The degree of the denominator of the Kahan map is at most and the degree of the numerator is at most . When the degree of the numerator is at most .

###### Proof.
1. Because the method is linearly covariant we can assume without loss of generality that is in its Darboux normal form

 K=(L000),L=(0I−I0)

where is . Numbering the blocks of 1 and 2, the denominator of is equal to

Each entry of the matrix is linear in so the final determinant has degree at most . Next, we write

where is the adjoint of . The first term in the numerator has degree at most , and have degree at most 2, and where each entry in is the determinant of a matrix whose entries are linear in . Hence the degree of the numerator of is at most .

Finally we consider the case . Since is even, must be even. First consider the case that is a homogeneous cubic. Then and . Thus in this case we have

Expanding the matrix inverse using Cramer’s rule now shows that the degree of the numerator is at most . The terms of degree come from terms in the minors of of degree in . Every in the matrix is multiplied by , thus these terms also have degree in . However, is an even function of and so these terms must sum to zero. Thus the degree of the numerator is at most .

When is a nonhomogeneous cubic, the terms of degree in the numerator of come from the cubic terms in only, and thus vanish as in the homogeneous case. The terms of degree are odd in and hence vanish as before.

2. The proof for the general case follows as above. For the case , we first consider the case that is a homogeneous cubic. Then and the Kahan map can be written

 x′=x+(I−h2f′(x))−1f(x)=(I−h2f′(x))−1x.

Expanding the matrix inverse using Cramer’s rule now shows that the degree of the numerator is at most . In the nonhomogeneous case, the terms of degree in the numerator come from the cubic terms in only, and thus vanish as in the homogeneous case.

Examples suggest that there are no other values of or other than which lead to a reduction in degree.

###### Proposition 5.

Kahan’s method preserves the measure

 dx1∧⋯∧dxndet(I−12hf′(x))

for all cubic Hamiltonians on symplectic vector spaces and on Poisson vector spaces with constant Poisson structure.

###### Proof.

Let be the Jacobian of the Kahan mapping. Differentiating the mapping (9) gives

because is linear in . Solving for ,

 A=(I−(14+12a)hf′(x′)−(14−12a)hf′(x))−1(I+(14+12a)hf′(x)+(14−12a)hf′(x′)).

Now where is symmetric. From Sylvester’s determinant theorem, . The sum of such matrices has the same property, so

 detA=det(I−(14+12a)hf′(x)−(14−12a)hf′(x′))det(I−(14+12a)hf′(x′)−(14−12a)hf′(x)).

This yields invariant measures in 3 cases:

1. when (Kahan’s method), ;

2. when (midpoint rule), ;

3. when (trapezoidal rule), .

By integrable symplectic map we adopt the definition of Bruschi et al. [1]: a symplectic map on a -dimensional symplectic manifold is integrable if it has functionally independent integrals in involution. We will say that leaf-preserving Poisson maps are integrable if the the map is integrable on each leaf.

###### Corollary 6.

Kahan’s method yields an integrable mapping of the plane when applied to any canonical Hamiltonian system in the plane with cubic Hamiltonian. Kahan’s method yields an integrable mapping of when applied to any Poisson system on with constant Poisson structure and any cubic Hamiltonian.

###### Proof.

A measure and a first integral are sufficient for integrability in the plane. The odd-dimensional case with constant has a linear Casimir which is conserved by the method, reducing the situation in this case to two dimensions on each level set of the Casimir. ∎

###### Corollary 7.

When and is a homogeneous cubic, and Kahan’s method preserves the -independent measure .

###### Proof.

We have

Any map that preserves a measure and an integral also preserves the measure . Taking and gives the -independent measure . ∎

## 4 Discussion

Level sets of are shown in Figure 1 for . Notice that the separatrices persist (and are unchanged) for all , but that the singular set moves in from infinity as increases and alters the topology of the level sets. For the topology of the bounded orbits is unaltered.

The bounded orbits of Figure 1 are symmetric, so for the following numerical experiments we used , which has bounded, nonsymmetric orbits, an elliptic fixed point at , and a separatrix meeting . Level sets of for this case are shown in Figure 2.

Numerical experiments strongly indicate that the following observations hold.

1. No other method of the family (9) has a modified Hamiltonian when is cubic, apart from the known cases , , and (see Figure 3).

2. The midpoint and trapezoidal rules do not have a first integral for all cubic (even though they do have a formal invariant close to ) (see Figure 4).

3. Simpson’s method is not measure-preserving for all cubic . (For , a numerical calculation finds eigenvalues of periodic orbits, with , contradicting measure preservation.)

4. Kahan’s method does not preserve any symplectic form in dimension for all cubic . (A numerical calculation of periodic points finds eigenvalues that do not occur in , pairs. Proposition 2 establishes this for a limited class of symplectic forms.)

5. Compositions of Kahan’s method with different step sizes do not have a modified Hamiltonian when is cubic (see Figure 5).

Our results are significant and novel for the study of both the integrability of the mappings produced by Kahan’s method and for the study of the geometric properties of Runge–Kutta methods:

• First, our results explain the integrability of the map obtained when Kahan’s method is applied to some of the examples of [10]: their Eq. (4.2) (); Eq. (5.4) (); Eq. (8.1) (Volterra chain in , , constant , integral in Eq. (8.6) is a function of our and the Casimir ); Eq. (9.1) (Dressing chain in , , constant ). Our results explain the invariant measure and cubic integral of their Eq. (11.1) (three wave system in , ), the invariant measure for the family of systems in their Prop. 1, and the linear integrals throughout [10].

• Second, our results (e.g. Corollary 6) systematically produce new integrable cases of Kahan’s method.

• Third, we have shown that Kahan’s method in dimension 4 and greater provides examples of maps with nonlinear integrals and conserved measures unrelated (in general) to integrability or obvious symmetries, again a novel feature. For example, our results imply that Kahan’s application of the method to the Korteweg–de Vries equation in [9] preserves a measure and a modified energy (but the higher order compositions of the method in [9] probably do not).

• Fourth, we have shown that Kahan’s method has novel properties previously unknown amongst Runge–Kutta methods, indeed amongst all B-series. It is known that B-series methods cannot conserve the measure even for linear vector fields [7]; Kahan’s method circumvents this by conserving a modified measure. It is a novel conjugate-to-energy preserving method for cubic . In the plane, it is also conjugate to symplectic. Thus, while no conjugate to symplectic methods are known that are also energy preserving in general, here we have one that preserves at least a modified energy, and preserves it exactly (not merely as a formal invariant).

On the other hand, there are open questions in all of these areas. While it was already suggested in [10] that there could be an underlying ‘integrability mechanism’ unifying the integrable cases, here we have unified only some of these. In addition the hoped-for unification should now be extended to include non-integrable cases preserving a measure and/or some integrals as well. On the numerical side, it is not known precisely which Runge–Kutta or B-series methods share the properties of Kahan’s method, if any are higher order integrators, or if any are conservative for nonquadratic (e.g. other polynomial) vector fields.

## Acknowledgements

This research was supported by a Marie Curie International Research Staff Exchange Scheme Fellowship within the 7th European Community Framework Programme, and by the Marsden Fund of the Royal Society of New Zealand and the Australian Research Council. We would like to thank the referees for their careful reading of the manuscript.

## References

### References

1. Bruschi M, Ragnisco O, Santini PM, and Gui-Zhang T 1991, Integrable symplectic maps, Physica D 49 273–294.
2. Celledoni E, McLachlan RI, McLaren DI, Owren B, Quispel GRW, and Wright W 2009, Energy-preserving Runge–Kutta methods, Mathematical Modelling and Numerical Analysis 43 645–649.
3. Duistermaat JJ 2010, Discrete Integrable Systems: QRT Maps and Elliptic Surfaces, Springer, Heidelberg.
4. Hairer E, Lubich C, and Wanner G 2006, Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations, 2nd ed., Springer, Berlin.
5. Hirota R and Kimura K 2000, Discretization of the Euler top, J. Phys. Soc. Jap. 69 627–630.
6. Hone ANW and Petrera M 2009, Three-dimensional discrete systems of Hirota–Kimura type and deformed Lie–Poisson algebras, J. Geom. Mech. 1 55–85.
7. Iserles A, Quispel GRW, and Tse PSP 2007, B-series methods cannot be volume-preserving, BIT 47 351–378.
8. Kahan W 1993, Unconventional numerical methods for trajectory calculations, Unpublished lecture notes.
9. Kahan W and Li R-C 1997, Unconventional schemes for a class of ordinary differential equations—with applications to the Korteweg–de Vries equation, J. Comput. Phys. 134 316–331.
10. Petrera M, Pfadler A, and Suris YB 2011, On integrability of Hirota–Kimura type discretizations, Regular and Chaotic Dynamics 16 245–289.
264136