Local linearization method

In numerical analysis, the local linearization (LL) method is a general strategy for designing numerical integrators for differential equations based on a local (piecewise) linearization of the given equation on consecutive time intervals. The numerical integrators are then iteratively defined as the solution of the resulting piecewise linear equation at the end of each consecutive interval. The LL method has been developed for a variety of equations such as the ordinary, delayed, random and stochastic differential equations. The LL integrators are key component in the implementation of inference methods for the estimation of unknown parameters and unobserved variables of differential equations given time series of (potentially noisy) observations. The LL schemes are ideals to deal with complex models in a variety of fields as neuroscience, finance, forestry management, control engineering, mathematical statistics, etc.

Background

Differential equations have become an important mathematical tool for describing the time evolution of several phenomenon, e.g., rotation of the planets around the sun, the dynamic of assets prices in the market, the fire of neurons, the propagation of epidemics, etc. However, since the exact solutions of these equations are usually unknown, numerical approximations to them obtained by numerical integrators are necessary. Currently, many applications in engineering and applied sciences focused in dynamical studies demand the developing of efficient numerical integrators that preserve, as much as possible, the dynamics of these equations. With this main motivation, the Local Linearization integrators have been developed.

High-order local linearization method

High-order local linearization (HOLL) method is a generalization of the Local Linearization method oriented to obtain high-order integrators for differential equations that preserve the stability and dynamics of the linear equations. The integrators are obtained by splitting, on consecutive time intervals, the solution x of the original equation in two parts: the solution z of the locally linearized equation plus a high-order approximation of the residual $\mathbf {r} =\mathbf {x} -\mathbf {z}$ .

Local linearization scheme

A Local Linearization (LL) scheme is the final recursive algorithm that allows the numerical implementation of a discretization derived from the LL or HOLL method for a class of differential equations.

LL methods for ODEs

Consider the d-dimensional Ordinary Differential Equation (ODE)

${\frac {d\mathbf {x} \left(t\right)}{dt}}=\mathbf {f} \left(t,\mathbf {x} \left(t\right)\right),\qquad t\in \left[t_{0},T\right],\qquad \qquad \qquad \qquad (4.1)$

with initial condition $\mathbf {x} (t_{0})=\mathbf {x} _{0}$ , where $\mathbf {f}$ is a differentiable function.

Let $\left(t\right)_{h}=\{t_{n}:n=0,..,N\}$ be a time discretization of the time interval $[t_{0},T]$ with maximum stepsize h such that $t_{n}<t_{n+1}$ and $h_{n}=t_{n+1}-t_{n}\leq h$ . After the local linearization of the equation (4.1) at the time step $t_{n}$ the variation of constants formula yields

$\mathbf {x} (t_{n}+h)=\mathbf {x} (t_{n})+\mathbf {\phi } (t_{n},\mathbf {x} (t_{n});h)+\mathbf {r} (t_{n},\mathbf {x} (t_{n});h),$

where

$\mathbf {\phi } (t_{n},\mathbf {z} _{n};h)=\int \limits _{0}^{h}e^{\mathbf {f} _{\mathbf {x} }(t_{n},\mathbf {z} _{n})(h-s)}(\mathbf {f} (t_{n},\mathbf {z} _{n})+\mathbf {f} _{t}(t_{n},\mathbf {z} _{n})s)\,ds\qquad$

results from the linear approximation, and

$\mathbf {r} (t_{n},\mathbf {z} _{n};h)=\int \limits _{0}^{h}e^{\mathbf {f} _{\mathbf {x} }(t_{n},\mathbf {z} _{n})(h-s)}\mathbf {g} _{n}(s,\mathbf {x} (t_{n}+s))\,ds,\qquad \qquad \qquad (4.2)$

is the residual of the linear approximation. Here, $\mathbf {f} _{\mathbf {x} }$ and $\mathbf {f} _{t}$ denote the partial derivatives of f with respect to the variables x and t, respectively, and $\mathbf {g} _{n}(s,\mathbf {u} )=\mathbf {f} (s,\mathbf {u} )-\mathbf {f} _{\mathbf {x} }(t_{n},\mathbf {z} _{n})\mathbf {u} -\mathbf {f} _{t}(t_{n},\mathbf {z} _{n})(s-t_{n})-\mathbf {f} (t_{n},\mathbf {z} _{n})+\mathbf {f} _{\mathbf {x} }(t_{n},\mathbf {z} _{n})\mathbf {z} _{n}.$

Local linear discretization

For a time discretization $\left(t\right)_{h}$ , the Local Linear discretization of the ODE (4.1) at each point $t_{n+1}\in \left(t\right)_{h}$ is defined by the recursive expression ^[1]^[2]

$\mathbf {z} _{n+1}=\mathbf {z} _{n}+\mathbf {\phi } (t_{n},\mathbf {z} _{n};h_{n}),\qquad {\text{ with }}\quad \mathbf {z} _{0}=\mathbf {x} _{0}.\qquad \qquad \qquad \qquad (4.3)$

The Local Linear discretization (4.3) converges with order 2 to the solution of nonlinear ODEs, but it match the solution of the linear ODEs. The recursion (4.3) is also known as Exponential Euler discretization.^[3]

High-order local linear discretizations

For a time discretization $(t)_{h},$ a high-order local linear (HOLL) discretization of the ODE (4.1) at each point $t_{n+1}\in (t)_{h}$ is defined by the recursive expression ^[1]^[4]^[5]^[6]

$\mathbf {z} _{n+1}=\mathbf {z} _{n}+\mathbf {\phi } (t_{n},\mathbf {z} _{n};h_{n})+{\widetilde {\mathbf {r} }}(t_{n},\mathbf {z} _{n};h_{n}),\qquad {\text{ with }}\quad \mathbf {z} _{0}=\mathbf {x} _{0},\qquad \qquad \qquad (4.4)$

where ${\tilde {r}}$ is an order $\alpha$ (> 2) approximation to the residual r $(i.e.,\left\vert \mathbf {r} (t_{n},\mathbf {z} _{n};h)-{\widetilde {\mathbf {r} }}(t_{n},\mathbf {z} _{n};h)\right\vert \propto h^{\alpha +1}).$ The HOLL discretization (4.4) converges with order $\alpha$ to the solution of nonlinear ODEs, but it match the solution of the linear ODEs.

HOLL discretizations can be derived in two ways:^[1]^[4]^[5]^[6] 1) (quadrature-based) by approximating the integral representation (4.2) of r; and 2) (integrator-based) by using a numerical integrator for the differential representation of r defined by

${\frac {d\mathbf {r} (t)}{dt}}=\mathbf {q} (t_{n},\mathbf {z} _{n};t\mathbf {,\mathbf {r} } (t)\mathbf {)} ,\qquad {\text{ with }}\qquad \mathbf {r} (t_{n})=\mathbf {0,} \qquad \qquad \qquad (4.5)$

for all $t\in \lbrack t_{k},t_{k+1}]$ , where

$\mathbf {q} (t_{n},\mathbf {z} _{n};s\mathbf {,\xi } )=\mathbf {f} (s,\mathbf {z} _{n}+\mathbf {\phi } \left(t_{n},\mathbf {z} _{n};s-t_{n}\right)+\mathbf {\xi } )-\mathbf {f} _{\mathbf {x} }(t_{n},\mathbf {z} _{n})\mathbf {\phi } (t_{n},\mathbf {z} _{n};s-t_{n})-\mathbf {f} _{t}(t_{n},\mathbf {z} _{n})(s-t_{n})-\mathbf {f} (t_{n},\mathbf {z} _{n}).$

HOLL discretizations are, for instance, the followings:

Locally Linearized Runge Kutta discretization^[6]^[4]

$\qquad \mathbf {z} _{n+1}=\mathbf {z} _{n}+\mathbf {\phi } (t_{n},\mathbf {z} _{n};h_{n})+h_{n}\sum _{j=1}^{s}b_{j}\mathbf {k} _{j},\quad {\text{ with }}\quad \mathbf {k} _{i}=\mathbf {q} (t_{n},\mathbf {z} _{n};{\text{ }}t_{n}+c_{i}h_{n}\mathbf {,} \mathbf {} h_{n}\sum _{j=1}^{i-1}a_{ij}\mathbf {k} _{j}),$

which is obtained by solving (4.5) via a s-stage explicit Runge–Kutta (RK) scheme with coefficients $\mathbf {c} =\left[c_{i}\right],\mathbf {A} =\left[a_{ij}\right]\quad and\quad \mathbf {b} =\left[b_{j}\right]$ .

Local linear Taylor discretization^[5]

$\mathbf {z} _{n+1}=\mathbf {z} _{n}+\mathbf {\phi } (t_{n},\mathbf {z} _{n};h_{n})+\int _{0}^{h_{n}}e^{(h_{n}-s)\mathbf {f} _{\mathbf {x} }\left(t_{n},\mathbf {z} _{n}\right)}\sum _{j=2}^{p}{\frac {\mathbf {c} _{n,j}}{j!}}s^{j}\,ds,{\text{ with }}\mathbf {c} _{n,j}=\left({\frac {d^{j+1}\mathbf {x} (t)}{dt^{j+1}}}-\mathbf {f} _{\mathbf {x} }(t_{n},\mathbf {z} _{n}){\frac {d^{j}\mathbf {x} (t)}{dt^{j}}}\right)\mid _{t=\mathbf {z} _{n}},$

which results from the approximation of $\mathbf {g} _{n}$ in (4.2) by its order-p truncated Taylor expansion.

Multistep-type exponential propagation discretization

${\textstyle \mathbf {z} _{n+1}=\mathbf {z} _{n}+\mathbf {\phi } (t_{n},\mathbf {z} _{n};h)+h\sum _{j=0}^{p-1}\gamma _{j}\nabla ^{j}\mathbf {g} _{n}(t_{n},\mathbf {z} _{n}),\quad with\quad \gamma _{j}=(-1)^{j}\int \limits _{0}^{1}e^{(1-\theta )h\mathbf {f} _{\mathbf {x} }(t_{n},\mathbf {z} _{n})}\left({\begin{array}{c}-\theta \\j\end{array}}\right)d\theta ,}$

which results from the interpolation of $\mathbf {g} _{n}$ in (4.2) by a polynomial of degree p on $t_{n},\ldots ,t_{n-p+1}$ , where $\nabla ^{j}\mathbf {g} _{n}(t_{m},\mathbf {z} _{m})$ denotes the j-th backward difference of $\mathbf {g} _{n}(t_{m},\mathbf {z} _{m})$ .

Runge Kutta type Exponential Propagation discretization ^[7]

${\textstyle \mathbf {z} _{n+1}=\mathbf {z} _{n}+\mathbf {\phi } (t_{n},\mathbf {z} _{n};h)+h\sum _{j=0}^{p-1}\gamma _{j,p}\nabla ^{j}\mathbf {g} _{n}(t_{n},\mathbf {z} _{n}),\quad {\text{ with }}\quad \gamma _{j,p}=\int \limits _{0}^{1}e^{(1-\theta )h\mathbf {f} _{\mathbf {x} }(t_{n},\mathbf {z} _{n})}\left({\begin{array}{c}\theta p\\j\end{array}}\right)d\theta ,}$

which results from the interpolation of $\mathbf {g} _{n}$ in (4.2) by a polynomial of degree p on $t_{n},\ldots ,t_{n}+(p-1)h/p$ ,

Linealized exponential Adams discretization^[8]

${\textstyle \mathbf {z} _{n+1}=\mathbf {z} _{n}+\mathbf {\phi } (t_{n},\mathbf {z} _{n};h)+h\sum _{j=1}^{p-1}\sum _{l=1}^{j}{\frac {\gamma _{j+1}}{l}}\nabla ^{l}\mathbf {g} _{n}(t_{n},\mathbf {z} _{n}),\quad {\text{ with }}\quad \gamma _{j+1}=(-1)^{j+1}\int \limits _{0}^{1}e^{(1-\theta )h\mathbf {f} _{\mathbf {x} }\left(t_{n},\mathbf {z} _{n}\right)}\theta \left({\begin{array}{c}-\theta \\j\end{array}}\right)d\theta ,}$

which results from the interpolation of $\mathbf {g} _{n}$ in (4.2) by a Hermite polynomial of degree p on $t_{n},\ldots ,t_{n-p+1}$ .

Local linearization schemes

All numerical implementation $\mathbf {y} _{n}$ of the LL (or of a HOLL) discretization $\mathbf {z} _{n}$ involves approximations ${\widetilde {\phi }}_{j}$ to integrals $\phi _{j}$ of the form

$\phi _{j}(\mathbf {A} ,h)=\int \limits _{0}^{h}e^{(h-s)\mathbf {A} }s^{j-1}\,ds,\qquad j=1,2\ldots ,$

where A is a d × d matrix. Every numerical implementation $\mathbf {y} _{n}$ of the LL (or of a HOLL) $\mathbf {z} _{n}$ of any order is generically called Local Linearization scheme.^[1]^[9]

Computing integrals involving matrix exponential

Among a number of algorithms to compute the integrals $\phi _{j}$ , those based on rational Padé and Krylov subspaces approximations for exponential matrix are preferred. For this, a central role is playing by the expression^[10]^[5]^[11]

$\sum \nolimits _{i=1}^{l}\phi _{i}(\mathbf {A} ,h)\mathbf {a} _{i}=\mathbf {L} e^{h\mathbf {H} }\mathbf {r} ,$

where $\mathbf {a} _{i}$ are d-dimensional vectors,

$\mathbf {H} ={\begin{bmatrix}\mathbf {A} &\mathbf {v} _{l}&\mathbf {v} _{l-1}&\cdots &\mathbf {v} _{1}\\\mathbf {0} &\mathbf {0} &1&\cdots &0\\\mathbf {0} &\mathbf {0} &0&\ddots &0\\\vdots &\vdots &\vdots &\ddots &1\\\mathbf {0} &\mathbf {0} &0&\cdots &0\end{bmatrix}}\in \mathbb {R} ^{(d+l)\times (d+l)},$

$\mathbf {L} =[\mathbf {I} \quad \mathbf {0} _{d\times l}]$ , $\mathbf {r} =[\mathbf {0} _{1\times (d+l-1)}\quad 1]^{\intercal },$ $\mathbf {v} _{i}=\mathbf {a} _{i}(i-1)!$ , being $\mathbf {I}$ the d-dimensional identity matrix.

If $\mathbf {P} _{p,q}(2^{-k}\mathbf {H} h)$ denotes the (p; q)-Padé approximation of $e^{2^{-k}\mathbf {H} h}$ and k is the smallest natural number such that $|2^{-k}\mathbf {H} h|\leq {\frac {1}{2}},then$ ^[12]^[9]

$\left\vert \sum \nolimits _{i=1}^{l}\phi _{i}(\mathbf {A} ,h)\mathbf {a} _{i}-\mathbf {L} \left(\mathbf {\mathbf {P} } _{p,q}(2^{-k}\mathbf {H} h)\right)^{2^{k}}\mathbf {r} \right\vert \varpropto h^{p+q+1}.$

If $\mathbf {\mathbf {k} } _{m,k}^{p,q}(h,\mathbf {H} ,\mathbf {r} )$ denotes the (m; p; q; k) Krylov-Padé approximation of $e^{h\mathbf {H} }\mathbf {r}$ , then ^[12]

$\left\vert \sum \nolimits _{i=1}^{l}\phi _{i}(\mathbf {A} ,h)\mathbf {a} _{i}-\mathbf {L\mathbf {k} } _{m,k}^{p,q}(h,\mathbf {H} ,\mathbf {r} )\right\vert \varpropto h^{\min({m,p+q+1})},$

where $m\leq d$ is the dimension of the Krylov subspace.

Order-2 LL schemes

$\mathbf {y} _{n+1}=\mathbf {y} _{n}+\mathbf {L} (\mathbf {P} _{p,q}(2^{-k_{n}}\mathbf {M} _{n}h_{n}))^{2^{k_{n}}}\mathbf {r,}$ ^[13]^[9] $\qquad \qquad (4.6)$

where the matrices $\mathbf {M} _{n}$ , L and r are defined as

$\mathbf {M} _{n}={\begin{bmatrix}\mathbf {f} _{\mathbf {x} }(t_{n},\mathbf {y} _{n})&\mathbf {f} _{t}(t_{n},\mathbf {y} _{n})&\mathbf {f} (t_{n},\mathbf {y} _{n})\\0&0&1\\0&0&0\end{bmatrix}}\in \mathbb {R} ^{(d+2)\times (d+2)},$

$\mathbf {L} =\left[{\begin{array}{ll}\mathbf {I} &\mathbf {0} _{d\times 2}\end{array}}\right]$ and $\mathbf {r} ^{\intercal }=\left[{\begin{array}{ll}\mathbf {0} _{1\times (d+1)}&1\end{array}}\right]$ with $p+q>1$ . For large systems of ODEs ^[3]

$\mathbf {y} _{n+1}=\mathbf {y} _{n}+\mathbf {L\mathbf {k} } _{m_{n},k_{n}}^{p,q}(h_{n},\mathbf {M} _{n},\mathbf {r} )\mathbf {,} \qquad {\text{ with }}\qquad m_{n}>2.$

Order-3 LL-Taylor schemes

$\mathbf {y} _{n+1}=\mathbf {y} _{n}+\mathbf {L} _{1}(\mathbf {P} _{p,q}(2^{-k_{n}}\mathbf {T} _{n}h_{n}))^{2^{k_{n}}}\mathbf {r} _{1}\mathbf {,}$ ^[5] $\qquad \qquad (4.7)$

where for autonomous ODEs the matrices $\mathbf {T} _{n},\mathbf {L} _{1}$ and $\mathbf {r} _{1}$ are defined as

$\mathbf {T} _{n}=\left[{\begin{array}{cccc}\mathbf {f} _{\mathbf {x} }(\mathbf {y} _{n})&(\mathbf {I} \otimes \mathbf {f} ^{\intercal }(\mathbf {y} _{n}))\mathbf {f} _{\mathbf {xx} }(\mathbf {y} _{n})\mathbf {f} (\mathbf {y} _{n})&\mathbf {0} &\mathbf {f} (\mathbf {y} _{n})\\0&0&0&0\\0&0&0&1\\0&0&0&0\end{array}}\right]\in \mathbb {R} ^{(d+3)\times (d+3)},$

$\mathbf {L} _{1}=\left[{\begin{array}{ll}\mathbf {I} &\mathbf {0} _{d\times 3}\end{array}}\right]\quad and\quad \mathbf {r} _{1}^{\intercal }=\left[{\begin{array}{ll}\mathbf {0} _{1\times (d+2)}&1\end{array}}\right]$ . Here, $\mathbf {f} _{\mathbf {xx} }$ denotes the second derivative of f with respect to x, and p + q > 2. For large systems of ODEs

$\mathbf {y} _{n+1}=\mathbf {y} _{n}+\mathbf {L\mathbf {k} } _{m_{n},k_{n}}^{p,q}(h_{n},\mathbf {T} _{n},\mathbf {r} )\mathbf {,} \qquad {\text{ with }}\qquad m_{n}>3.$

Order-4 LL-RK schemes

$\mathbf {y} _{n+1}=\mathbf {y} _{n}+\mathbf {u} _{4}+{\frac {h_{n}}{6}}(2\mathbf {k} _{2}+2\mathbf {k} _{3}+\mathbf {k} _{4}),\quad$ ^[4]^[6] $\qquad \qquad (4.8)$

where

$\mathbf {u} _{j}=\mathbf {L} (\mathbf {P} _{p,q}(2^{-\kappa _{j}}\mathbf {M} _{n}c_{j}h_{n}))^{2^{\kappa _{j}}}\mathbf {r}$

and

$\mathbf {k} _{j}=\mathbf {f} \left(t_{n}+c_{j}h_{n},\mathbf {y} _{n}+\mathbf {u} _{j}+c_{j}h_{n}\mathbf {k} _{j-1}\right)-\mathbf {f} \left(t_{n},\mathbf {y} _{n}\right)-\mathbf {f} _{\mathbf {x} }\left(t_{n},\mathbf {y} _{n}\right)\mathbf {u} _{j}\ -\mathbf {f} _{t}\left(t_{n},\mathbf {y} _{n}\right)c_{j}h_{n},$

with $\mathbf {k} _{1}\equiv \mathbf {0} ,c=\left[{\begin{array}{cccc}0&{\frac {1}{2}}&{\frac {1}{2}}&1\end{array}}\right],$ and p + q > 3. For large systems of ODEs, the vector $\mathbf {u} _{j}$ in the above scheme is replaced by $\mathbf {u} _{j}=\mathbf {L\mathbf {k} } _{m_{j},k_{j}}^{p,q}(c_{j}h_{n},\mathbf {M} _{n},\mathbf {r} )$ with $m_{j}>4.$

Locally linearized Runge–Kutta scheme of Dormand and Prince

$\mathbf {y} _{n+1}=\mathbf {y} _{n}+\mathbf {u} _{s}+h_{n}\sum _{j=1}^{s}b_{j}\mathbf {k} _{j}\qquad {\text{ and }}\qquad {\widehat {\mathbf {y} }}_{n+1}=\mathbf {y} _{n}+\mathbf {u} _{s}+h_{n}\sum _{j=1}^{s}{\widehat {b}}_{j}\mathbf {k} _{j},\quad$ ^[14]^[15] $\qquad \qquad (4.9)$

where s = 7 is the number of stages,

$\mathbf {k} _{j}=\mathbf {f(} t_{n}+c_{j}h_{n},\mathbf {y} _{n}+\mathbf {u} _{j}+h_{n}\sum _{i=1}^{s-1}a_{j,i}\mathbf {k} _{i})-\mathbf {f} \left(t_{n},\mathbf {y} _{n}\right)-\mathbf {f} _{\mathbf {x} }\left(t_{n},\mathbf {y} _{n}\right)\mathbf {u} _{j}\ -\mathbf {f} _{t}\left(t_{n},\mathbf {y} _{n}\right)c_{j}h_{n},$

with $\mathbf {k} _{1}\equiv \mathbf {0}$ , and $a_{j,i},b_{j},{\widehat {b}}_{j}\quad and\quad c_{j}$ are the Runge–Kutta coefficients of Dormand and Prince and p + q > 4. The vector $\mathbf {u} _{j}$ in the above scheme is computed by a Padé or Krylor–Padé approximation for small or large systems of ODE, respectively.

Stability and dynamics

By construction, the LL and HOLL discretizations inherit the stability and dynamics of the linear ODEs, but it is not the case of the LL schemes in general. With $p\leq q\leq p+2$ , the LL schemes (4.6)-(4.9) are A-stable.^[4] With q = p + 1 or q = p + 2, the LL schemes (4.6)–(4.9) are also L-stable.^[4] For linear ODEs, the LL schemes (4.6)-(4.9) converge with order p + q.^[4]^[9] In addition, with p = q = 6 and $m_{n}$ = d, all the above described LL schemes yield to the ″exact computation″ (up to the precision of the floating-point arithmetic) of linear ODEs on the current personal computers.^[4]^[9] This includes stiff and highly oscillatory linear equations. Moreover, the LL schemes (4.6)-(4.9) are regular for linear ODEs and inherit the symplectic structure of Hamiltonian harmonic oscillators.^[5]^[13] These LL schemes are also linearization preserving, and display a better reproduction of the stable and unstable manifolds around hyperbolic equilibrium points and periodic orbits that other numerical schemes with the same stepsize.^[5]^[13] For instance, Figure 1 shows the phase portrait of the ODEs

{\begin{aligned}&{\frac {dx_{1}}{dt}}=-2x_{1}+x_{2}+1-\mu f(x_{1},\lambda )\qquad \qquad (4.10)\\[6pt]&{\frac {dx_{2}}{dt}}=x_{1}-2x_{2}+1-\mu f(x_{2},\lambda )\qquad \qquad \quad (4.11)\end{aligned}}

with $f(u,\lambda )=u(1+u+\lambda u^{2})^{-1}$ , $\mu =15$ and $\lambda =57$ , and its approximation by various schemes. This system has two stable stationary points and one unstable stationary point in the region $0\leq x_{1},x_{2}\leq 1$ .

LL methods for DDEs

Consider the d-dimensional Delay Differential Equation (DDE)

${\frac {d\mathbf {x} (t)}{dt}}=\mathbf {f} (t,\mathbf {x} (t),\mathbf {x} _{t}(-\tau _{1}),\ldots ,\mathbf {x} _{t}(-\tau _{m})),\qquad t\in [t_{0},T],\qquad \qquad (5.1)$

with m constant delays $\tau _{i}>0$ and initial condition $\mathbf {x} _{t_{0}}(s)=\mathbf {\varphi } (s)$ for all $s\in [-\tau ,0],$ where f is a differentiable function, $\mathbf {x} _{t}:[-\tau ,0]\longrightarrow \mathbb {R} ^{d}$ is the segment function defined as

$\mathbf {x} _{t}(s):=\mathbf {x} (t+s),{\text{ }}s\in [-\tau ,0],$

for all $t\in [t_{0},T],\mathbf {\varphi } :[-\tau ,0]\longrightarrow \mathbb {R} ^{d}$ is a given function, and $\tau =\max \left\{\tau _{1},\ldots ,\tau _{m}\right\}.$

Local linear discretization

For a time discretization $(t)_{h}$ , the Local Linear discretization of the DDE (5.1) at each point $t_{n+1}\in (t)_{h}$ is defined by the recursive expression ^[11]

$\mathbf {z} _{n+1}=\mathbf {z} _{n}+\Phi (t_{n},\mathbf {z} _{n},h_{n};{\widetilde {\mathbf {z} }}_{t_{n}}^{1},\ldots ,{\widetilde {\mathbf {z} }}_{t_{n}}^{m}),\qquad \qquad (5.2)$

where

$\Phi (t_{n},\mathbf {z} _{n},h_{n};{\widetilde {\mathbf {z} }}_{t_{n}}^{1},\ldots ,{\widetilde {\mathbf {z} }}_{t_{n}}^{m})=\int \limits _{0}^{h_{n}}e^{\mathbf {A} _{n}(h_{n}-u)}\left[\sum \limits _{i=1}^{m}\mathbf {B} _{n}^{i}({\widetilde {\mathbf {z} }}_{t_{n}}^{i}(u-\tau _{i})-{\widetilde {\mathbf {z} }}_{t_{n}}^{i}(-\tau _{i}))+\mathbf {d} _{n}\right]\,du+\int \limits _{0}^{h_{n}}\int \limits _{0}^{u}e^{\mathbf {A} _{n}(h_{n}-u)}\mathbf {c} _{n}\,dr\,du$

${\widetilde {\mathbf {z} }}_{t_{n}}^{i}:\left[-\tau _{i},0\right]\longrightarrow \mathbb {R} ^{d}$ is the segment function defined as

${\widetilde {\mathbf {z} }}_{t_{n}}^{i}(s):={\widetilde {\mathbf {z} }}^{i}(t_{n}+s),{\text{ }}s\in [-\tau _{i},0],$

and ${\widetilde {\mathbf {z} }}^{i}:\left[t_{n}-\tau _{i},t_{n}\right]\longrightarrow \mathbb {R} ^{d}$ is a suitable approximation to $\mathbf {x} (t)$ for all $t\in \lbrack t_{n}-\tau _{i},t_{n}]$ such that ${\widetilde {\mathbf {z} }}^{i}(t_{n})=\mathbf {z} _{n}.$ Here,

$\mathbf {A} _{n}=\mathbf {f} _{x}(t_{n},\mathbf {z} _{n},{\widetilde {\mathbf {z} }}_{t_{n}}^{1}(-\tau _{1}),\ldots ,{\widetilde {\mathbf {z} }}_{t_{n}}^{m}(-\tau _{d})),{\text{ }}\mathbf {B} _{n}^{i}=\mathbf {f} _{x_{t}(-\tau _{i})}(t_{n},\mathbf {z} _{n},{\widetilde {\mathbf {z} }}_{t_{n}}^{1}(-\tau _{1}),\ldots ,{\widetilde {\mathbf {z} }}_{t_{n}}^{m}(-\tau _{d}))$

are constant matrices and

$\mathbf {c} _{n}=\mathbf {f} _{t}(t_{n},\mathbf {z} _{n},{\widetilde {\mathbf {z} }}_{t_{n}}^{1}(-\tau _{1}),\ldots ,{\widetilde {\mathbf {z} }}_{t_{n}}^{m}(-\tau _{d})){\text{ and }}\mathbf {d} _{n}=\mathbf {f(} t_{n},\mathbf {z} _{n},{\widetilde {\mathbf {z} }}_{t_{n}}^{1}(-\tau _{1}),\ldots ,{\widetilde {\mathbf {z} }}_{t_{n}}^{m}(-\tau _{d}))$

are constant vectors. $\mathbf {f} _{t},\mathbf {f} _{x}\quad and\quad \mathbf {f} _{x_{t}(-\tau _{i})}$ denote, respectively, the partial derivatives of f with respect to the variables t and x, and $\mathbf {x} _{t}(-\tau _{i})$ . The Local Linear discretization (5.2) converges to the solution of (5.1) with order $\alpha =\min\{2,r\},$ if ${\widetilde {\mathbf {z} }}_{t_{n}}^{i}$ approximates $\mathbf {z} _{t_{n}}^{i}$ with order $r\quad (i.e.,\left\vert \mathbf {z} _{t_{n}}^{i}\mathbf {(} u-\tau _{i}\mathbf {)} -{\widetilde {\mathbf {z} }}_{t_{n}}^{i}\mathbf {(} u-\tau _{i}\mathbf {)} \right\vert \propto h_{n}^{r}$ for all $u\in \lbrack 0,h_{n}])$ .

Local linearization schemes

Depending on the approximations ${\widetilde {\mathbf {z} }}_{t_{n}}^{i}$ and on the algorithm to compute $\mathbf {\phi }$ different Local Linearizations schemes can be defined. Every numerical implementation $\mathbf {y} _{n}$ of a Local Linear discretization $\mathbf {z} _{n}$ is generically called local linearization scheme.

Order-2 polynomial LL schemes

$\mathbf {y} _{n+1}=\mathbf {y} _{n}+\mathbf {L} (\mathbf {P} _{p,q}(2^{-k_{n}}\mathbf {M} _{n}h_{n}))^{2^{k_{n}}}\mathbf {r,} \quad$ ^[11] $\qquad (5.3)$

where the matrices $\mathbf {M} _{n},\mathbf {L}$ and $\mathbf {r}$ are defined as

$\mathbf {M} _{n}={\begin{bmatrix}\mathbf {A} _{n}&\mathbf {c} _{n}+\sum \limits _{i=1}^{m}\mathbf {B} _{n}^{i}\mathbf {\alpha } _{n}^{i}&\mathbf {d} _{n}\\0&0&1\\0&0&0\end{bmatrix}}\in \mathbb {R} ^{(d+2)\times (d+2)},$

$\mathbf {L} =\left[{\begin{array}{ll}\mathbf {I} &\mathbf {0} _{d\times 2}\end{array}}\right]$ and $\mathbf {r} ^{\intercal }=\left[{\begin{array}{ll}\mathbf {0} _{1\times (d+1)}&1\end{array}}\right],h_{n}\leq \tau$ , and $p+q>1$ . Here, the matrices $\mathbf {A} _{n}$ , $\mathbf {B} _{n}^{i}$ , $\mathbf {c} _{n}$ and $\mathbf {d} _{n}$ are defined as in (5.2), but replacing $\mathbf {z}$ by $\mathbf {y}$ and $\mathbf {\alpha } _{n}^{i}=(\mathbf {y} (t_{n+1}-\tau _{i})-\mathbf {y} (t_{n}-\tau _{i}))/h_{n},$ where

$\mathbf {y} \left(t\right)=\mathbf {y} _{n_{t}}+\mathbf {L} (\mathbf {P} _{p,q}(2^{-k_{n}}\mathbf {M} _{n_{t}}(t-t_{n_{t}})))^{2^{k_{n}}}\mathbf {r} ,$

with $n_{t}=\max\{n=0,1,2,...,:t_{n}\leq t{\text{ and }}t_{n}\in \left(t\right)_{h}\}$ , is the Local Linear Approximation to the solution of (5.1) defined through the LL scheme (5.3) for all $t\in \lbrack t_{0},t_{n}]$ and by $\mathbf {y} \left(t\right)=\mathbf {\varphi } \left(t\right)$ for $t\in \left[t_{0}-\tau ,t_{0}\right]$ . For large systems of DDEs

$\mathbf {y} _{n+1}=\mathbf {y} _{n}+\mathbf {L\mathbf {k} } _{m_{n},k_{n}}^{p,q}(h_{n},\mathbf {M} _{n},\mathbf {r} )\quad and\quad \mathbf {y} \left(t\right)=\mathbf {y} _{n_{t}}+\mathbf {L\mathbf {k} } _{m_{n_{t}},k_{n_{t}}}^{p,q}(t-t_{n_{t}},\mathbf {M} _{n_{t}},\mathbf {r} ),$

with $p+q>1$ and $m_{n}>2$ . Fig. 2 Illustrates the stability of the LL scheme (5.3) and of that of an explicit scheme of similar order in the integration of a stiff system of DDEs.

LL methods for RDEs

Consider the d-dimensional Random Differential Equation (RDE)

{\frac {d\mathbf {x} \left(t\right)}{dt}}=\mathbf {f} (\mathbf {x} (t),\mathbf {\xi } (t)),\quad t\in \left[t_{0},T\right],\qquad \qquad \qquad (6.1)

with initial condition $\mathbf {x} (t_{0})=\mathbf {x} _{0},$ where $\mathbf {\xi }$ is a k-dimensional separable finite continuous stochastic process, and f is a differentiable function. Suppose that a realization (path) of $\mathbf {\xi }$ is given.

Local Linear discretization

For a time discretization $\left(t\right)_{h}$ , the Local Linear discretization of the RDE (6.1) at each point $t_{n+1}\in \left(t\right)_{h}$ is defined by the recursive expression ^[16]

$\mathbf {z} _{n+1}=\mathbf {z} _{n}+\mathbf {\phi } (t_{n},\mathbf {z} _{n};h_{n}),\qquad {\text{ with }}\qquad \mathbf {z} _{0}=\mathbf {x} _{0},$

where

$\mathbf {\phi } (t_{n},\mathbf {z} _{n};h_{n})=\int \limits _{0}^{h_{n}}e^{\mathbf {f} _{\mathbf {x} }(\mathbf {z} _{n},\mathbf {\xi } (t_{n}))(h_{n}-u)}(\mathbf {f(z} _{n},\mathbf {\xi } (t_{n}))+\mathbf {f} _{\mathbf {\xi } }(\mathbf {z} _{n},\mathbf {\xi } (t_{n}))({\widetilde {\mathbf {\xi } }}(t_{n}+u)-{\widetilde {\mathbf {\xi } }}(t_{n})))\,du$

and ${\widetilde {\mathbf {\xi } }}$ is an approximation to the process $\mathbf {\xi }$ for all $t\in \left[t_{0},T\right].$ Here, $\mathbf {f} _{x}$ and $\mathbf {f} _{\xi }$ denote the partial derivatives of $\mathbf {f}$ with respect to $\mathbf {x}$ and $\xi$ , respectively.

Local linearization schemes

Depending on the approximations ${\widetilde {\mathbf {\xi } }}$ to the process $\mathbf {\xi }$ and of the algorithm to compute $\mathbf {\phi }$ , different Local Linearizations schemes can be defined. Every numerical implementation $\mathbf {y} _{n}$ of the local linear discretization $\mathbf {z} _{n}$ is generically called local linearization scheme.

LL schemes

\mathbf {y} _{n+1}=\mathbf {y} _{n}+\mathbf {L} (\mathbf {P} _{p,q}(2^{-k_{n}}\mathbf {M} _{n}h_{n}))^{2^{k_{n}}}\mathbf {r,} \quad

^[16]^[17]

where the matrices $\mathbf {M} _{n},\quad \mathbf {L} \quad and\quad \mathbf {r}$ are defined as

$\mathbf {M} _{n}=\left[{\begin{array}{ccc}\mathbf {f} _{\mathbf {x} }\left(\mathbf {y} _{n},\mathbf {\xi } (t_{n})\right)&\mathbf {f} _{\mathbf {\xi } }(\mathbf {y} _{n},\mathbf {\xi } (t_{n})(\mathbf {\xi } (t_{n+1})-\mathbf {\xi } (t_{n}))/h_{n}&\mathbf {f} \left(\mathbf {y} _{n},\mathbf {\xi } (t_{n})\right)\\0&0&1\\0&0&0\end{array}}\right]$

$\mathbf {L} =\left[{\begin{array}{ll}\mathbf {I} &\mathbf {0} _{d\times 2}\end{array}}\right]$ , $\mathbf {r} ^{\intercal }=\left[{\begin{array}{ll}\mathbf {0} _{1\times (d+1)}&1\end{array}}\right]$ , and p+q>1. For large systems of RDEs,^[17]

$\mathbf {y} _{n+1}=\mathbf {y} _{n}+\mathbf {L\mathbf {k} } _{m_{n},k_{n}}^{p,q}(h_{n},\mathbf {M} _{n},\mathbf {r} ),\quad p+q>1\quad and\quad m_{n}>2.$

The convergence rate of both schemes is $min\{2,2\gamma \}$ , where is $\gamma$ the exponent of the Holder condition of $\mathbf {\xi }$ .

Figure 3 presents the phase portrait of the RDE

${\frac {dx_{1}}{dt}}=-x_{2}+\left(1-x_{1}^{2}-x_{2}^{2}\right)x_{1}\sin(w^{H}(t))^{2},\quad \qquad x_{1}(0)=0.8\qquad (6.2)$

${\frac {dx_{2}}{dt}}=x_{1}+(1-x_{1}^{2}-x_{2}^{2})x_{2}\sin(w^{H}(t))^{2},\qquad \qquad x_{2}(0)=0.1,\qquad (6.3)$

and its approximation by two numerical schemes, where $w^{H}$ denotes a fractional Brownian process with Hurst exponent H=0.45.

Strong LL methods for SDEs

Consider the d-dimensional Stochastic Differential Equation (SDE)

$d\mathbf {x} (t)=\mathbf {f} (t,\mathbf {x} (t))dt+\sum \limits _{i=1}^{m}\mathbf {g} _{i}(t)d\mathbf {w} ^{i}(t),\quad t\in \left[t_{0},T\right],\qquad \qquad \qquad (7.1)$

with initial condition $\mathbf {x} (t_{0})=\mathbf {x} _{0}$ , where the drift coefficient $\mathbf {f}$ and the diffusion coefficient $\mathbf {g} _{i}$ are differentiable functions, and $\mathbf {w=(\mathbf {w} } ^{1},\ldots ,\mathbf {w} ^{m}\mathbf {)}$ is an m-dimensional standard Wiener process.

Local linear discretization

For a time discretization $\left(t\right)_{h}$ , the order- $\mathbb {\gamma }$ (=1,1.5) Strong Local Linear discretization of the solution of the SDE (7.1) is defined by the recursive relation ^[18]^[19]

$\mathbf {z} _{n+1}=\mathbf {z} _{n}+\mathbf {\phi } _{\mathbb {\gamma } }(t_{n},\mathbf {z} _{n};h_{n})+\mathbf {\xi } (t_{n},\mathbf {z} _{n};h_{n}),\quad with\quad \mathbf {z} _{0}=\mathbf {x} _{0},$

where

$\mathbf {\phi } _{\mathbb {\gamma } }(t_{n},\mathbf {z} _{n};\delta )=\int _{0}^{\delta }e^{\mathbf {f} _{\mathbf {x} }(t_{n},\mathbf {y} _{n})(\delta -u)}(\mathbf {f(} t_{n},\mathbf {z} _{n})+\mathbf {a} ^{\mathbb {\gamma } }(t_{n},\mathbf {z} _{n})u)du$

and

$\mathbf {\xi } \left(t_{n},\mathbf {z} _{n};\delta \right)=\sum \limits _{i=1}^{m}\int \nolimits _{t_{n}}^{t_{n}+\delta }e^{\mathbf {f} _{\mathbf {x} }(t_{n},\mathbf {z} _{n})(t_{n}+\delta -u)}\mathbf {g} _{i}(u)\,d\mathbf {w} ^{i}(u).$

Here,

$\mathbf {a} ^{\mathbb {\gamma } }(t_{n},\mathbf {z} _{n})=\left\{{\begin{array}{cl}\mathbf {f} _{t}(t_{n},\mathbf {z} _{n})&{\text{for }}\qquad \mathbb {\gamma } =1\\\mathbf {f} _{t}(t_{n},\mathbf {z} _{n})+{\frac {1}{2}}\sum \limits _{j=1}^{m}(\mathbf {I} \otimes \mathbf {g} _{j}^{\intercal }(t_{n}))\mathbf {f} _{\mathbf {xx} }(t_{n},\mathbf {z} _{n})\mathbf {g} _{j}(t_{n})&{\text{for }}\quad \mathbb {\gamma } =1.5,\end{array}}\right.$

$\mathbf {f} _{\mathbf {x} },\mathbf {f} _{t}$ denote the partial derivatives of $\mathbf {f}$ with respect to the variables $\mathbf {x}$ and t, respectively, and $\mathbf {f} _{\mathbf {xx} }$ the Hessian matrix of $\mathbf {f}$ with respect to $\mathbf {x}$ . The strong Local Linear discretization $\mathbf {z} _{n+1}$ converges with order $\mathbb {\gamma }$ (= 1, 1.5) to the solution of (7.1).

High-order local linear discretizations

After the local linearization of the drift term of (7.1) at $(t_{n},\mathbf {z} _{n})$ , the equation for the residual $\mathbf {r}$ is given by

$d\mathbf {r} (t)=\mathbf {q} _{\gamma }(t_{n},\mathbf {z} _{n};t\mathbf {,\mathbf {r} } (t))\,dt+\sum \limits _{i=1}^{m}\mathbf {g} _{i}(t)\,d\mathbf {w} ^{i}(t)\mathbf {,} \qquad \mathbf {r} (t_{n})=\mathbf {0}$

for all $t\in \lbrack t_{n},t_{n+1}]$ , where

$\mathbf {q} _{\gamma }(t_{n},\mathbf {z} _{n};s\mathbf {,\xi } )=\mathbf {f} (s,\mathbf {z} _{n}+\mathbf {\phi } _{\gamma }(t_{n},\mathbf {z} _{n};s-t_{n})+\mathbf {\xi } )-\mathbf {f} _{\mathbf {x} }(t_{n},\mathbf {z} _{n})\mathbf {\phi } _{\gamma }(t_{n},\mathbf {z} _{n};s-t_{n})-\mathbf {a} ^{\gamma }(t_{n},\mathbf {z} _{n})(s-t_{n})-\mathbf {f} (t_{n},\mathbf {z} _{n}).$

A high-order local linear discretization of the SDE (7.1) at each point $t_{n+1}\in (t)_{h}$ is then defined by the recursive expression ^[20]

$\mathbf {z} _{n+1}=\mathbf {z} _{n}+\mathbf {\phi } _{\gamma }(t_{n},\mathbf {z} _{n};h_{n})+{\widetilde {\mathbf {r} }}(t_{n},\mathbf {z} _{n};h_{n}),\qquad {\text{ with }}\qquad \mathbf {z} _{0}=\mathbf {x} _{0},$

where ${\widetilde {\mathbf {r} }}$ is a strong approximation to the residual $\mathbf {r}$ of order $\alpha$ higher than 1.5. The strong HOLL discretization $\mathbf {z} _{n+1}$ converges with order $\alpha$ to the solution of (7.1).

Local linearization schemes

Depending on the way of computing $\mathbf {\phi } _{\mathbb {\gamma } }$ , $\mathbf {\xi }$ and ${\widetilde {\mathbf {r} }}$ different numerical schemes can be obtained. Every numerical implementation $\mathbf {y} _{n}$ of a strong Local Linear discretization $\mathbf {z} _{n}$ of any order is generically called Strong Local Linearization (SLL) scheme.

Order 1 SLL schemes

$\mathbf {y} _{n+1}=\mathbf {y} _{n}+\mathbf {L} (\mathbf {P} _{p,q}(2^{-k_{n}}\mathbf {M} _{n}h_{n}))^{2^{k_{n}}}\mathbf {r+} \sum \limits _{i=1}^{m}\mathbf {g} _{i}(t_{n})\Delta \mathbf {w} _{n}^{i},\quad$ ^[21] $\qquad \qquad (7.2)$

where the matrices $\mathbf {M} _{n}$ , $\mathbf {L}$ and $\mathbf {r}$ are defined as in (4.6), $\Delta \mathbf {w} _{n}^{i}$ is an i.i.d. zero mean Gaussian random variable with variance $h_{n}$ , and p + q > 1. For large systems of SDEs,^[21] in the above scheme $(\mathbf {P} _{p,q}(2^{-k_{n}}\mathbf {M} _{n}h_{n}))^{2^{k_{n}}}\mathbf {r}$ is replaced by $\mathbf {\mathbf {k$

[1]

[2]

[3]

[4]

[5]

[6]

[7]

[8]

[9]

[10]

[11]

[12]

[13]

[14]

[15]

[16]

[17]

[18]

[19]

[20]

[21]