Explicit Methods for Time Integration with Extended Stability Domain

Abstract

Initial value problems (IVPs) for systems of ordinary differential equations (ODEs) appear in many realworld applications. This thesis is motivated by the need for real-time simulation of large-scale mechanical systems to support control of robots in hardware-in-the-loop environments. Classical IVP solvers include time marching schemes such as the Euler method, available in both explicit and implicit variants. While explicit methods are computationally inexpensive, they suffer from stability issues when applied to stiff systems -- especially during contact between stiff materials. Implicit methods handle stiffness more robustly but require solving large systems of equations at every time step, which is impractical under real-time constraints with high-dimensional systems. To meet these constraints, we investigate explicit methods that scale linearly in computational cost with system size. However, classical explicit solvers require prohibitively small time steps on stiff systems, negating their e#ciency. Therefore, the focus of this thesis is on developing explicit solvers that remain stable at larger time steps. Modern high-performance computing -- especially GPU architectures-favors algorithms that are parallelizable rather than sequential. This motivates the use of multiderivative methods, which allow parallel computation of higher-order time derivatives and are well-suited to modern hardware. Chapter 1 introduces key concepts related to IVPs, stiffness, and solver requirements, and explains how stiffness is characterized through the eigenvalues of the system's Jacobian. Chapter 2 reviews classical solvers such as Runge-Kutta and general linear methods, with a focus on their stability properties. Chapter 3 introduces multiderivative methods such as Taylor and multiderivative Runge-Kutta (MDRK) integrators, and outlines how to compute multiple time derivatives in parallel. MDRK methods can be used to approximate higher time-derivatives, where their calculation is not yet implemented in the model. Chapter 4 presents the central contribution of this thesis: a Krylov subspace-based technique for approximating the system Jacobian with linear computational complexity. Because the Jacobian determines the system's local linear dynamics, this approximation enables the use of advanced methods that incorporate Jacobian information without incurring the full computational cost of evaluating or storing it. We focus on two classes of such methods: the time-accurate stabilized explicit (TASE) methods proposed by Bassenne, Fu, and Mani [2], and the exponential Rosenbrock integrators introduced by Hochbruck, Ostermann, and Schweitzer [15]. In chapter 5, we explore these methods in combination with the Krylov Jacobian approximation, referring to the resulting schemes as Krylov Jacobian methods. The stability properties of Krylov Jacobian methods depend on the accuracy of the Jacobian approximation and differ significantly from classical methods such as Runge-Kutta schemes. They are determined by the global distribution of eigenvalues, rather than their individual positions. Therefore, the notion of the stability domain does not apply to these methods in the same sense as it does to classical methods. Krylov Jacobian methods support large time steps in systems with a few dominant stiff modes but require smaller steps on large systems where stiffness is more evenly distributed. Although theoretical guarantees are limited in nonlinear cases, our results indicate that small nonlinear terms, while leading to large operator-norm errors, have limited practical effect. However, care must be taken in system modelling to avoid discontinuities or singularities in higher derivatives used for the Krylov approximation. When these conditions are met, Krylov Jacobian methods provide a powerful means to enhance the stability of explicit solvers without incurring the cost of full Jacobian evaluation. In many test cases, they permit significantly larger time steps than classical methods like the fourth-order Runge-Kutta (RK4) scheme.