A Randomized Runge-Kutta Method for time-irregular delay differential equations

We deal with approximation of solutions to the following delay differential equations (DDEs)

with a constant time-lag $\tau\in(0,+\infty)$, a fixed time horizon $n\in\mathbb{N}$, a right-hand side function $f:[0,(n+1)\tau]\times\mathbb{R}^{d}\times\mathbb{R}^{d}\mapsto\mathbb{R}^{d}$, and initial-value function: $\varphi(t):[-\tau,0)\mapsto\mathbb{R}^{d}$.

We assume that the function $f=f(t,x,z)$ is integrable with respect to $t$ and (at least) continuous with respect to $(x,z)$.

Building upon the concepts presented in [9], our objective in this paper is to introduce a Randomized Runge-Kutta scheme tailored specifically for time-irregular delay differential equations. This novel scheme differentiates itself from existing methods, including those analyzed in the literature.

The motivation behind exploring these Delay Differential Equations (DDEs) is multifaceted. One aspect stems from the need to model switching systems with memory, as expounded upon in [11] and [12]. Another facet is rooted in practical engineering applications, as evidenced by discussions in [7] and [6]. For instance, in scenarios like the phase change of metallic materials, a DDE becomes crucial due to the time delay in response to changes in processing conditions.

Additionally, inspiration is drawn from delayed differential equations involving rough paths, as represented by the form:

Here, $Z$ symbolizes an integrable perturbation, potentially of stochastic nature, with paths that might exhibit discontinuities. Consequently, if $x(t)=U(t)-Z(t)$, it satisfies the (possibly random) DDE (1.1) with $f(t,x,z)=a(x+Z(t),z+Z(t-\tau))$, and $x(t)=U_{0}$, assuming $Z(t)=0$ for $t\in[-\tau,0]$. In this context, the function $f$ inherits from $Z$ its low smoothness concerning the variable $t$. Noteworthy is the fact that equation (1.2) is a generalization of an ODE with rough paths discussed in [16].

When classical assumptions, such as $C^{r}$-regularity of $f=f(t,x,z)$ concerning all variables $t,x,z$, are imposed on the right-hand side function, errors for deterministic schemes have been documented in [2]. Furthermore, in [7], the error of the Euler scheme has been explored for a certain class of nonlinear DDEs under nonstandard assumptions like the one-side Lipschitz condition and local Hölder continuity. However, limited knowledge exists regarding the approximation of solutions for DDEs with less regular Carathéodory right-hand side functions.

For Carathéodory ODEs, deterministic algorithms encounter convergence issues, necessitating the use of randomized algorithms. In this context, we propose a randomized version of the Runge-Kutta scheme tailored for DDEs of the form (1.1).
In the landscape of existing literature, several relevant contributions are worth noting. [10] provides an analysis of multistep Runge-Kutta methods for delay differential equations, while [18] offers comprehensive insights into numerical methods for delay differential equations. The work by [5] introduces a general-purpose implicit-explicit Runge-Kutta integrator for delay differential equations, and [21] proposes a novel explicit two-stage Runge-Kutta method for delay differential equations with constant delay. Each of these works contributes valuable perspectives to the broader understanding of numerical methods for DDEs.

Despite the extensive exploration of randomized algorithms for ODEs in the literature (see, for instance, [4], [3], [8], [13], [14], [15], [16]), to the best of our knowledge, this paper represents a pioneering effort in defining a randomized Runge-Kutta scheme and rigorously investigating its error for Carathéodory-type DDEs.

The structure of the paper is as follows. In section 2 we give basic notions, definitions, and provide detailed construction and output of the randomized two-stage Runge-Kutta method. Section 3 is devoted to properties of solutions to CDDEs under assumptions stated in Section 2. Section 4 contains detailed error analysis of the randomized R-K method for CDDEs. Finally, details of the Python implementation together with extensive numerical experiments are reported in Section 5.

By $|\cdot|$ we mean the Euclidean norm in $\mathbb{R}^{d}$. We consider a complete probability space $(\Omega,\Sigma,\mathbb{P})$. For a random variable $X:\Omega\to\mathbb{R}$ we denote by $\|X\|_{p}=(\mathbb{E}|X|^{p})^{1/p}$, where $p\in[2,+\infty)$.

Let us fix the horizon parameter $n\in\mathbb{N}$. On the right-hand side function $f$ and initial-value function: $\varphi(t)$, we impose the following assumptions:

1. (A1)

   $f(t,\cdot,\cdot)\in C(\mathbb{R}^{d}\times\mathbb{R}^{d};\mathbb{R}^{d})$ for all $t\in[0,(n+1)\tau]$ and $\varphi\in C([-\tau,0);\mathbb{R}^{d})$ for $t\in[-\tau,0]$.

2. (A2)

   $f(\cdot,x,z)$ is Borel measurable for all $(x,z)\in\mathbb{R}^{d}\times\mathbb{R}^{d}$,

3. (A3)

   There exists $K\in(0,\infty)$ such that for all $t\in[0,(n+1)\tau]$, $x,z\in\mathbb{R}^{d}$

   (2.1)

   $${|f(t,x,z)|\leq K(1+|x|)(1+|z|),}$$

4. (A4)

   There exists $L\in(0,\infty)$ such that for all $t\in[0,(n+1)\tau]$, $x_{1},x_{2},z\in\mathbb{R}^{d}$

   (2.2)

   $${|f(t,x_{1},z)-f(t,x_{2},z)|\leq L(1+|z|)|x_{1}-x_{2}|,}$$

In Section 3, under the assumptions (A1)-(A4), we investigate existence and uniqueness of solution for (1.1). Next, in Section 4 we investigate error of the randomized two-stage Runge-Kutta scheme under slightly stronger assumptions. Namely, we impose Hölder continuity assumptions on $f=(t,x,z)$ with respect to $(t,z)$.

The definition of the two-stage randomized Runge-Kutta method for DDEs goes as follows. $(\gamma_{k}^{j})_{j\in\mathbb{N}_{0},k\in\mathbb{N}}$ iid from $\mathcal{U}(0,1)$, $N\in\mathbb{N}$, $h=\tau/N$, $t_{k}^{j}=j\tau+kh$, $\theta_{k+1}^{j}=t_{k}^{j}+h\gamma_{k+1}^{j}$ for $k=0,1,\ldots,N$, $j=-1,0,1,\ldots,n$. Also let $h_{k+1}^{j}:=h\gamma_{k+1}^{j}$ and define $y_{k}^{-1}=\varphi(t_{k}^{-1})$.

For $j=0$, define

and for $j\geq 1$,

Note that for $j=1$, the delay term $y^{j-2}_{k}$ in the second line of (2.4) is exactly the evaluation of the initial condition $\varphi(t^{-1}_{k})$.

Let us call $\tilde{y}_{k+1}^{j-1,j}$ and $\tilde{y}_{k+1}^{j}$ the intermediate delay term and the intermediate term respectively. A key component here is the intermediate delay term $\tilde{y}_{k+1}^{j-1,j}$: we do not recycle the intermediate term $\tilde{y}_{k+1}^{j-1}$ computed from preceding $\tau$-interval $[(j-1)\tau,j\tau]$ to replace $\tilde{y}_{k+1}^{j-1,j}$ for a simplified computation, because $\tilde{y}_{k+1}^{j-1}$ is generated on $\gamma^{j-1}_{k+1}$, a different random resource from $\gamma^{j}_{k+1}$. It is easy to see that each random vector $y_{k}^{j}$, $j=0,1,\ldots,n$, $k=0,1,\ldots,N$, is measurable with respect to the $\sigma$-field generated by the following family of independent random variables

As the horizon parameter $n$ is fixed, the randomized two-stage Runge-Kutta scheme uses $O(N)$ evaluations of $f$ (with a constant in the $O$ notation that only depends on $n$ but not on $N$).

To help readers better understand the randomized RK method, we present two figures 1(a) and 1(b). Figure 1(a) depicts the scenario when $j=0$, while Figure 1(b) showcases the situation for $j\geq 1$, using $j=1$ as a specific example. Within these figures, different markers are employed to convey distinct meanings: dashed lines represent assignments; black solid lines and solid circles signify computation methods; colored solid circles indicate the values of different terms at the time $t$; the black dots represent the value of time $t$. For clarity, text is color-coded to match the terms they represent. Structurally, the figures are divided into two primary sections vertically along the coordinate axis. The upper section encompasses all computations excluding the generation of $y_{i+1}^{j}$, while the lower section specifically illustrates the operations involved in generating $y_{i+1}^{j}$. In both figures, for any $i$ in $[0,1,\ldots,N-1]$, we aim to compute the value of $y_{i+1}^{j}$ at the time point $t_{i+1}^{j}$. We assume that the values of $y_{i}^{j}$ for all time points before $t_{i+1}^{j}$ are known. Utilizing these known values, we can calculate the value of $y_{i+1}^{j}$ at the time point $t_{i+1}^{j}$ using Eqn. (2.3) if $j=0$ or Eqn. (2.4) if $j\geq 1$.

In this section we investigate the issue of existence and uniqueness of the solution of (1.1) under the assumptions (A1)-(A4).

In the sequel we use the following equivalent representation of the solution of (1.1), that is very convenient when proving its properties and when estimating the error of the randomized RK scheme. For $j=0,1,\ldots,n$ and $t\in[0,\tau]$ it holds

Hence, we take $\phi_{-1}(t):=\varphi(t-\tau)$, $\phi_{j}(t)\mathrel{\mathop{:}}=x(t+j\tau)$ and for $j=0,1,\ldots,n$ we consider the following sequence of initial-value problems

with $g_{j}(t,x)=f(t+j\tau,x,\phi_{j-1}(t))$, $(t,x)\in[0,\tau]\times\mathbb{R}^{d}$. Then the solution of (1.1) can be written as

We prove the following result about existence, uniqueness and Hölder regularity of the solution of the delay differential equation (1.1). We will use this theorem in the next section when proving error estimate for the randomized RK algorithm.

The proof of this theorem follows Theorem 3.1 from [9]. However, due to different initial conditions and assumptions, there are some subtle variations in the proof and conclusions. Specifically, the differences are as follows:

1. (1)

   [9] considers a constant initial condition case, ie, $\varphi(t)=x_{0}$ for $t\in[-\tau,0]$, while in our DDE, a general initial condition $\varphi(t)$ is taken, allowing the initial condition to vary with time.

2. (2)

   Compared to [9, Assumptions A3 and A4], we impose stronger assumptions (A3) and (A4) on $f$ such that the coefficients of linear growth and Lipschitz condition are time-invariant.

In this section we perform detailed error analysis for the randomized RK scheme. As mentioned in Section 1, for the error analysis we impose global Lipschitz assumption on $f=f(t,x,z)$ with respect to $x$ together with global Hölder condition with respect to $z$. Namely, instead of (A3) and (A4), we assume

1. (A3’)

   there exist $\bar{K}\in[0,+\infty)$ such that for all $t\in[0,(n+1)\tau]$,

   (4.1)

   $$|f(t,0,0)|\leq\bar{K},$$

   and there exist $L\in[0,+\infty)$, $\alpha,\gamma\in(0,1]$ such that:
   1. for all $t\in[0,(n+1)\tau]$, $x_{1},x_{2},z_{1},z_{2}\in\mathbb{R}^{d}$,

   (4.2)

   $$|f(t,x_{1},z_{1})-f(t,x_{2},z_{2})|\leq L\Bigl{(}|x_{1}-x_{2}|+|z_{1}-z_{2}|^{\alpha}\Bigr{)},$$

   2. for all $t_{1},t_{2}\in[0,(n+1)\tau]$, $x,z\in\mathbb{R}^{d}$,

   (4.3)

   $$|f(t_{1},x,z)-f(t_{2},x,z)|\leq L(1+|x|+|z|)|t_{1}-t_{2}|^{\gamma}.$$

   3. for all $s,t\in[-\tau,0]$,

   (4.4)

   $$|\varphi(t)-\varphi(s)|\leq L|t-s|,\quad s,t\in[-\tau,0].$$

The main result of this section is as follows.