Local linearization method

Local Linearization Method

In numerical analysis, the Local Linearization (LL) method is a general strategy for designing numerical inte-

grators 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 development for a variety of

equations such that 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 deals with complex models in a variety of fields as neuroscience, finance, forestry

management, control engineering, mathematical statistics, etc.

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 an high order approximation of the

residual ${\displaystyle \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.

This sandbox is in the article namespace. Either move this page into your userspace, or remove the {{User sandbox}} template.

LL methods for ODEs

Consider the d-dimensional Ordinary Differential Equation (ODE).

${\displaystyle {\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 (1).}$ 

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

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

${\displaystyle \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

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

results from the linear approximation, and

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

is the residual of the linear approximation. Here, ${\displaystyle \mathbf {f} _{\mathbf {x} }}$ and ${\displaystyle \mathbf {f} _{t}}$  denote the partial derivatives of f with respect to the variables x and t, respectively, and ${\displaystyle \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}\left(t_{n},\mathbf {z} _{n}\right)(s-t_{n})-\mathbf {f} \left(t_{n},\mathbf {z} _{n}\right)+\mathbf {f} _{\mathbf {x} }(t_{n},\mathbf {z} _{n})\mathbf {z} _{n}}$ 

Local Linear discretization

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

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

The Local Linear discretization (3) converges with order 2 to the solution of nonlinear ODEs, but it match the

solution of the linear ODEs. The recursion (3) is also known as Exponential Euler discretization.

High Order Local Linear discretizations

For a time discretization ${\displaystyle \left(t\right)_{h},}$  a High Order Local Linear (HOLL) discretization of the ODE (1) at each point ${\displaystyle t_{n+1}\in \left(t\right)_{h}}$  is deffined by the recursive expression

${\displaystyle \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 with\quad \mathbf {z} _{0}=\mathbf {x} _{0},\qquad \qquad \qquad (4)}$ 

where ${\displaystyle {\tilde {r}}}$  is an approximation to the residual r of order ${\displaystyle \alpha }$  higher than 2 ${\displaystyle (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 }).}$  The HOLL discretization (4) converges with order to the solution of nonlinear ODEs, but it match the solution of the linear ODEs.

HOLL discretizations can be derived in two ways: 1) by approximating the integral representation (2) of r; and 2) by using a numerical integrator for the di§erential representation of r deffined by

${\displaystyle {\frac {d\mathbf {r} \left(t\right)}{dt}}=\mathbf {q} (t_{n},\mathbf {z} _{n};t\mathbf {,\mathbf {r} } \left(t\right)\mathbf {),} \qquad with\qquad \mathbf {r} \left(t_{n}\right)=\mathbf {0,} \qquad \qquad \qquad (5)}$ 

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

${\displaystyle \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 } \left(t_{n},\%\mathbf {z} _{n};s-t_{n}\right)-\mathbf {f} _{t}\left(t_{n},\mathbf {z} \%_{n}\right)(s-t_{n})-\mathbf {f} \left(t_{n},\mathbf {z} _{n}\right).}$ 

The resulting approximation is often called Locally Linearized discretization.

Known HOLL discretizations are the following.

Locally Linearized Runge Kutta discretization

${\displaystyle \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 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 (5) via a s-stage RK scheme with coefficients ${\displaystyle \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

${\textstyle \mathbf {z} _{n+1}=\mathbf {z} _{n}+\mathbf {\phi } (t_{n},\mathbf {z} \%_{n};h_{n})+\int _{0}^{h_{n}}e^{\left(h_{n}-s\right)\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} \left(t\right)}{dt^{j+1}}}-\mathbf {f} _{\mathbf {x} }\left(t_{n},\mathbf {z} \%_{n}\right){\frac {d^{j}\mathbf {x} \left(t\right)}{dt^{j}}}\right)\mid _{t=\%\mathbf {z} _{n}},}$ 

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

Exponential Rosembrock discretization (poner link) is obtained by approximating the integral (2) by aquadrature rule.

Linealized Exponential Adams discretization

${\textstyle \mathbf {z} _{n+1}=\mathbf {z} _{n}+\mathbf {\phi } (t_{n},\mathbf {z} _{n};h_{n})+h_{n}\sum _{j=1}^{p}\sum _{l=1}^{j}{\frac {\gamma _{j+1}}{l}}\nabla ^{l}\mathbf {g} _{n}(t_{n},\mathbf {z} _{n}),\quad with\quad \gamma _{j+1}=(-1)^{j+1}\int \limits _{0}^{1}e^{(1-\theta )h_{n}\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 ${\displaystyle \mathbf {g} _{n}}$ in (2) by a Hermite polynomial of degree p, where ${\displaystyle \nabla ^{l}\mathbf {g} \%_{n}(t_{m},\mathbf {z} _{m})}$  denotes the l-th backward di§erence of ${\displaystyle \mathbf {g} _{n}(t_{m},\mathbf {z} _{m})}$ .

Local Linearization schemes

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

integrals ${\displaystyle \phi _{j}}$  of the form

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

where A is an d ${\displaystyle \times }$  d matrix. Every numerical implementation ${\displaystyle \mathbf {y} _{n}}$  of a Local Linear discretization ${\displaystyle \mathbf {z} _{n}}$  of any order is generically called Local Linearization scheme.

Computing integrals involving matrix exponential

Among a number of algorithms to compute the integrals ${\displaystyle \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

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

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

${\displaystyle \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+lr)\times (d+lr)},}$ 

${\displaystyle \mathbf {L} =[\mathbf {I} \mathbf {0} _{d\times l)}],\mathbf {r} =[\mathbf {0} _{1\times (d+l-1)}\quad 1]^{\intercal }\quad ,and\quad \mathbf {v} _{i}=\mathbf {a} _{i}(i-1)!}$ 

If ${\displaystyle \mathbf {P} _{p,q}(2^{-k}\mathbf {H} h)}$  denotes the (p; q)-Padé approximation of ${\displaystyle e^{2^{-k}\mathbf {H} h}}$  and k is the smallest integer number such that ${\displaystyle \left\Vert 2^{-k}\mathbf {H} h\right\Vert \leq {\frac {1}{2}},}$ 

${\displaystyle \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 ${\displaystyle \mathbf {\mathbf {k} } _{m,k}^{p,q}(h,\mathbf {H} ,\mathbf {r} )}$  denotes the (m; p; q; k) Krylov-Padé approximation of ${\displaystyle e^{h\mathbf {H} }\mathbf {r} }$ 

${\displaystyle \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}}.}$ 

Order 2 LL schemes

${\displaystyle \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,} \qquad \qquad (6)}$ 

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

${\displaystyle \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)},}$ 

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

${\displaystyle \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 with\qquad m_{n}>1.}$ 

Order 3 LL-Taylor schemes

${\displaystyle \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 {,} \qquad \qquad (7)}$ 

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

${\displaystyle \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)},}$ 

${\displaystyle \mathbf {L} _{1}=\left[{\begin{array}{ll}\mathbf {I} _{d}&\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, ${\displaystyle \mathbf {f} _{\mathbf {xx} }}$  denotes the second derivative of f with respect to x,

I the d-dimensional identity matrix, and p + q > 2. For large systems of ODEs.

${\displaystyle \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 with\qquad m_{n}>2.}$ 

Order 4 LL-RK schemes

${\displaystyle \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}),\qquad \qquad (8)}$ 

where

${\displaystyle \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

${\displaystyle \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 ${\displaystyle \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 ${\displaystyle \mathbf {u} _{j}}$  in the above

scheme is replaced by ${\displaystyle \mathbf {u} _{j}=\mathbf {L\mathbf {k} } _{m_{j},k_{j}}^{p,q}(c_{j}h_{n},\mathbf {M} _{n},\mathbf {r} ).}$ 

Locally Linearized Runge-Kutta of Dormand & Prince

${\displaystyle \mathbf {y} _{n+1}=\mathbf {y} _{n}+\mathbf {u} _{s}+h_{n}\sum _{j=1}^{s}b_{j}\mathbf {k} _{j}\qquad 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},\qquad \qquad (9)}$ 

where s = 6 is the number of the stages,

${\displaystyle \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 ${\displaystyle \mathbf {k} _{1}\equiv \mathbf {0} }$ , and ${\displaystyle 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. For large systems of ODEs, the vector ${\displaystyle \mathbf {u} _{j}}$  in the above scheme is replaced by ${\displaystyle \mathbf {u} _{j}=\mathbf {L\mathbf {k} } _{m_{j},k_{j}}^{p,q}(c_{j}h_{n},\mathbf {M} _{n},\mathbf {r} ).}$ 

 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${\displaystyle p\leq q\leq p+2}$ , the LL schemes (6)-(9) are A-stable.

With q = p + 1 or q = p + 2, the LL schemes (6)-(9) are also L-stable. For linear ODEs, the LL schemes

(6)-(9) converge with order p + q. Moreover, with p = q = 6 and mn = 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. This includes stiff and highly oscillatory linear equations. Moreover, the LL

schemes (6)-(9) are regular for linear ODEs and inherit the symplectic structure of Hamiltonian harmonic

oscillators. 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.

LL methods for DDEs

Consider the d-dimensional Delay Differential Equation (DDE)

${\displaystyle {\frac {d\mathbf {x} \left(t\right)}{dt}}=\mathbf {f} \left(t,\mathbf {x} \left(t\right),\mathbf {x} _{t}(-\tau _{1}),\cdots ,\mathbf {x} _{t}(-\tau _{m})\right),\qquad t\in \left[t_{0},T\right],\qquad \qquad (10)}$ 

with m constant delays ${\displaystyle \tau _{i}>0}$  and initial condition ${\displaystyle \mathbf {x} _{t_{0}}(s)=\mathbf {\varphi } (s)}$  for all ${\displaystyle s\in \left[-\tau ,0\right],}$  where f is a differentiable

function, ${\displaystyle \mathbf {x} _{t}:\left[-\tau ,0\right]\longrightarrow \mathbb {R} ^{d}}$  is the segment function deffined as

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

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

Local Linear discretization

For a time discretization ${\displaystyle \left(t\right)_{h}}$  , the Local Linear discretization of the DDE (10) at each point ${\displaystyle t_{n+1}\in \left(t\right)_{h}}$  is deffined by the recursive expression

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

where

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

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

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

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

Here,

${\displaystyle \mathbf {A} _{n}=\mathbf {f} _{x}(t_{n},\mathbf {z} _{n},{\widetilde {\mathbf {z} }}_{t_{n}}^{1}(-\tau _{1}),...,{\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}),...,{\widetilde {\mathbf {z} }}_{t_{n}}^{m}(-\tau _{d}))}$ 

are constant matrices and

${\displaystyle \mathbf {c} _{n}=\mathbf {f} _{t}(t_{n},\mathbf {z} _{n},{\widetilde {\mathbf {z} }}_{t_{n}}^{1}(-\tau _{1}),...,{\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}),...,{\widetilde {\mathbf {z} }}_{t_{n}}^{m}(-\tau _{d}))}$ 

are constant vectors. ${\displaystyle \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, x and ${\displaystyle \mathbf {x} _{t}(-\tau _{i})}$  The Local Linear discretization (11) converges to the solution of (10) with order ${\displaystyle \alpha =\min\{2,r\},}$  if ${\displaystyle {\widetilde {\mathbf {z} }}_{t_{n}}^{i}}$  approximates ${\displaystyle \mathbf {z} _{t_{n}}^{i}}$  with order ${\displaystyle r(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 ${\displaystyle u\in \lbrack 0,h_{n}])}$ .

Local Linearization schemes

Depending of the approximations ${\displaystyle {\widetilde {\mathbf {z} }}_{t_{n}}^{i}}$  and of the algorithm to compute ${\displaystyle \mathbf {\phi } }$  di§erent Local Linearizations schemes

can be deÖned. Every numerical implementation ${\displaystyle \mathbf {y} _{n}}$  of a Local Linear discretization ${\displaystyle \mathbf {z} _{n}}$  is generically called Local Linearization scheme.

Order 2 Polynomial LL schemes

${\displaystyle \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,} \qquad \qquad (12)}$ 

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

${\displaystyle \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_{1\times d}&0&1\\0_{1\times d}&0&0\end{bmatrix}}\in \mathbb {R} ^{(d+2)\times (d+2)},}$ 

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

${\displaystyle \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} ,{\text{ 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 (10) deÖned through the LL scheme (12) for all ${\displaystyle t\in \lbrack t_{0},t_{n}]}$  and by ${\displaystyle \mathbf {y} \left(t\right)=\mathbf {\varphi } \left(t\right)}$  for ${\displaystyle t\in \left[t_{0}-\tau ,t_{0}\right]}$ 

For large systems of DDEs

${\displaystyle \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 ${\displaystyle p+q>1\quad and\quad m_{n}>2.}$ 

LL methods for RDEs

Consider the d-dimensional Random Differential Equation (RDE)

${\displaystyle {\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 (13).}$ 

with initial condition ${\displaystyle \mathbf {x} (t_{0})=\mathbf {x} _{0},}$  where ${\displaystyle \mathbf {\xi } }$  is a k-dimensional separable finite continuous stochastic process, and

f is a differentiable function. Suppose that a realization (path) of ${\displaystyle \mathbf {\xi } }$  is given.

Local Linear discretization

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

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

where

${\displaystyle \mathbf {\phi } (t_{n},\mathbf {z} _{n};h_{n})=\int \limits _{0}^{h_{n}}e^{\mathbf {f} _{\mathbf {x} }\left(\mathbf {z} _{n},\mathbf {\xi } (t_{n})\right)(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 ${\displaystyle {\widetilde {\mathbf {\xi } }}}$  is an approximation to the process ${\displaystyle \mathbf {\xi } }$  for all ${\displaystyle t\in \left[t_{0},T\right].}$  Here, ${\displaystyle \mathbf {f} _{x}}$  and ${\displaystyle \mathbf {f} _{\xi }}$  denote the partial derivatives with respect to ${\displaystyle \mathbf {x} }$  and ${\displaystyle \xi }$ , respectively.

Local Linearization schemes

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

LL schemes

${\displaystyle \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,} }$  where the matrices ${\displaystyle \mathbf {M} _{n},\quad \mathbf {L} \quad and\quad \mathbf {r} }$  are defined as

${\displaystyle \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]\in \mathbb {R} ^{(d+2)\times (d+2)},}$ 

${\displaystyle \mathbf {L} =\left[{\begin{array}{ll}\mathbf {I} _{d}&\mathbf {0} _{d\times 2}\end{array}}\right]}$ , ${\displaystyle \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,

${\displaystyle \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 ${\displaystyle min\{2,2\gamma \}}$ , where is ${\displaystyle \gamma }$  the exponent of the Holder condition of ${\displaystyle \mathbf {\xi } }$ .

Strong LL methods for SDEs

Consider the d-dimensional Stochastic Differential Equation (SDE)

${\displaystyle 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 (14)}$ 

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

Local Linear discretization

For a time discretization ${\displaystyle \left(t\right)_{h}}$  , the order-${\displaystyle \mathbb {\gamma } }$  (=1,1.5) Strong Local Linear discretization of the solution of

the SDE (14) is defined by the recursive relation.

${\displaystyle \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

${\displaystyle \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

${\displaystyle \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,

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

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

High Order Local Linear discretizations

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

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

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

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

High Order Local Linear discretization of the SDE (14) at each point ${\displaystyle t_{n+1}\in \left(t\right)_{h}}$  is then defined by the

recursive expression.

${\displaystyle \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 with\qquad \mathbf {z} _{0}=\mathbf {x} _{0},}$ 

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

Local Linearization schemes

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

Order 1 SLL schemes

${\displaystyle \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},(15)}$ 

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

Order 1.5 SLL schemes

${\displaystyle \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}\left(\mathbf {g} _{i}(t_{n})\Delta \mathbf {w} _{n}^{i}+\mathbf {f} _{\mathbf {x} }(t_{n},{\widetilde {\mathbf {y} }}_{n})\mathbf {g} _{i}(t_{n})\Delta \mathbf {z} _{n}^{i}+{\frac {d\mathbf {g} _{i}(t_{n})}{dt}}(\Delta \mathbf {w} _{n}^{i}h_{n}-\Delta \mathbf {z} _{n}^{i})\right),(16)}$ 

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

${\displaystyle \mathbf {M} _{n}={\begin{bmatrix}\mathbf {f} _{\mathbf {x} }(t_{n},\mathbf {y} _{n})&\mathbf {f} _{t}(t_{n},\mathbf {y} _{n})\quad +{\frac {1}{2}}\sum \limits _{j=1}^{m}\left(\mathbf {I} _{d\times d}\otimes \mathbf {g} _{j}^{\intercal }\left(t_{n}\right)\right)\mathbf {f} _{\mathbf {xx} }(t_{n},\mathbf {y} _{n})\mathbf {g} _{j}\left(t_{n}\right)&\mathbf {f} (t_{n},\mathbf {y} _{n})\\0&0&1\\0&0&0\end{bmatrix}}\in \mathbb {R} ^{(d+2)\times (d+2)},}$ 

${\displaystyle \mathbf {L} =\left[{\begin{array}{ll}\mathbf {I} _{d}&\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]}$ , ${\displaystyle \Delta \mathbf {z} _{n}^{i}}$  is a i.i.d. zero mean Gaussian random variable with variance ${\displaystyle E\left((\Delta \mathbf {z} _{n}^{i})^{2}\right)=\%{\frac {1}{3}}h_{n}^{3}}$  and covariance ${\displaystyle E(\Delta \mathbf {w} _{n}^{i}\Delta \mathbf {z} _{n}^{i})={\frac {1}{2}}h_{n}^{2}}$  and p+q>1. For large systems of SDEs, in the above scheme ${\displaystyle (\mathbf {P} _{p,q}(2^{-k_{n}}\mathbf {M} _{n}h_{n}))^{2^{k_{n}}}\mathbf {r} }$  is replaced by ${\displaystyle \mathbf {\mathbf {k} } _{m_{n},k_{n}}^{p,q}(h_{n},\mathbf {M} _{n}^{\gamma },\mathbf {r} )}$ .

Order 2 SLL-Taylor schemes

${\displaystyle \mathbf {y} _{t_{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 _{j=1}^{m}\mathbf {g} _{j}\left(t_{n}\right)\Delta \mathbf {w} _{n}^{j}+\sum \limits _{j=1}^{m}\mathbf {f} _{\mathbf {x} }(t_{n},\mathbf {y} _{n})\mathbf {g} _{j}\left(t_{n}\right){\widetilde {J}}_{\left(j,0\right)}+\sum \limits _{j=1}^{m}{\frac {d\mathbf {g} _{_{j}}}{dt}}\left(t_{n}\right){\widetilde {J}}_{\left(0,j\right)}}$ 

${\displaystyle \qquad \qquad +\sum \limits _{j_{1},j_{2}=1}^{m}\left(\mathbf {I} _{d\times d}\otimes \mathbf {g} _{j_{2}}^{\intercal }\left(t_{n}\right)\right)\mathbf {f} _{\mathbf {xx} }(t_{n},\mathbf {y} _{n})\mathbf {g} _{j_{1}}\left(t_{n}\right){\widetilde {J}}_{\left(j_{1},j_{2},0\right)},}$  (17)

where ${\displaystyle \mathbf {M} _{n}}$ , ${\displaystyle \mathbf {L} }$ , ${\displaystyle \mathbf {r} }$  and ${\displaystyle \Delta \mathbf {w} _{n}^{i}}$  are defined as in the order-1 SLL schemes, and ${\displaystyle {\widetilde {J}}_{\alpha }}$  is order-2 approximation to the multiple Stratonovish integral ${\displaystyle J_{\alpha }}$ .

Order 2 SLL-RK schemes

For SDEs with a single Wiener noise (m=1)

${\displaystyle \mathbf {y} _{t_{n+1}}=\mathbf {y} _{n}+{\widetilde {\mathbf {\phi } }}(t_{n},\mathbf {y} _{n};h_{n})+{\frac {h_{n}}{2}}\left(\mathbf {k} _{1}+\mathbf {k} _{2}\right)+\mathbf {g} \left(t_{n}\right)\Delta w_{n}+{\frac {\left(\mathbf {g} \left(t_{n+1}\right)-\mathbf {g} \left(t_{n}\right)\right)}{h_{n}}}J_{\left(0,1\right)}.\qquad \qquad (18)}$ 

where

${\displaystyle \mathbf {k} _{1}=\mathbf {f} (t_{n}+{\frac {h_{n}}{2}},\mathbf {y} _{n}+{\widetilde {\mathbf {\phi } }}(t_{n},\mathbf {y} _{n};{\frac {h_{n}}{2}})+\gamma _{+})-\mathbf {f} _{\mathbf {x} }(t_{n},\mathbf {y} _{n}){\widetilde {\mathbf {\phi } }}(t_{n},\mathbf {y} _{n};{\frac {h_{n}}{2}})-\mathbf {f} \left(t_{n},\mathbf {y} _{n}\right)-\mathbf {f} _{t}\left(t_{n},\mathbf {y} _{n}\right){\frac {h_{n}}{2}},}$ 

${\displaystyle \mathbf {k} _{2}=\mathbf {f} (t_{n}+{\frac {h_{n}}{2}},\mathbf {y} _{n}+{\widetilde {\mathbf {\phi } }}(t_{n},\mathbf {y} _{n};{\frac {h_{n}}{2}})+\gamma _{-})-\mathbf {f} _{\mathbf {x} }(t_{n},\mathbf {y} _{n}){\widetilde {\mathbf {\phi } }}(t_{n},\mathbf {y} _{n};{\frac {h_{n}}{2}})-\mathbf {f} \left(t_{n},\mathbf {y} _{n}\right)-\mathbf {f} _{t}\left(t_{n},\mathbf {y} _{n}\right){\frac {h_{n}}{2}},}$ 

with

${\displaystyle \gamma _{\pm }={\frac {1}{h_{n}}}\mathbf {g} \left(t_{n}\right)\left\{{\widetilde {J}}_{\left(1,0\right)}\pm {\sqrt {2{\widetilde {J}}_{\left(1,1,0\right)}h_{n}-{\widetilde {J}}_{\left(1,0\right)}^{2}}}\right\}.}$ 

Here, ${\displaystyle {\widetilde {\mathbf {\phi } }}(t_{n},\mathbf {y} _{n};h_{n})=\mathbf {L} (\mathbf {P} _{p,q}(2^{-k_{n}}\mathbf {M} _{n}h_{n}))^{2^{k_{n}}}\mathbf {r} }$  for for low dimensional SDEs, and ${\displaystyle {\widetilde {\mathbf {\phi } }}(t_{n},\mathbf {y} _{n};h_{n})=\mathbf {L\mathbf {k} } _{m_{n},k_{n}}^{p,q}(h_{n},\mathbf {M} ,\mathbf {r} )}$  for large systems of SDEs, where ${\displaystyle \mathbf {M} _{n}}$ , ${\displaystyle \mathbf {L} }$ , ${\displaystyle \mathbf {r} }$ , ${\displaystyle \Delta \mathbf {w} _{n}^{i}}$  and ${\displaystyle {\widetilde {J}}_{\alpha }}$  are defined as in the order-2 SLL-Taylor schemes, p+q>1 and ${\displaystyle m_{n}>2}$ .

Stability and dynamics

By construction the strong LL and HOLL discretizations inherit the stability and dynamics of the linear ODEs, but it is not the case of the strong LL schemes in general. LL schemes (15)-(18) with ${\displaystyle p\leq q\leq p+2}$  are A-stable, which includes stiff and highly oscillatory linear equations. Moreover, for linear SDEs with random attractors, these schemes also have a random attractor that converges in probability to the exact one as stepsizes decrease and preserve the ergodicity of these equations for any stepsize. These schemes also reproduce essential dynamical properties of simple and coupled harmonic oscillators such as the linear growth of energy along the paths, the oscillatory behavior around 0, the symplectic structure of Hamiltonian oscillators, and the mean of the paths. For nonlinear SDEs with small noise (e.g., (14) with ${\displaystyle \mathbf {g} _{i}(t)\approx 0}$ ), the paths of these SLL schemes are basically the nonrandom paths of the LL scheme (6) for ODEs plus a small disturbance related to the small noise. In this situation, the dynamical properties of that deterministic scheme, such as the linearization preserving and the preservation of the exact solution dynamics around hyperbolic equilibrium points and periodic orbits, become relevant for the paths of the SLL scheme.

Weak LL methods for SDEs