Research projects

Splitting and composition methods

Some of the most efficient schemes in Geometric Numerical Integration are based on the related ideas of splitting and composition. If the vector field of the differential equation can be decomposed (split) into two or more parts so that each subproblem can be exactly solved, then a splitting method, possibly of high order, can be obtained by composing these solutions with appropriately chosen weights. If the system to be integrated has a qualitative property which is important to preserve under discretization (e.g., it is a Hamiltonian system, so that the exact solution is symplectic) and each part also has this property, then the numerical approximation furnished by the splitting method, by construction, shares this feature with the exact solution. Another possibility to construct high-order schemes preserving qualitative properties is to take a low-order geometric integrator and compose it with its adjoint, again with appropriate coefficients. In this way, one has composition methods. In any case, the overall scheme still preserves the desirable properties the basic method has, as long as the geometric property is preserved by composition. Splitting and composition methods are essentially equivalent when the vector field can be separated into two parts. There are, however, many instances where the system can be split naturally into three or more solvable parts, and it is in this context where new and more efficient composition methods can play a relevant role. Methods of this class are also very useful for problems originating from partial differential equations of evolution, such as the time dependent Gross-Pitaevskii and the complex Ginzubrg-Landau equations. In this last case, the use of splitting methods with complex coefficients constitutes an interesting alternative. As a matter of fact, we have initiated a systematic analysis of methods of this class, with special emphasis on their preservation properties.

References

Computational problems in celestial mechanics and astrodynamics

The numerical computation of the evolution of a planetary system formed around a star, such as the Solar System, or the study of the dynamics of certain extrasolar planets are just two current problems of great interest in Celestial Mechanics requiring efficient numerical techniques for their treatment. A related domain in this area is the high precision simulation of the trajectories of artificial satellites and other bodies (such as space debris) that are orbiting around the Earth. The number of relevant bodies orbiting around the Earth is greater than the number of sensors that are available for tracking their trajectories. As a consequence, the efficient and precise propagation of orbits of such objects is essential to guarantee the space security. In both cases, the high precision numerical integration of perturbed Keplerian orbits is required. The special structure of the involved differential equations makes it possible to design specially tailored integrators that perform substantially better for the numerical simulation for long (and even short) time spans than existing integrators. In this area we have developed new time-regularization techniques for N-bocy problems with close encounters, as well as new integrators for orbit propagation by means of high order implicit Runge-Kutta methods with Gauss-Legendre collocation nodes.

References

Exponential approximations in quantum mechanics

One of the salient features of quantum mechanics is the fact that the evolution operator, solution of the Schrödinger equation, is unitary. A particular consequence is the preservation of the total transition probabilities between states. Standard perturbation theory, however, fails to provide unitary approximations, and therefore the qualitative description it renders is fundamentally flawed. For this reason, several alternative formulations were proposed in the physics literature since the 1950s to remedy this situation, particularly since the publication by Magnus of his seminal paper. We can mention also in this context the pioneering contributions by Fer and Wilcox. In this approach, the evolution operator is approximated by the exponential or product of exponentials of linear combinations of iterated integrals involving nested commutators of the Hamiltonian operator (or coefficient matrix) evaluated at different times. This proposal has the attractive property of leading to approximate solutions which exhibit, at any order of approximation, the qualitative physical properties of the exact solution, much in the spirit of Geometric Numerical Integration.
In this area we have carried out several studies trying to clarify the algebraic structure of different expansions, and formulating a unifying framework for obtaining all of them. In addition, we have proposed new integration schemes incorporating into their formulation the special features of the problems involved.
On the other hand, and since all these schemes involve in one way or another the evaluation of the exponential of a matrix (or its action on a gven vector), we have also proposed new and optimized schemes for this problem, leading to better efficiencies than the procedures incorporated in commercial suites such as Matlab.

References

Monte Carlo simulations

Sampling from a probability distribution, possibly in spaces of very high dimension, by the Hybrid Monte Carlo (HMC) algorithm is extremely useful in molecular dynamics, Bayesian in statistics,due to its potential for a quick exploration of the state space. Most computational work in HMC goes into integrating numerically an auxiliary Hamiltonian dynamics and therefore improving the efficiency of the numerical integrator is of key importance. The resulting Hamiltonian can be decomposed into two parts, and thus splitting methods constitute a natural option in this setting. Our group has done some work in this area, especially by introducing new numerical schemes combining the ideas of splitting and processing. Another line of much interest consists in developing new integrators for distributions that may be written as a perturbation of a distribution easily sampled with, for instance a Gaussian. This kind of structure arises in many applications, for instance in sampling from distributions defined on a functional space. Also to develop new adaptive splitting integration schemes for HMC is very promising, as well as the extension of such schemes to statistical applications, together with the development of higher stage adaptive integrators.
One essential drawback of Hybrid Monte Carlo is its low acceptance probability when sampling high dimensional spaces. This hampers its application to simulation of very large realistic systems. The modified Hamiltonian Monte Carlo (MHMC) methods, i.e., importance sampling methods that use modified Hamiltonians within a Hybrid Monte Carlo (HMC) framework, often outperform in sampling efficiency HMC, especially for high dimensional problems. Such methods take advantage of two properties of modified Hamiltonian in order to enhance the sampling efficiency of HMC. First, the closeness of modified Hamiltonian to true Hamiltonian guarantees low variability among importance weights, implying faster convergence to the target distribution. Second, the symplectic integrator preserves modified Hamiltonian better than it preserves true Hamiltonian that leads to a more favorable value of the acceptance probability in the MHMC algorithms.