Explicit Runge-Kutta algorithm to solve non-local equations with memory effects: case of the Maxey-Riley-Gatignol equation

TL;DR

The paper presents an explicit Runge-Kutta framework for solving non-local equations with memory, using the Maxey-Riley-Gatignol (MRG) equation as a representative testbed. By embedding the non-Markovian dynamics into an extended Markovian system through a spectral representation, the authors introduce a history function H(k,t) that evolves locally and enables efficient time stepping. Two concrete RK implementations are provided: a two-stage, first-order method and a four-stage, second-order method, both operating on the inflated state and employing a Clenshaw-Curtis quadrature to compute the history integral. The analysis demonstrates convergence properties and constant memory cost, and numerical experiments in 1D oscillatory forcing and 2D vortex flows validate the expected rates and scalability, while an alternative embedding with Hermite basis is shown to be less efficient. The approach generalizes to other memory-bearing equations possessing a suitable spectral representation and offers restart capability and linear-in-time computational costs, addressing typical memory-growth challenges in non-local-in-time simulations.

Abstract

A standard approach to solve ordinary differential equations, when they describe dynamical systems, is to adopt a Runge-Kutta or related scheme. Such schemes, however, are not applicable to the large class of equations which do not constitute dynamical systems. In several physical systems, we encounter integro-differential equations with memory terms where the time derivative of a state variable at a given time depends on all past states of the system. Secondly, there are equations whose solutions do not have well-defined Taylor series expansion. The Maxey-Riley-Gatignol equation, which describes the dynamics of an inertial particle in nonuniform and unsteady flow, displays both challenges. We use it as a test bed to address the questions we raise, but our method may be applied to all equations of this class. We show that the Maxey-Riley-Gatignol equation can be embedded into an extended Markovian system which is constructed by introducing a new dynamical co-evolving state variable that encodes memory of past states. We develop a Runge-Kutta algorithm for the resultant Markovian system. The form of the kernels involved in deriving the Runge-Kutta scheme necessitates the use of an expansion in powers of $t^{1/2}$. Our approach naturally inherits the benefits of standard time-integrators, namely a constant memory storage cost, a linear growth of operational effort with simulation time, and the ability to restart a simulation with the final state as the new initial condition.

Figures (5)

• Figure 1: A schematic for the Markovian embedding procedure for the MRG equation (Section \ref{['sec3']}). (a) Non-local inter-state interaction in the original representation Eq. \ref{['Eq:compactMRG']}. Interaction between no two states is the same. (b) Local interaction via Eq. \ref{['newstatetransformed']} due to inflated state description with the introduction of new 'history function', $H(k,t)$. All local interactions are identical.
• Figure 2: Numerical solution of 1D scalar MRG equation forced with $N = \sin(\omega t)$ and subject to initial condition, $w_0=1$, for the parameters $(\alpha,\gamma,\omega,T) = (0.33,1,5,5)$ using the 2- and 4-stage schemes (\ref{['oc:2f']} and \ref{['oc:4f']}). $M=51$ (Chebyshev) quadrature points were used to compute the history integral. (a) Comparison of the exact solution with the numerical solutions for $\Delta t = 2^{-3}$. (b) Error as a function of time for $\Delta t = 2^{-3}$. (c) Scaling of error with time-step $\Delta t$. The errors are measured against the analytical solution (\ref{['exactsinsol']}) in the $l_2-$norm. The slopes are consistent with the estimated orders of accuracy of schemes. (d) Scaling of operational cost with error. CPU time is used as a proxy for the operational cost. The slopes of the curves indicate the scheme's cost to improve in accuracy. As per the scheme's algorithm, $l_2$ error $\sim \Delta t^p$, where $p=\{1,2\}$ is the order of accuracy, and operational cost $\sim \Delta t^{-1}$. The slopes of the curves verify this estimate. The codes for the schemes used in this example are available on GitHub at https://github.com/jagannathan-divya/rk4mrg.
• Figure 3: Numerical advection of a particle in 2D stationary Lamb-Oseen vortex starting at $\textbf{y}_0=(1,0)$, with non-zero initial slip velocity $\textbf{w}_0=(1,0)$, for the parameters $(\alpha,\gamma) = (1,1)$ using our 2-stage and 4-stage schemes. $M=51$ (Chebyshev) quadrature points were used to compute the history integral. (a) Trajectory of the particle evolved up to $T=200$ with $\Delta t=2^{-3}$. (b) Scaling of error with time-step $\Delta t$ for simulations run up to $T=5$. Here, the errors are defined for the magnitude of slip velocity vector, $|\textbf{w}|$, and measured in the $l_2-$norm against a fine-resolution numerical solution computed using our 4-stage scheme with $\Delta t = 2^{-16}$ and $M=101$. The codes for the schemes used in this example are available on GitHub at https://github.com/jagannathan-divya/rk4mrg.
• Figure 4: Repeat of numerical experiment \ref{['sec6p1']} under the alternative Markovian embedding procedure (Section \ref{['sec6p3']}) for parameters $(\alpha,\gamma,\omega,T)=(0.33,1,5,5)$ using ETD2RK (Eq. 81, 82 in CM2002). $M$ is the number of basis Hermite functions used to compute the history integral. Compare the red curve here with the red curve in Fig \ref{['fig:sinusoid:err_v_dt']} both of which correspond to $M=51$. The codes used for the example are available on GitHub at https://github.com/jagannathan-divya/rk4mrg.
• Figure 5: Contour $\partial D^-$ (continuous red).