%format:amsplain \input brimacs \input fonts \input footmacs %\baselineskip=16pt \def\d{\displaystyle} \overfullrule=0pt \TagsOnRight \nopagenumbers \cen{\bf SOLVING ODE's NUMERICALLY} \cen{\bf WHILE PRESERVING A FIRST INTEGRAL} \gsk \cen{G.R.W. Quispel} \cen{School of Mathematics, La Trobe University} \cen{Bundoora, Victoria 3083, Australia.} \cen{(email: R.QUISPEL\@LATROBE.EDU.AU)} \bsk \cen{and} \bsk \cen{H.W. Capel} \cen{Instituut voor Theoretische Fysica, Universiteit van Amsterdam} \cen{Valckenierstraat 65, 1018XE, Amsterdam} \cen{The Netherlands} \gsk \baselineskip=16pt \nin{\bfss Abstract} We give general algorithms for the numerical integration of ordinary differential equations (ODE's) that possess a first integral $I(x)$. Our discrete algorithms preserve the integral $I$ exactly. Our method works both for dissipative and for all Hamiltonian ODE's. For non-Hamiltonian systems a sufficient (but not necessary) condition for our method to work is that $|\nabla I| \not= 0$ on the isosurface we are integrating on. \vfill \eject \pageno=1 \nin{\bfss 1. Introduction} In recent years much effort has been devoted to the construction of numerical integration schemes for ordinary differential equations (ODE's) in such a way that some {\it qualitative} property of the ODE is exactly preserved. This has resulted in symplectic integrators for Hamiltonian ODE's [1-3], volume-preserving integrators for divergence-free ODE's [4-6], and symmetry-preserving integrators for ODE's with symmetries and/or time-reversing symmetries [7-9]. In this Letter we are interested in the numerical integration of an ODE that possesses a first integral, and we explicitly give numerical integration methods that preserve this first integral. Special cases of this problem have been studied before. We mention the work by Cooper on preserving quadratic integrals [10], the work by Itoh and Abe\footnote{We mention these authors in particular because their integrators have the nice property of being free of singularities.} on preserving energy in Hamiltonian systems [11], and the work by Greenspan and collaborators on preserving the first integrals of the $n$-body problems of celestial mechanics [12]. We stress that our method as described in this Letter is not restricted to such special cases, and in particular is not restricted to Hamiltonian systems (although Hamiltonian systems are included as a very special case). In this Letter our treatment is confined to preserving {\it one} first integral. The simultaneous preservation of two or more first integrals will be treated in a future paper [13]. \bsk \nin{\bfss 2. \ Integral-preserving integrators} Consider the autonomous ODE $$ {dx\over dt} = f(x), \qquad\qquad x \in \Bbb R^n, \tag 1 $$ where $f$ is a function from $\Bbb R^n$ to $\Bbb R^n$. We are interested in the case that the ODE (1) has a first integral $I$: $$ {dI(x) \over dt} = 0, \tag 2 $$ where $I$ is a function from $\Bbb R^n$ to $\Bbb R$ (in the sequel we will assume $f$ and $I$ possess suitable regularity properties). In general it is of course impossible to solve the initial value problem for the ODE (1) analytically. We therefore pose the following problem: \ {\it Find a discrete approximation to the ODE (1) that preserves the first integral $I(x)$ exactly.} In this paper this problem is solved\footnote{A sufficient, but not necessary, condition for our method to work is that $|\nabla I| \not= 0$ on the isosurface we are integrating on.}. Our solution proceeds in 4 steps: \item{(i)} We rewrite the ODE (1) in the form of a ``skew-gradient system", i.e. $$ {dx\over dt} = S(x). \nabla I(x), \tag 3 $$ where $S(x)$ is some skew-symmetric matrix (i.e. $S^t = -S$).\footnote{Note that eq.(3) is in general {\it not} Hamiltonian, see example III of section 3.} This can always be done (generally in infinitely many different ways) subject to a mild technical condition. Below we will give an explicit construction of $S(x)$ for a given vectorfield $f(x)$ and integral $I(x)$. \item{(ii)} We use (3) to split\footnote{A general discussion of splitting methods is given in [14].} the $n$-dimensional vectorfield into essentially 2-dimensional vectorfields that all possess the original integral $I$. \item{(iii)} We integrate each 2-dimensional $(2-D)$ vectorfield, using a discrete approximation that preserves the integral $I$ exactly. \item{(iv)} We construct an $n$-dimensional integrator for the ODE (1) by composition of the 2-dimensional integrators obtained in (iii). Since each $2-D$ integrator preserves $I$, so does the ensuing $n$-dimensional integrator. \msk \nin STEP (i): \ WRITE THE ODE IN THE FORM ${dx\over dt} = S.\nabla I$ Combining (1) and (2) it follows that $$ f.\nabla I = 0. \tag 4 $$ We now want to find a skew-symmetric matrix $S$ such that $$ S.\nabla I = f. \tag 5 $$ This will automatically ensure that (4) holds. For given $f$ and $\nabla I$, equation (5) is a system of linear equations for the entries of the matrix $S$. For dimension $n \geq 3$ this system of equations is under-determined. It is easy to check that a particular solution $S^P$ to equation (5) is given by $$ S^P_{i,j} = { f_i {\partial I\over \partial x_j} - f_j {\partial I\over \partial x_i} \over \sum^n_{k=1} ({\partial I \over \partial x_k})^2} \tag 6 $$ The general solution of (5) follows by adding this particular solution to the solution $S^H$ of the homogeneous equation $$ S^H. \nabla I = 0. \tag 7 $$ Assuming $\d{{\partial I\over \partial x_j}}$ is not identically zero\footnote{The case that ${\partial I\over \partial x_j}$ is identically zero for some $j$ is treated in section 4.} for any $j$, the solution of (7) is given by $$ S^H_{i,j} = \alpha_{i,j} \prod\Sb k\endSb^* \left( {\partial I \over \partial x_k}\right) \tag 8 $$ where the product symbol $\d{ \prod\Sb k\endSb^*}$ means $k \not= i, k \not= j$, i.e. $k$ takes on the values $1, \dots, i-1, i+1, \dots, j-1, j+1, \dots, n$. The coefficient functions $\alpha_{i,j}(x)$ in (8) can be chosen arbitrarily for $1 \leq i < j < n$. The remaining $\alpha_{i,j}$ are then given by $$ \align \alpha_{i,n}(x) &= - \sum^{n-1}_{j=1} \alpha_{i,j}(x), \qquad i = 1, \dots, n-1, \\ \alpha_{j,i}(x) &= - \alpha_{i,j}(x). \tag 9 \endalign $$ The general solution $S^G$ to (5) is then obtained as the sum of the particular solution (6) and the homogeneous solution (8)-(9), i.e. $$ S^G = S^P + S^H. $$ STEP (ii): SPLITTING INTO $2-D$ VECTOR FIELDS ALL POSSESSING THE INTEGRAL $I$ For convenience we use here the equivalence of the ODE (1) with the vectorfield $$ \align v &:= \sum_i f_i {\partial \over \partial x_i} \\ &= \sum_{i < j} v_{i,j} \tag 10 \endalign $$where the vector fields $v_{i,j}$ are defined by $$ v_{i,j} := S_{i,j} {\partial I\over \partial x_j} {\partial \over \partial x_i} + S_{j,i} {\partial I \over \partial x_i} {\partial \over \partial x_j}. \tag 11 $$ Using the skew-symmetry of $S$, these vectorfields $v_{i,j}$ can be rewritten $$ v_{i,j} = S_{i,j} \left( {\partial I \over \partial x_j} {\partial \over \partial x_i} - {\partial I \over \partial x_i} {\partial \over \partial x_j} \right). \tag 12 $$ From (12) it follows immediately that each vectorfield $v_{i,j}$ possesses the original integral $I$, i.e. $$ \align v_{i,j} I(x) &= S_{i,j} \left( {\partial I\over \partial x_j} {\partial I \over \partial x_i} - {\partial I\over \partial x_i} {\partial I \over \partial x_j}\right) \\ &= 0. \tag 13 \endalign $$ (It is for this reason that we have insisted that the matrix $S$ in (3) must be skew-symmetric). Note that the above leads to a splitting into a maximum of $\d{\pmatrix n\\2\endpmatrix}$ vectorfields. \msk \nin STEP (iii): INTEGRAL-PRESERVING INTEGRATION OF $2-D$ VECTORFIELDS The vectorfield $v_{i,j}$ (12) is equivalent to the following essentially 2-dimensional ODE: $$ \align {dx_i\over dt} &= S_{i,j} (x) {\partial I \over \partial x_j}, \\ {dx_j\over dt} &= - S_{i,j} (x) {\partial I \over \partial x_i}, \tag 14\\ {dx_k\over dt} &= 0. \qquad\qquad\quad k \not= i, j \endalign $$ (Note that the ODE (14) is ``symplectic in the $x_i, x_j$ plane"). A first-order integral-preserving integrator for the ``2-dimensional" ODE is (cf. [15]) $$ \align x'_i &= x_i + \tau \tilde S_{i,j} (x,x',\tau) {I(x'_i,x'_j)-I(x'_i,x_j) \over x'_j - x_j} \\ \phi_{i,j}(x, \tau): x'_j &= x_j-\tau \tilde S_{i,j}(x,x',\tau) {I(x'_i,x_j) - I(x_i,x_j) \over x'_i - x_i} \tag 15 \\ x'_k &= x_k, \qquad \qquad \qquad\qquad \qquad \qquad\qquad \qquad \qquad k \not= i, j. \endalign $$ A second-order integral-preserving integrator for the ``2-dimensional" ODE (14) is $$ \align x'_i &= x_i + {\tau\over 2} \hat S_{i,j} (x,x',\tau) {I(x'_i,x'_j)- I(x'_i,x_j)+I(x_i,x'_j)-I(x_i,x_j) \over x'_j - x_j} \\ \psi_{i,j}(x,\tau): x'_j &= x_j - {\tau\over 2} \hat S_{i,j} (x,x',\tau) {I(x'_i,x'_j)-I(x_i,x'_j)+I(x'_i,x_j)-I(x_i,x_j) \over x'_i - x_i} \tag 16 \\ x'_k &= x_k, \qquad \qquad \qquad\qquad \qquad \qquad\qquad \qquad \qquad k \not= i,j. \endalign $$ In these expressions for the integrators (15) and (16) we have used the shorthand notation $$ \align I(x'_i,x'_j) &:= I(x_1, \dots, x_{i-1}, x'_i, x_{i+1}, \dots x_{j-1}, x'_j, x_{j+1} \dots, x_n),\\ I(x'_i,x_j) &:= I(x_1, \dots, x_{i-1}, x'_i, x_{i+1}, \dots, x_n),\\ I(x_i,x'_j) &:= I(x_1, \dots, x_{j-1}, x'_j, x_{j+1}, \dots, x_n), \tag 17\\ I(x_i,x_j) &:= I(x_1, \ldots\ldots\ldots\ldots\ldots\ldots, x_n), \endalign $$ and $\tilde S_{i,j}$ and $\hat S_{i,j}$ are approximations to $S_{i,j}$, i.e. they satisfy respectively $$ \tilde S_{i,j}(x, x', \tau) = S_{i,j} (x) + \Cal O(\tau). \tag 18 $$ and $$ \hat S_{i,j}(x,x', \tau) = \hat S_{i,j}(x', x, -\tau) := S_{i,j}(x) + \Cal O(\tau). \tag 19 $$ (This freedom in the choice of $\tilde S_{i,j}$ and $\hat S_{i,j}$ is sometimes useful). It is not difficult to check that both integrators (15) and (16) preserve the integral $I$, i.e. they satisfy $$ I(x') = I(x) \tag 20 $$ The integrator $\psi_{i,j}$ moreover satisfies $$ \chi(x,\tau) = \chi^{-1}(x,-\tau) \tag 21 $$ which makes it suitable as a building block for constructing integrators of arbitrary order, using Yoshida's method [1,16]. \msk STEP (iv): CONSTRUCTION OF $n$-DIMENSIONAL INTEGRAL-PRESERVING INTEGRATOR FROM $2-D$ INTEGRATORS Here we will give two ways to construct $n$-dimensional integral-preserving integrators for the ODE (1). These integrators are constructed from the ``2-dimensional" integral-preserving integrators $\phi_{i,j}$ resp. $\psi_{i,j}$ given in equations (15) and (16). The simplest integral-preserving integrator for the ODE (1) is $$ \phi(x,\tau) := \prod_{i