Section: Research Program

Structure-preserving numerical schemes for solving ordinary differential equations

Participants : François Castella, Philippe Chartier, Erwan Faou.

ordinary differential equation, numerical integrator, invariant, Hamiltonian system, reversible system, Lie-group system

In many physical situations, the time-evolution of certain quantities may be written as a Cauchy problem for a differential equation of the form

For a given $y_0$, the solution $y\left(t\right)$ at time $t$ is denoted ${\varphi }_t\left(y_0\right)$. For fixed $t$, ${\varphi }_t$ becomes a function of $y_0$ called the flow of (1 ). From this point of view, a numerical scheme with step size $h$ for solving (1 ) may be regarded as an approximation ${\Phi }_h$ of ${\varphi }_h$. One of the main questions of geometric integration is whether intrinsic properties of ${\varphi }_t$ may be passed on to ${\Phi }_h$.

This question can be more specifically addressed in the following situations:

Reversible ODEs

The system (1 ) is said to be $\rho$-reversible if there exists an involutive linear map $\rho$ such that

It is then natural to require that ${\Phi }_h$ satisfies the same relation. If this is so, ${\Phi }_h$ is said to be symmetric. Symmetric methods for reversible systems of ODEs are just as much important as symplectic methods for Hamiltonian systems and offer an interesting alternative to symplectic methods.

ODEs with an invariant manifold

The system (1 ) is said to have an invariant manifold $g$ whenever

is kept globally invariant by ${\varphi }_t$. In terms of derivatives and for sufficiently differentiable functions $f$ and $g$, this means that

As an example, we mention Lie-group equations, for which the manifold has an additional group structure. This could possibly be exploited for the space-discretisation. Numerical methods amenable to this sort of problems have been reviewed in a recent paper [56] and divided into two classes, according to whether they use $g$ explicitly or through a projection step. In both cases, the numerical solution is forced to live on the manifold at the expense of some Newton's iterations.

Hamiltonian systems

Hamiltonian problems are ordinary differential equations of the form:

with some prescribed initial values $(p\left(0\right),q\left(0\right))=(p_0,q_0)$ and for some scalar function $H$, called the Hamiltonian. In this situation, $H$ is an invariant of the problem. The evolution equation (4 ) can thus be regarded as a differential equation on the manifold

Besides the Hamiltonian function, there might exist other invariants for such systems: when there exist $d$ invariants in involution, the system (4 ) is said to be integrable. Consider now the parallelogram $P$ originating from the point $(p,q)\in {ℝ}^{2d}$ and spanned by the two vectors $\xi \in {ℝ}^{2d}$ and $\eta \in {ℝ}^{2d}$, and let $\omega (\xi ,\eta )$ be the sum of the oriented areas of the projections over the planes $(p_i,q_i)$ of $P$,

where $J$ is the canonical symplectic matrix

A continuously differentiable map $g$ from ${ℝ}^{2d}$ to itself is called symplectic if it preserves $\omega$, i.e. if

A fundamental property of Hamiltonian systems is that their exact flow is symplectic. Integrable Hamiltonian systems behave in a very remarkable way: as a matter of fact, their invariants persist under small perturbations, as shown in the celebrated theory of Kolmogorov, Arnold and Moser. This behavior motivates the introduction of symplectic numerical flows that share most of the properties of the exact flow. For practical simulations of Hamiltonian systems, symplectic methods possess an important advantage: the error-growth as a function of time is indeed linear, whereas it would typically be quadratic for non-symplectic methods.

Differential-algebraic equations

Whenever the number of differential equations is insufficient to determine the solution of the system, it may become necessary to solve the differential part and the constraint part altogether. Systems of this sort are called differential-algebraic systems. They can be classified according to their index, yet for the purpose of this expository section, it is enough to present the so-called index-2 systems

where initial values $(y\left(0\right),z\left(0\right))=(y_0,z_0)$ are given and assumed to be consistent with the constraint manifold. By constraint manifold, we imply the intersection of the manifold

and of the so-called hidden manifold

This manifold $ℳ={ℳ}_1\bigcap {ℳ}_2$ is the manifold on which the exact solution $\left(y\right(t),z(t\left)\right)$ of (5 ) lives.

There exists a whole set of schemes which provide a numerical approximation lying on ${ℳ}_1$. Furthermore, this solution can be projected on the manifold $ℳ$ by standard projection techniques. However, it it worth mentioning that a projection destroys the symmetry of the underlying scheme, so that the construction of a symmetric numerical scheme preserving $ℳ$ requires a more sophisticated approach.