Local linearization method: Difference between revisions

Revision as of 15:13, 7 July 2017

Latest revision as of 21:00, 14 April 2025

(532 intermediate revisions by 35 users not shown)

Line 1:

{{short description|Numerical method for differential equations}}

= Local Linearization Method =

In [[numerical analysis]], the Local'''local Linearizationlinearization (LL) method''' is a general strategy for designing [[Numerical integration|numerical inte]][[Numerical integration|gratorsintegrators]] 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 developmentdeveloped for a variety of equations such thatas the [[Ordinary differential equation|ordinary]], [[Delay differential equation|delayed]], random and [[Stochastic differential equation|stochastic]] differential equations. The LL integrators are key component in the implementation of [[Estimation theory|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 dealsdeal with complex models in a variety of fields as [[neuroscience]], [[finance]], [[Forestry|forestry]] [[Forestry|management]], [[Control engineering|control engineering,]], [[mathematical statistics]], etc.

== Background ==

= 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 [[Dynamical system|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 <math>\mathbf{r}=

\mathbf{x}-\mathbf{z}</math>.

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.

= Local Linearization scheme =

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

== 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 theory|stability]] and [[Dynamical system|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 <math>\mathbf{r}=

\mathbf{x}-\mathbf{z}</math>.

== Local linearization scheme ==

{{User sandbox|Lilian|Local Linearization Method}}

A '''''Local Linearization (LL) scheme''''' is the final [[Recursion (computer science)|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|Ordinary Differential Equation]] (ODE).

<div style="text-align: center;">

<math>

\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)

</math>

</div>

Line 23 ⟶ 25:

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

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

<div style="text-align: center;">

<math>\mathbf{x}(t_{n}t_n+h)=\mathbf{x}(t_{n}t_n)+\mathbf{\phi }(t_t_n,\mathbf{nx}(t_n);h)+\mathbf{r}(t_n,\mathbf{x}(t_n);h),

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

</math>

</div>

Line 34 ⟶ 35:

<div style="text-align: center;">

<math>\mathbf{\phi }(t_{n}t_n,\mathbf{z}_{n}_n;h)=\int\limits_{0}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

\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</math>

</div>

Line 43 ⟶ 42:

<div style="text-align: center;">

<math>\mathbf{r}(t_{n}t_n,\mathbf{z}_{n}_n;h)=\int\limits_{0}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)</math>

x}}\left( t_{n},\mathbf{z}_{n}\right) (h-s)}\mathbf{g}_{n}(s,\mathbf{x}%

(t_{n}+s))ds, \qquad \qquad \qquad (2)</math>

</div>

is the residual of the linear approximation. Here, <math>\mathbf{f}_{\mathbf{x}}</math> and <math>\mathbf{f}_{t}</math> denote the partial derivatives of '''f''' with respect to the variables '''x''' and '''t''', respectively, and <math>\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. </math>

_{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}</math>.

=== Local Linearlinear discretization ===

For a time discretization <math>\left( t\right) _{h}</math>, the ''Local Linear discretization'' of the ODE (4.1) at each point <math>t_{n+1}\in \left(

t\right) _{h}</math> is defined by the recursive expression <ref name=":8">Jimenez J.C. (2009). "Local Linearization methods for the numerical integration of ordinary differential equations: An overview". [https://inis.iaea.org/search/search.aspx?orig_q=RN:40101978 ICTP Technical Report]. 035: 357–373.</ref><ref name=":18" />

t\right) _{h}</math> is deffined by the recursive expression

<div style="text-align: center;">

<math>\mathbf{z}_{n+1}=\mathbf{z}_{n}_n+\mathbf{\phi }(t_{n}t_n,\mathbf{z}_{n}_n;h_{n}h_n),

\qquad \text{ with } \quad \mathbf{z}_{0}_0=\mathbf{x}_{0}\text{_0.} \qquad \qquad \qquad \qquad (4.3)</math>

</div>

The Local Linear discretization (4.3) [[Rate of convergence|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.<ref name=":22" />

=== High-order Orderlocal Local Linearlinear discretizations ===

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

For a time discretization <math>(t)_h,</math> a ''high-order local linear (HOLL)'' discretization of the ODE (4.1) at each point <math>t_{n+1} \in (t)_h</math> is defined by the recursive expression <ref name= ":8"/><ref name= ":3"/><ref name= ":2"/><ref name=":1" />

<div style="text-align: center;">

<math>\mathbf{z}_{n+1}=\mathbf{z}_{n}_n+\mathbf{\phi }(t_t_n,\mathbf{nz}_n;h_n) + \widetilde{\mathbf{r}}(t_n,\mathbf{z}__n;h_n),\qquad \text{n with };h_ \quad \mathbf{nz}_0=\mathbf{x}_0, \qquad \qquad \qquad(4.4)+</math>

\widetilde{\mathbf{r}}(t_{n},\mathbf{z}_{n};h_{n}),\qquad with \quad

\mathbf{z}_{0}=\mathbf{x}_{0}, \qquad \qquad \qquad(4)</math>

</div>

where <math>\tilde{r}</math> is an order <math>\alpha </math> (>&nbsp;'''2''') approximation to the residual '''r''' of order <math>\alpha </math> higher than 2 <math>(i.e., \left\vert \mathbf{r}(t_t_n, \mathbf{nz}_n;h)-\widetilde{\mathbf{r}}(t_n,\mathbf{z}_n;h)\right\vert \propto h^{\alpha+1 }).</math> The HOLL discretization (4.4) [[Rate of convergence|converges]] with order <math>\alpha</math> to the solution of nonlinear ODEs, but it match the solution of the linear ODEs.

\mathbf{z}_{n};h)-\widetilde{\mathbf{r}}(t_{n},\mathbf{z}_{n};h)\right\vert

\propto h^{\alpha }).</math> The HOLL discretization (4) [[Rate of convergence|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:<ref name= ":8"/><ref name= ":3"/><ref name= ":2"/><ref name=":1" /> 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''' deffineddefined by

<div style="text-align: center;">

<math>\frac{d\mathbf{r}\left( t\right) }{dt} = \mathbf{q}(t_{n}t_n,\mathbf{z}__n;t \mathbf{n,\mathbf{r};}(t) \mathbf), \qquad \text{ with } \qquad \mathbf{r}(t_n) =\mathbf{0,} \qquad \qquad \qquad (4.5)</math>

\mathbf{,\mathbf{r}}\left( t\right) \mathbf{),}\qquad with \qquad \mathbf{r}

\left( t_{n}\right) =\mathbf{0,} \qquad \qquad \qquad (5)</math>

</div>

for all <math>t\in \lbrack t_{k}t_k,t_{k+1}]</math>, where

<div style="text-align: center;">

<small><math>\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}t_n,

\mathbf{z}_{n}_n;s-t_{n}\rightt_n) -\mathbf{f}__t( t_n,\mathbf{tz}\left(_n) t_(s-t_n)-\mathbf{nf}( t_n,\mathbf{z}_n).</math>

</div>

_{n}\right) (s-t_{n})-\mathbf{f}\left( t_{n},\mathbf{z}_{n}\right) .</math></small>

HOLL discretizations are, for instance, the followings:

* ''Locally Linearized Runge Kutta discretization''<ref name=":1">de la Cruz H.; Biscay R.J.; Carbonell F.; Jimenez J.C.; Ozaki T. (2006). "Local Linearization-Runge Kutta (LLRK) methods for solving ordinary differential equations". Lecture Note in Computer Sciences 3991: 132–139, Springer-Verlag. [[doi:10.1007/11758501 22]]. {{ISBN|978-3-540-34379-0}}.</ref><ref name=":3">de la Cruz H.; Biscay R.J.; Jimenez J.C.; Carbonell F. (2013). "Local Linearization - Runge Kutta Methods: a class of A-stable explicit integrators for dynamical systems". Math. Comput. Modelling. 57 (3–4): 720–740. [[doi:10.1016/j.mcm.2012.08.011]].</ref>

<div style="text-align: left;">

<math>\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), </math>

</div>

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

The resulting approximation is often called Locally Linearized discretization.

* ''Local linear Taylor discretization''<ref name= ":2" />

<div style="text-align: center;">

<math display="block">\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}, </math>

</div>

which results from the approximation of <math>\mathbf{g}_{n}</math> in (4.2) by its order-''p'' truncated [[Taylor series|Taylor expansion]].

* ''Multistep-type exponential propagation discretization''

Known HOLL discretizations are the following:

* ''Locally Linearized Runge Kutta discretization''

<div style="text-align: center;">

<math display="inline">\mathbf{z}_{n+1}=\mathbf{z}_{n}_n+\mathbf{\phi }(t_t_n,\mathbf{nz}_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(

_{n};h_{n})+h_{n}\sum_{j=1}^{s}b_{j}\mathbf{k}_{j},\quad with \quad \mathbf{k}

\begin{array}{c}

_{i}=\mathbf{q}(t_{n},\mathbf{z}_{n};\text{ }t_{n}+c_{i}h_{n}\mathbf{,}\mathbf{

-\theta \\

}h_{n}\sum_{j=1}^{i-1}a_{ij}\mathbf{k}_{j}), </math>

j

\end{array}

\right) d\theta , </math>

</div>

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

which is obtained by solving (5) via a s-stage [[Runge–Kutta methods|RK scheme]] with coefficients <math>\mathbf{c}=\left[ c_{i}\right] , \mathbf{A}=\left[ a_{ij}\right] \quad and \quad \mathbf{b}=\left[ b_{j}\right]</math>

* ''Local Linear Taylor discretization''

* ''Runge Kutta type Exponential Propagation discretization'' <ref name=":17">Tokman M. (2006). "Efficient integration of large stiff systems of ODEs with exponential propagation iterative (EPI) methods". J. Comput. Physics. 213 (2): 748–776. [[doi:10.1016/j.jcp.2005.08.032]].</ref>

<div style="text-align: center;">

<small><math display="blockinline">\mathbf{z}_{n+1}=\mathbf{z}__n+\mathbf{n\phi}(t_n,\mathbf{z}_n ;h) + h\sum_{j=0}^{p-1} \gamma _{j,p}\nabla^j \mathbf{g}_n (t_n,\phimathbf{z}_n),\quad \text{ with } \quad \gamma_{j,p} = \int\limits_0^1 e^{(t_1-\theta )h\mathbf{nf}_{\mathbf{x}} (t_n,\mathbf{z}_n) } \left(

\begin{array}{c}

_{n};h_{n})+\int_{0}^{h_{n}}e^{\left( h_{n}-s\right) \mathbf{f}_{\mathbf{x}

\theta p\\

}\left( t_{n},\mathbf{z}_{n}\right) }\sum_{j=2}^{p}\frac{\mathbf{c}_{n,j}}{j!

j

}s^{j}ds,\text{ with }\mathbf{c}_{n,j}=\left( \frac{d^{j+1}\mathbf{x}\left(

\end{array}

t\right) }{dt^{j+1}}-\mathbf{f}_{\mathbf{x}}\left( t_{n},\mathbf{z}

\right) d\theta , </math>

_{n}\right) \frac{d^{j}\mathbf{x}\left( t\right) }{dt^{j}}\right) \mid _{t=

\mathbf{z}_{n}}, </math></small>

</div>

which results from the approximationinterpolation of <math>\mathbf{g}_{n}</math> in (4.2) by itsa polynomial of degree order-''p'' truncatedon [[Taylor<math>t_{n},\ldots, series|Taylor expansion]].t_{n}+(p-1)h/p</math>,

* ''[[Exponential integrator|Exponential Rosembrock discretization]]'' (poner link) is obtained by approximating the integral (2) by [[Numerical integration|aquadrature rule]].

* ''Linealized exponential Adams discretization''<ref name=":7">M. Hochbruck.; A. Ostermann. (2011). "Exponential multistep methods of Adams-type". BIT Numer. Math. 51 (4): 889–908. [[doi:10.1007/s10543-011-0332-6]].</ref>

* ''Linealized Exponential Adams discretization''

<div style="text-align: center;">

<small><math display="inline">\mathbf{z}_{n+1}=\mathbf{z}_{n}_n+\mathbf{\phi }(t_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^1e^{(1-\theta )h\mathbf{f}_{\mathbf{x}} \left( t_n,\mathbf{z}_n\right) }\theta \left(

_{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 , </math></small>

</div>

which results from the interpolation of <math>\mathbf{g}_{n}</math> in (4.2) by a [[Hermite polynomials|Hermite polynomial]] of degree ''p'', whereon <math>\nabla ^t_{ln},\mathbfldots, t_{gn-p+1}</math>.

_{n}(t_{m},\mathbf{z}_{m})</math> denotes the ''l''-th backward difference of <math>\mathbf{g}_{n}(t_{m},\mathbf{z}_{m})</math>.

=== Local Linearizationlinearization schemes ===

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

<div style="text-align: center;">

<math>\phi _{j}phi_j(\mathbf{A},h)=\int\limits_{0}limits_0^{h} e^{(h-s)\mathbf{A}} s^{j-1} \, ds,\qquad j=1,2...\ldots, </math>

</div>

where '''A''' is ana ''d <math>\times</math> ''&nbsp;×&nbsp;''d'' matrix. Every numerical implementation <math>\mathbf{y}_{n}_n</math> of a the LocalLL Linear(or of a discretizationHOLL) <math>\mathbf{z}_{n}</math> of any order is generically called ''Local Linearization scheme''.<ref name= ":8"/><ref name=":25">Jimenez, J. C., & Carbonell, F. (2005). "Rate of convergence of local linearization schemes for initial-value problems". Appl. Math. Comput., 171(2), 1282-1295. [[doi:10.1016/j.amc.2005.01.118]].</ref>

==== Computing integrals involving matrix exponential ====

Among a number of algorithms to compute the integrals <math>\phi _{j}</math>, those based on rational Padé and Krylov subspaces approximations for exponential matrix are preferred. For this, a central role is playing by the expression<ref>Carbonell F.; Jimenez J.C.; Pedroso L.M. (2008). "Computing multiple integrals involving matrix exponentials". J. Comput. Appl. Math. 213: 300–305. [https://doi.org/10.1016%2Fj.cam.2007.01.007 doi:10.1016/j.cam.2007.01.007].</ref><ref name=":2">de la Cruz H.; Biscay R.J.; Carbonell F.; Ozaki T.; Jimenez J.C. (2007). "A higher order Local Linearization method for solving ordinary differential equations". Appl. Math. Comput. 185: 197–212. [[doi:10.1016/j.amc.2006.06.096]].</ref><ref name=":13">Jimenez J.C.; Pedroso L.; Carbonell F.; Hernandez V. (2006). "Local linearization method for numerical integration of delay differential equations". SIAM J. Numer. Analysis. 44 (6): 2584–2609. [[doi:10.1137/040607356]].</ref>

<div style="text-align: center;">

<math>\sum\nolimits_{i=1}^{l}\phi _{i}\phi_i(\mathbf{A},h)\mathbf{a}_{i}_i = \mathbf{L}e^{h\mathbf{H}}\mathbf{r}, </math>

\mathbf{H}}\mathbf{r,} </math>

</div>

where <math>\mathbf{a}_{i}_i</math> are ''d''-dimensional vectors,

<div style="text-align: center;">

Line 166 ⟶ 163:

\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+lrl)\times (d+lrl)},

</math>

</div>

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

<div style="text-align: center;">

<math>\mathbf{L}=[\mathbf{I}\mathbf{0}_{d\times l)}],\mathbf{r}=[\mathbf{0}

_{1\times (d+l-1)}\quad1]^{\intercal }\quad, and \quad \mathbf{v}_{i}=\mathbf{a}

_{i}(i-1)!

</math>

</div>

If <math>\mathbf{P}_{p,q}(2^{-k}\mathbf{H}h)

</math> denotes the (''p''; &nbsp;''q'')-[[Padé approximant|Padé approximation]] of <math>e^{2^{-k}\mathbf{H}h}

</math> and ''k'' is the smallest integernatural number such that <math>|2^{-k}\mathbf{H}h|\leq \frac{1}{2}, then</math> <ref name=":12">Jimenez J.C.; de la Cruz H. (2012). "Convergence rate of strong Local Linearization schemes for stochastic differential equations with additive noise". BIT Numer. Math. 52 (2): 357–382. [[doi:10.1007/s10543-011-0360-2]].</ref><ref name=":25" />

</math>

<div style="text-align: center;">

<math>\left\vert \sum\nolimits_{i=1}^{l} \phiphi_i _{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}.</math>

</div>

If <math>\mathbf{\mathbf{k}}_{m,k}^{p,q}(h,\mathbf{H},\mathbf{r})</math> denotes the ''(m; p; q; k)'' [[Krylov subspace|Krylov-Padé approximation]] of <math>e^{h\mathbf{H}}\mathbf{r}</math>, then <ref name=":12"/>

<div style="text-align: center;">

<math>\left\vert \sum\nolimits_{i=1}^{l}\phi _{i}(\mathbf{A},h)\mathbf{a}_i - \mathbf{L\mathbf{k}}_{im,k}-^{p,q}(h,\mathbf{H},\mathbf{r})\right\vert \varpropto h^{\min ({m,p+q+1})},</math>

\mathbf{L\mathbf{k}}_{m,k}^{p,q}(h,\mathbf{H},\mathbf{r})\right\vert

\varpropto h^{\min {m,p+q+1}}.</math>

</div>

where <math>m \leq d</math> is the dimension of the Krylov subspace.

==== Order 2 LL schemes ====

==== Order-2 LL schemes ====

<div style="text-align: center;">

<math>\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,} </math> <ref name=":9">Jimenez J.C.; Biscay R.; Mora C.; Rodriguez L.M. (2002). "Dynamic properties of the Local Linearization method for initial-value problems". Appl. Math. Comput. 126: 63–68. [[doi:10.1016/S0096-3003(00)00100-4]].</ref><ref name=":25"/> <math>\qquad \qquad (4.6)</math>

</div>

where the matrices <math>\mathbf{M}_{n}_n</math>, '''L''' and '''r''' are deffineddefined as

<div style="text-align: center;">

<math>\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) \\

<math>\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

Line 220 ⟶ 208:

<math>\mathbf{L}=\left[

\begin{array}{ll}

\mathbf{I}_{d} & \mathbf{0}_{d\times 2}

\end{array}

\right]</math> and <math>\mathbf{r}^{\intercal }=\left[

\begin{array}{ll}

\mathbf{0}_{1\times (d+1)} & 1

\end{array}

\right] </math> with <math>p+q>1</math> . For large systems of ODEs <ref name=":22" />

<div style="text-align: center;">

<math>\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}>12. </math>

</div>

==== Order -3 LL-Taylor schemes ====

<div style="text-align: center;">

<math>\mathbf{y}_{n+1}=\mathbf{y}_{n}_n+\mathbf{L}_{1}(\mathbf{P}_{p,q}(2^{-k_{n}k_n}

\mathbf{T}_{n}h_{n}_n h_n))^{2^{k_{n}k_n}}\mathbf{r}_{1}_1 \mathbf{,}</math> <ref name=":2" /> <math> \qquad \qquad (4.7)</math>

</div>

where for [[Autonomous system (mathematics)|autonomous]] ODEs the matrices <math>\mathbf{T}_{n}, \mathbf{L}_{1}</math> and <math>\mathbf{r}_{1}</math> are deffineddefined as

<div style="text-align: center;">

Line 250 ⟶ 238:

0 & 0 & 0 & 0 \\

0 & 0 & 0 & 1 \\

0 & 0 & 0 & 0%

\end{array}%

\right] \in \mathbb{R}^{(d+3)\times (d+3)}, </math>

</div>

Line 257 ⟶ 245:

<math>\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[

Line 263 ⟶ 251:

\mathbf{0}_{1\times (d+2)} & 1

\end{array}

\right]</math>. Here, <math>\mathbf{f}_{\mathbf{xx}}</math> denotes the second derivative of '''f''' with respect to '''x''', and ''p + q > 2''. For large systems of ODEs

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

<div style="text-align: center;">

<math>\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}>23. </math>

</div>

==== Order -4 LL-RK schemes ====

<div style="text-align: center;">

<math>\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</math> <ref name=":3" /><ref name=":1" /> <math>\qquad \qquad (4.8)</math>

</div>

where

<div style="text-align: center;">

Line 286 ⟶ 272:

</div>

and

<div style="text-align: center;">

<math>\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)

Line 300 ⟶ 286:

0 & \frac{1}{2} & \frac{1}{2} & 1

\end{array}

\right] , </math> and ''p + q > 3''. For large systems of ODEs, the vector <math>\mathbf{u}_{j} </math> in the above scheme is replaced by <math>\mathbf{u}_{j}=\mathbf{L\mathbf{k}

}_{m_{j},k_{j}}^{p,q}(c_{j}h_{n},\mathbf{M}_{n},\mathbf{r})</math> with <math> m_j > 4.</math>

==== Locally linearized Runge–Kutta scheme of Dormand and Prince ====

scheme is replaced by <math>\mathbf{u}_{j}=\mathbf{L\mathbf{k}

}_{m_{j},k_{j}}^{p,q}(c_{j}h_{n},\mathbf{M}_{n},\mathbf{r}).</math>

==== Locally Linearized Runge-Kutta of Dormand & Prince ====

<div style="text-align: center;">

<math>\mathbf{y}_{n+1}=\mathbf{y}__n+\mathbf{nu}_s+h_n\sum_{j=1}^s b_j \mathbf{uk}_j \qquad \text{ and } \qquad \widehat{\mathbf{y}}_{sn+1}=\mathbf{y}_n+h_\mathbf{nu}_s+h_n\sum_{j=1}^{s \widehat{b}b__j \mathbf{jk}_j,\quad </math>

<ref name=":15">Jimenez J.C.; Sotolongo A.; Sanchez-Bornot J.M. (2014). "Locally Linearized Runge Kutta method of Dormand and Prince". Appl. Math. Comput. 247: 589–606. [[doi:10.1016/j.amc.2014.09.001]].</ref><ref name=":16">Naranjo-Noda, Jimenez J.C. (2021) "Locally Linearized Runge_Kutta method of Dormand and Prince for large systems of initial value problems." J.Comput. Physics. 426: 109946. [[doi:10.1016/j.jcp.2020.109946]].</ref>

\mathbf{k}_{j}\qquad and \qquad \widehat{\mathbf{y}}_{n+1}=\mathbf{y}

<math>\qquad \qquad (4.9)</math>

_{n}+\mathbf{u}_{s}+h_{n}\sum_{j=1}^{s}\widehat{b}_{j}\mathbf{k}_{j},

\qquad \qquad (9)</math>

</div>

where ''s'' = 67 is the number of the stages,

<div style="text-align: center;">

<math>\mathbf{k}_{j}_j=\mathbf{f(}t_t_n+c_j h_n,\mathbf{ny}_n+c_\mathbf{ju}h__j + h_n \sum_{ni=1}^{s-1} a_{j,i} \mathbf{yk}_{n}+_i)-\mathbf{uf}\left( t_n,

_\mathbf{jy}+h__n\right) -\mathbf{nf}_{\sum_mathbf{i=1x}^{s-1}a_{j\left( t_n,i}\mathbf{ky}__n\right) \mathbf{iu})_j\ -\mathbf{f}_t\left( t_t_n,\mathbf{ny}_n \right) c_j h_n, </math>

\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}, </math>

</div>

with <math>\mathbf{k}_{1}\equiv \mathbf{0}</math>, and <math>a_{j,i}, b_{j}, \widehat{b}_{j}_j \quad and \quad c_{j}c_j</math> are the [[Dormand–Prince method|Runge-KuttaRunge–Kutta coefficients of Dormand and Prince]] and ''p'' + ''q'' > 4. For large systems of ODEs, theThe vector <math>\mathbf{u}_{j}_j</math> in the above scheme is replacedcomputed by <math>\mathbf{u}a Padé or Krylor–Padé approximation for small or large systems of ODE, respectively.

_{j}=\mathbf{L\mathbf{k}}_{m_{j},k_{j}}^{p,q}(c_{j}h_{n},\mathbf{M}_{n},

\mathbf{r}).</math>

====  Stability and dynamics ====

[[File:Figure ODE.jpg|thumb|488x488px|'''Fig. 1''' Phase portrait (dashed line) and approximate phase portrait (solid line) of the nonlinear ODE (4.10)-(4.11) computed by the order-2 LL scheme (4.2), the order-4 classical Rugen-Kutta scheme [[Runge–Kutta methods|''RK''4]], ''and the order-4 LLRK''4 schemes (4.8) with step size h=1/2, and p=q=6.]] 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 <math>p\leq q\leq p+2</math>, the LL schemes (4.6)-(4.9) are [[Stiff equation|''A''-stable]].<ref name=":3" /> With ''q'' = ''p'' + 1 or ''q'' = ''p'' + 2, the LL schemes (4.6)–(4.9) are also [[L-stability|''L''-stable]].<ref name=":3" /> For linear ODEs, the LL schemes (4.6)-(4.9) converge with order ''p'' + ''q''.<ref name=":3" /><ref name=":25" /> In addition, with ''p = q = 6'' and <math>m_{n}</math> = ''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.<ref name=":3" /><ref name=":25" /> This includes [[Stiff equation|stiff]] and highly oscillatory linear equations. Moreover, the LL schemes (4.6)-(4.9) are regular for linear ODEs and inherit the [[Symplectic geometry|symplectic structure]] of [[Hamiltonian mechanics|Hamiltonian]] [[harmonic oscillator]]s.<ref name=":2" /><ref name=":9" /> These LL schemes are also linearization preserving, and display a better reproduction of the [[Stable manifold|stable and unstable manifolds]] around [[hyperbolic equilibrium point]]s and [[Limit cycle|periodic orbits]] that [[Numerical methods for ordinary differential equations|other numerical schemes]] with the same stepsize.<ref name=":2" /><ref name=":9" /> For instance, Figure 1 shows the [[phase portrait]] of the ODEs

: <math>

By construction the LL and HOLL discretizations inherit the stability and dynamics of the linear ODEs, but

\begin{align}

& \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{align}

</math>

with <math>f(u,\lambda) =u(1+u+\lambda u^2)^{-1}</math>, <math>\mu =15</math> and <math>\lambda =57</math>, and its approximation by various schemes. This system has two [[Equilibrium point|stable stationary points]] and one [[Equilibrium point|unstable stationary point]] in the region <math>0\leq x_1,x_2\leq 1</math>.

it is not the case of the LL schemes in general. With<math>p\leq q\leq p+2</math>, the LL schemes (6)-(9) are [[Stiff equation|A-stable]].

With ''q = p + 1 or q = p + 2'', the LL schemes (6)-(9) are also [[L-stability|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 equation|stiff]] and highly oscillatory linear equations. Moreover, the LL

schemes (6)-(9) are regular for linear ODEs and inherit the [[Symplectic geometry|symplectic structure]] of [[Hamiltonian mechanics|Hamiltonian]] [[Harmonic oscillator|harmonic]]

[[Harmonic oscillator|oscillators]]. These LL schemes are also linearization preserving, and display a better reproduction of the [[Stable manifold|stable]]

[[Stable manifold|and unstable manifolds]] around [[Hyperbolic equilibrium point|hyperbolic equilibrium points]] and [[Limit cycle|periodic orbits]] that [[Numerical methods for ordinary differential equations|other numerical schemes]]

with the same stepsize.

= LL methods for DDEs =

== LL methods for DDEs ==

Consider the ''d''-dimensional [[Delay differential equation|Delay Differential Equation]] (DDE)

<div style="text-align: center;">

<math>\frac{d\mathbf{x}\left( t\right) }{dt} = \mathbf{f}(t,\leftmathbf{x}( t) ,\mathbf{x}_t(-\lefttau_1),\ldots ,\mathbf{x}_t(-\tau_m)),\qquad t\in[t_0,T] , \qquad\qquad (5.1) </math>

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)</math>

</div>

with ''m'' constant delays <math>\tautau_i>0</math> and initial condition <math>\mathbf{x}_{it_0}(s)=\mathbf{\varphi}(s)</math> for all <math>s\in[-\tau,0], </math> andwhere initial'''f''' conditionis a differentiable function, <math>\mathbf{x}_t:[-\tau,0] \longrightarrow \mathbb{R}^d</math> is the segment function defined as

_{t_{0}}(s)=\mathbf{\varphi }(s)</math> for all <math>s\in \left[ -\tau ,0\right], </math> where ''f'' is a differentiable

<div style="text-align: center;">

function, <math>\mathbf{x}_{t}:\left[

<math>\mathbf{x}_t(s):=\mathbf{x}(t+s),\text{ } s\in[-\tau,0], </math>

-\tau ,0\right] \longrightarrow

\mathbb{R}

^{d}</math> is the segment function deffined as

<div style="text-align: center;">

<math>\mathbf{x}_{t}(s):=\mathbf{x}(t+s),\text{ }s\in \left[ -\tau ,0\right] </math>

</div>

for all <math>t\in \left[ t_{0}t_0,T\right], \mathbf{\varphi }:\left[ -\tau ,0] \longrightarrow \mathbb{R}^d</math> is a given function, and <math>\tau = \max \left\{ \tau_1,\ldots,\tau_m\right\} .</math>

\right] \longrightarrow

\mathbb{R}

^{d}</math> is a given function, and <math>\tau =\max \left\{ \tau _{1,}...,\tau

_{m}\right\} .</math>

=== Local Linearlinear discretization ===

For a time discretization <math>\left( t\right) _{h}_h</math> , the ''Local Linear discretization'' of the DDE (105.1) at each point <math>t_{n+1} \in \left(t)_h</math> is defined by the recursive expression <ref name=":13" />

t\right) _{h}</math> is deffined by the recursive expression

<div style="text-align: center;">

<math>\mathbf{z}_{n+1}=\mathbf{z}_{n}_n+\Phi (t_t_n,\mathbf{nz}_n,h_n;\widetilde{\mathbf{z}}_{nt_n}^1, \ldots,h_{n};\widetilde{\mathbf{z}}_{t_n}^m),

\qquad \qquad (5.2) </math>

\mathbf{z}}_{t_{n}}^{1},..,\widetilde{\mathbf{z}}_{t_{n}}^{m}),

\qquad \qquad \qquad (11)</math>

</div>

Line 392 ⟶ 349:

<div style="text-align: center;">

<math>\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 </math>

<small><math>\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. </math></small>

</div>

where <math>\widetilde{\mathbf{z}}_{t_{n}t_n}^{i}:\left[ -\tau _{i}tau_i,0\right] \longrightarrow \mathbb{R}^d</math> is the segment function defined as

\longrightarrow

\mathbb{R}

^{d}</math> is the segment function deffined as

<div style="text-align: center;">

<math>\widetilde{\mathbf{z}}_{t_{n}t_n}^{i}(s):=\widetilde{\mathbf{z}}^{i}(t_{n}t_n+s),% \text{ } s\in [-\tau_i,0] ,</math>

\text{ }s\in \left[ -\tau _{i},0\right] ,</math>

</div>

and <math>\widetilde{\mathbf{z}}^i:\left[ t_n-\tau_i,t_n\right] \longrightarrow \mathbb{R}^d</math> is a suitable approximation to <math>\mathbf{x}(t)</math> for all <math>t\in \lbrack t_n-\tau_i,t_n]</math> such that <math>\widetilde{\mathbf{z}}^i(t_n)=\mathbf{z}_n.</math> Here,<div style="text-align: center;">

and <math>\widetilde{\mathbf{z}}^{i}:\left[ t_{n}-\tau _{i},t_{n}\right]

<math>\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)),

\longrightarrow

\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)) </math>

\mathbb{R}

^{d}</math>is a suitable approximation to <math>\mathbf{x}(t)</math> for all <math>t\in \lbrack

t_{n}-\tau _{i},t_{n}]</math> such that <math>\widetilde{\mathbf{z}}^{i}(t_{n})=\mathbf{z}_{n}.</math>

Here,

<div style="text-align: center;">

<math>\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})) </math>

</div>

Line 430 ⟶ 368:

<div style="text-align: center;">

<math>\mathbf{c}_{n}_n=\mathbf{f}_{t}_t(t_t_n,\mathbf{nz}_n,\widetilde{\mathbf{z}}_{nt_n}^1 (-\tau_1),\ldots,\widetilde{\mathbf{z}}_{t_n}^m(-\tau_d))

_\text{t_ and }\mathbf{nd}_n=\mathbf{f(}^t_n,\mathbf{1z}(-_n,\tau widetilde{\mathbf{z}}_{1t_n}^1(-\tau_1),...\ldots,\widetilde{\mathbf{z}}_{t_{n}t_n}^{m}(-\tau _{d}tau_d)) </math>

\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})) </math>

</div>

are constant vectors. <math>\mathbf{f}_{t}, \mathbf{f}_{x} \quad and \quad \mathbf{f}

_{x_{t}(-\tau _{i})}</math> denote, respectively, the partial derivatives of '''''f''''' with respect to the variables '''''t''''', and '''''x''''', and <math>\mathbf{x}_{t}(-\tau _{i})</math>. The Local Linear discretization (115.2) converges to the solution of (105.1) with order <math>\alpha =\min \{2,r\},</math> if <math>\widetilde{\mathbf{z}}_{t_{n}}^{i}</math> approximates <math>\mathbf{z}_{t_{n}}^{i}</math> with order <math>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}</math> for all <math>u\in \lbrack 0,h_{n}])</math>.

=== Local Linearizationlinearization schemes ===

[[File:Figure DDE.png|thumb|450x450px|'''Fig. 2''' Approximate paths of the [https://doi.org/10.1016/S0022-5193(05)80142-0 Marchuk et al. (1991)] antiviral immune model described by a stiff system of ten-dimensional nonlinear DDEs with five time delays: top, [[doi:10.1016/S0168-9274(00)00055-6|continuous Runge–Kutta (2,3)]] scheme; bottom, LL scheme (5.3). Step-size ''h''&nbsp;=&nbsp;0.01 fixed, and ''p''&nbsp;=&nbsp;''q''&nbsp;=&nbsp;6.]]

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

Depending on the approximations <math>\widetilde{\mathbf{z}}_{t_{n}}^{i}</math> and on the algorithm to compute <math>\mathbf{\phi }</math> different Local Linearizations schemes can be deÖneddefined. Every numerical implementation <math>\mathbf{y}_{n}</math> of a Local Linear discretization <math>\mathbf{z}_{n}</math> is generically called ''Locallocal Linearizationlinearization scheme''.

==== Order-2 polynomial LL schemes ====

==== Order 2 Polynomial LL schemes ====

<div style="text-align: center;">

<math>\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</math><ref name=":13" /> <math> \qquad (5.3) </math>

\mathbf{M}_{n}h_{n}))^{2^{k_{n}}}\mathbf{r,} \qquad \qquad (12)</math>

</div>

where the matrices <math>\mathbf{M}_{n}, \mathbf{L}\quad</math> and \quad <math>\mathbf{r}</math> are deffineddefined as

<div style="text-align: center;">

Line 460 ⟶ 397:

\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 & 0 & 1 \\

0_{1\times d}0 & 0 & 0

\end{bmatrix}

\in \mathbb{R}^{(d+2)\times (d+2)},</math>

</div>

<math>\mathbf{L}=\left[

\begin{array}{ll}

\mathbf{I}_{d} & \mathbf{0}_{d\times 2}%

\end{array}%

\right]</math> and <math>\mathbf{r}^{\intercal }=\left[

\begin{array}{ll}

\mathbf{0}_{1\times (d+1)} & 1

\end{array}

\right], h_{n}\leq \tau </math>,\quad and \quad <math>p+q>1.</math>. Here, the matrices <math>\mathbf{A}_{n}</math>, <math>\quadmathbf{B}_{n}^{i}</math>, <math>\mathbf{c}_{n}</math> and <math>\mathbf{d}_{n}</math> are defined as in (5.2), but replacing <math>\mathbf{z}</math> by <math>\mathbf{y}</math> and <math>\mathbf{\alpha }_{n}^{i}=(\mathbf{y}(t_{n+1}-\tau _{i})-\mathbf{y}

(t_{n}-\tau _{i}))/h_{n},</math> where

\mathbf{B}_{n}^{i}, \quad \mathbf{c}_{n}\quad and \quad \mathbf{d}_{n}</math> are defined as in '''(11),''' but replacing <math>\mathbf{z}</math> by <math>\mathbf{y}</math> and <math>\mathbf{\alpha }_{n}^{i}=(\mathbf{y}(t_{n+1}-\tau _{i})-\mathbf{y}

(t_{n}-\tau _{i}))/h_{n},</math> where

<div style="text-align: center;">

<small><math>\mathbf{y}\left( t\right) =\mathbf{y}_{n_{t}}+\mathbf{L}(\mathbf{P}%_{p,q}(2^{-k_{n}}

_{p,q}(2^{-k_{n}}\mathbf{M}_{n_{t}}(t-t_{n_{t}})))^{2^{k_{n}}}\mathbf{r},</math>

\text{ with }n_{t}=\max \{n=0,1,2,...,:t_{n}\leq t\text{ and }t_{n}\in

\left( t\right) _{h}\},</math></small>

</div>

with <math>n_{t}=\max \{n=0,1,2,...,:t_{n}\leq t\text{ and }t_{n}\in \left( t\right) _{h}\}</math>, is the ''Local Linear Approximation'' to the solution of (105.1) deÖneddefined through the LL scheme '''(125.3)''' for all <math>t\in \lbrack

t_{0},t_{n}]</math> and by <math>\mathbf{y}\left( t\right) =\mathbf{\varphi }\left(

t\right)</math> for <math>t\in \left[ t_{0}-\tau ,t_{0}\right]</math>. For large systems of DDEs

For large systems of DDEs

<div style="text-align: center;">

Line 496 ⟶ 428:

_{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}),</math>

</div>

with <math>p+q>1</math> and <math>m_{n}>2</math>. 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.

with <math>p+q>1 \quad and \quad m_{n}>2.</math>

== LL methods for RDEs ==

Consider the ''d-''dimensional Random Differential Equation (RDE)

Consider the ''d-dimensional'' Random Differential Equation (RDE)

<div style="text-align: center;"><math>\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)6.1)

</math>

</div>

with initial condition <math>\mathbf{x}(t_{0}t_0)=\mathbf{x}_{0}_0,</math> where <math>\mathbf{

\xi }</math> is a ''k''-dimensional [[Stochastic process|separable finite continuous stochastic process]], and '''f''' is a differentiable function. Suppose that a [[Realization (probability)|realization]] (path) of <math>\mathbf{\xi }</math> is given.

'''f''' is a differentiable function. Suppose that a realization (path) of <math>\mathbf{\xi }</math> is given.

=== Local Linear discretization ===

For a time discretization <math>\left( t\right) _{h}</math>, the ''Local Linear discretization'' of the RDE (13) at each point <math>t_{n+1}\in \left(

For a time discretization <math>\left( t\right) _{h}</math>, the ''Local Linear discretization'' of the RDE (6.1) at each point <math>t_{n+1}\in \left(

t\right) _{h}</math> is defined by the recursive expression

t\right) _{h}</math> is defined by the recursive expression <ref name=":24" />

<div style="text-align: center;">

<math>\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},</math>

</div>

Line 526 ⟶ 457:

<div style="text-align: center;">

<math>\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 </math>

<math>\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 </math>

</div>

and <math>\widetilde{\mathbf{\xi }}</math> is an approximation to the process <math>\mathbf{

\xi }</math> for all <math>t\in \left[ t_{0},T\right]. </math> Here, <math>\mathbf{f}_{x}</math> and <math>\mathbf{f}_{\xi }</math> denote the partial derivatives of <math>\mathbf{f}</math> with respect to <math>\mathbf{x}</math> and <math>\xi </math>, respectively.

=== Local Linearizationlinearization schemes ===

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

[[File:FigureRDE1.png|thumb|377x377px|'''Fig. 3''' Phase portrait of trajectories of the ''Euler'' and ''LL'' schemes in the integration of the nonlinear RDE (6.2)–(6.3) with step size ''h''&nbsp;=&nbsp;1/32, and ''p''&nbsp;=&nbsp;''q''&nbsp;=&nbsp;6.]]

==== LL schemes ====

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

==== LL schemes ====

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

<div style="text-align: center;">

<math>\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 </math> <ref name=":24" /><ref name=":10">Jimenez J.C.; Carbonell F. (2009). "Rate of convergence of local linearization schemes for random differential equations". BIT Numer. Math. 49 (2): 357–373. [[doi:10.1007/s10543-009-0225-0]].</ref> </div>

where the matrices <math>\mathbf{M}_{n}, \quad \mathbf{L} \quad and \quad \mathbf{r}</math> are defined as

<math>\mathbf{M}_{n}=\left[

\begin{array}{ccc}

Line 553 ⟶ 484:

\mathbf{\xi }(t_{n})\right) \\

0 & 0 & 1 \\

0 & 0 & 0%

\end{array}%

\right]</math>

\right] \in \mathbb{R}^{(d+2)\times (d+2)},

</math>

</div>

<math>\mathbf{L}=\left[

\begin{array}{ll}

\mathbf{I}_{d} & \mathbf{0}_{d\times 2}

\end{array}

\right] </math>, <math>\mathbf{r}^{\intercal }=\left[

Line 567 ⟶ 496:

\mathbf{0}_{1\times (d+1)} & 1

\end{array}

\right]</math>, and '''''p+q>1'''''. For large systems of RDEs,<ref name=":10" />

<div style="text-align: center;">

Line 573 ⟶ 502:

_{m_{n},k_{n}}^{p,q}(h_{n},\mathbf{M}_{n},\mathbf{r}),\quad p+q>1

\quad and \quad m_{n}>2.</math>

</div>

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

Figure 3 presents the phase portrait of the RDE

= Strong LL methods for SDEs =

<math>\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)

</math>

<math>\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)

</math>

and its approximation by two numerical schemes, where <math>w^{H}</math> denotes a [[Fractional Brownian motion|fractional Brownian process]] with [[Hurst exponent]] ''H=0.45''.

== Strong LL methods for SDEs ==

Consider the ''d''-dimensional [[Stochastic differential equation|Stochastic Differential Equation]] (SDE)

<div style="text-align: center;">

<math>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 (147.1)</math>

</div>

with initial condition <math>\mathbf{x}(t_{0})=\mathbf{x}_{0}</math>, where the drift coefficient <math>\mathbf{f}</math> and the diffusion coefficient <math>\mathbf{g}_{i}</math> are differentiable functions, and <math>\mathbf{w=(\mathbf{w}}^{1},\ldots ,\mathbf{w}

^{m}\mathbf{)}</math> is an ''m''-dimensional standard [[Wiener process]].

=== Local Linearlinear discretization ===

For a time discretization <math>\left( t\right) _{h}</math> , the order-<math>\mathbb{\gamma }</math> (=1,1.5) ''Strong Local Linear discretization'' of the solution of the SDE (7.1) is defined by the recursive relation <ref name=":14">Jimenez J.C, Shoji I., Ozaki T. (1999) "Simulación of stochastic differential equation through the local linearization method. A comparative study". J. Statist. Physics. 99: 587-602, [[doi:10.1023/A:1004504506041]].</ref><ref name=":20" />

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

<div style="text-align: center;">

Line 597 ⟶ 540:

\mathbf{z}_{n};h_{n})+\mathbf{\xi }(t_{n},\mathbf{z}_{n};h_{n}),\quad with \quad

\mathbf{z}_{0}=\mathbf{x}_{0},</math>

</div>

Line 606 ⟶ 550:

-u)}(\mathbf{f(}t_{n},\mathbf{z}_{n})+\mathbf{a}^{\mathbb{\gamma }}(t_{n},

\mathbf{z}_{n})u)du </math>

</div>

Line 611 ⟶ 556:

<div style="text-align: center;">

<math>\mathbf{\xi } \left( t_{n}t_n,\mathbf{z}_n;\delta \right) = \sum\limits_{i=1}^m \int\nolimits_{t_n}^{t_n+\delta} e^{\mathbf{f}_{n\mathbf{x};} (t_n,\mathbf{z}_n)(t_n+\delta -u)}\rightmathbf{g}_i(u) \, d\mathbf{w}^i(u). </math>

=\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). </math>

</div>

Line 620 ⟶ 563:

<div style="text-align: center;">

<math>
\mathbf{a}^{\mathbb{\gamma }}(t_{n}t_n,\mathbf{z}_{n}_n)=

\left\{

\begin{matrixarray}{cl}

\mathbf{f}_{t}_t(t_{n}t_n,\mathbf{z}_{n}_n) & \qquadtext{for for} \qquad \mathbb{\gamma }=1 \\

\mathbf{f}_{t}_t(t_{n}t_n,\mathbf{z}_{n}_n)\quad +\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}

\left( \mathbf{I}_{d\times d}\times \mathbf{g}_{j}^{\intercal }

\right.

\left( t_{n}\right) \right) \mathbf{f}_{\mathbf{xx}}(t_{n},

</math>

\mathbf{z}_{n})\mathbf{g}_{j}\left( t_{n}\right) & \quad for \quad \mathbb{\gamma }=1.5,

\end{matrix}

\right.</math>

</div>

<math>\mathbf{f}_{\mathbf{x}}, \mathbf{f}_{t}_t</math> denote the partial derivatives of <math>\mathbf{f}</math> with respect to the variables <math>\mathbf{x}</math> and '''''t''''', respectively, and <math>\mathbf{f}_{\mathbf{xx}}</math> the Hessian matrix of <math>\mathbf{f}</math> with respect to <math>\mathbf{x}</math>. The strong Local Linear discretization <math>\mathbf{z}_{n+1}</math> [[Convergence of random variables|converges]] with order <math>\mathbb{\gamma } \quad</math> (=&nbsp;1,&nbsp;1.5)</math> to the solution of '''(147.1)'''.

=== High-order Orderlocal Local Linearlinear discretizations ===

After the local linearization of the drift term of '''(14)''' at <math>(t_{n},

After the local linearization of the drift term of (7.1) at <math>(t_n, \mathbf{z}_{n}_n)</math>, the equation for the residual <math>\mathbf{r}</math> is given by

<div style="text-align: center;">

<math>d\mathbf{r}\left( t\right) =\mathbf{q}_{\gamma }(t_t_n,\mathbf{nz}_n;t \mathbf{,\mathbf{zr}_}(t)) \, dt + \sum\limits_{ni=1};^m \mathbf{g}_i(t) \, d\mathbf{w}^i(t)\mathbf{,}\qquad \mathbf{r}(t_n) = \mathbf{0} </math>

\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} </math>

</div>

for all <math>t\in \lbrack t_{n}t_n,t_{n+1}]</math>, where

<div style="text-align: center;">

<small><math>\mathbf{q}_{\gamma }(t_t_n,\mathbf{nz}_n;s\mathbf{,\xi })=\mathbf{f}(s,\mathbf{z}_n+\mathbf{\phi}_\gamma (t_n,\mathbf{nz}_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). </math>

\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) . </math></small>

</div>

A ''Highhigh-order Order Locallocal Linearlinear discretization'' of the SDE '''''(147.1)''''' at each point <math>t_{n+1}\in \left( t\right) _{h}_h </math> is then defined by the recursive expression <ref name=":4">de la Cruz H.; Biscay R.J.; Jimenez J.C.; Carbonell F.; Ozaki T. (2010). "High Order Local Linearization methods: an approach for constructing A-stable high order explicit schemes for stochastic differential equations with additive noise". BIT Numer. Math. 50 (3): 509–539. [[doi:10.1007/s10543-010-0272-6]].</ref>

<div style="text-align: center;">

<math>\mathbf{z}_{n+1}=\mathbf{z}_{n}_n+\mathbf{\phi }_{\gamma (t_n,\mathbf{z}_n;h_n)+\widetilde{\mathbf{r}}(t_t_n,\mathbf{nz}_n;h_n),\qquad \text{ with } \qquad \mathbf{z}_0=\mathbf{x}_0, </math>

_{n};h_{n})+\widetilde{\mathbf{r}}(t_{n},\mathbf{z}_{n};h_{n}),\qquad with \qquad

\mathbf{z}_{0}=\mathbf{x}_{0}, </math>

</div>

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

=== Local Linearizationlinearization schemes ===

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

==== Order 1 SLL schemes ====

<div style="text-align: center;">

<math>\mathbf{y}_{n+1}=\mathbf{y}_{n}_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 </math> <ref name=":23" /> <math>\qquad \qquad (7.2)</math>

\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)</math>

</div>

where the matrices <math>\mathbf{M}_{n}</math>, <math>\mathbf{L}</math> and <math>\mathbf{r}</math> are defined as in '''(4.6)''', <math>\Delta \mathbf{w}_{n}^{i}</math> is aan [[Independent and identically distributed random variables|i.i.d.]] zero mean [[Normal distribution|Gaussian random variable]] with variance <math>h_{n}</math>, and '''p''&nbsp;+q>1&nbsp;''q''&nbsp;>&nbsp;1. For large systems of SDEs,<ref name=":23" /> in the above scheme <math>(\mathbf{P}_{p,q}(2^{-k_k_n}\mathbf{M}_{n}h_{n}))^{2^{k_n}}\mathbf{r}</math> is replaced by <math>\mathbf{\mathbf{k}}_{m_{n}, k_n}^{p,q}(h_n,\mathbf{M}_n,\mathbf{r})</math>.

\mathbf{M}_{n}h_{n}))^{2^{k_{n}}}\mathbf{r}</math> is replaced by <math>\mathbf{\mathbf{

k}}_{m_{n},k_{n}}^{p,q}(h_{n},\mathbf{M}_{n}^{\gamma },\mathbf{r})</math>.

==== Order 1.5 SLL schemes ====

<div style="text-align: center;">

<math>\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) , \qquad \qquad (7.3)</math>

<small><math>\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)</math></small>

</div>

where the matrices <math>\mathbf{M}_{n}</math>, <math>\mathbf{L}</math> and <math>\mathbf{r}</math> are defined as

<div style="text-align: center;">

<small><math>\mathbf{M}_{n}_n=

\begin{bmatrix}

\mathbf{f}_{\mathbf{x}}(t_t_n,\mathbf{ny}_n) & \mathbf{f}_{t}(t_n,\mathbf{y}__n)+\frac{n1}{2}\sum\limits_{j=1}^{m}\left( \mathbf{I}\otimes \mathbf{g}_j^\intercal (t_n) &\right) \mathbf{f}_{t \mathbf{xx}}(t_t_n,\mathbf{ny},_n)\mathbf{g}_j(t_n) &

\mathbf{f}(t_{n},\mathbf{y}_n) \\

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)}, </math></small>

</div>

<div style="text-align: center;">

<math>\mathbf{L}=\left[

\begin{array}{ll}

\mathbf{I}_{d} & \mathbf{0}_{d\times 2}

\end{array}

\right] , \mathbf{r}^{\intercal }=\left[

Line 718 ⟶ 642:

\mathbf{0}_{1\times (d+1)} & 1

\end{array}

\right]</math>, <math>\Delta \mathbf{z}_{n}^{i}</math> is a i.i.d. zero mean Gaussian random variable with variance <math>E\left( (\Delta \mathbf{z}_{n}^{i})^{2}\right) =% \frac{1}{3}h_{n}^{3}</math> and covariance <math>E(\Delta \mathbf{w}_{n}^{i}\Delta

\mathbf{z}_{n}^{i})=\frac{1}{32}h_{n}^{32}</math> and covariance''p+q>1 <ref name=":12" />''. For large systems of SDEs,<ref name=":12" /> in the above scheme <math>E(\Delta mathbf{P}_{p,q}(2^{-k_{n}}\mathbf{wM}_{n}h_{n}))^{i2^{k_{n}}}\Deltamathbf{r}</math> is replaced by <math>\mathbf{\mathbf{k}}

_{m_{n},k_{n}}^{p,q}(h_{n},\mathbf{M}_{n},\mathbf{r})</math>.

\mathbf{z}_{n}^{i})=\frac{1}{2}h_{n}^{2}</math> and '''p+q>1.''' For large systems of SDEs, in the above scheme <math>(\mathbf{P}_{p,q}(2^{-k_{n}}\mathbf{M}_{n}h_{n}))^{2^{k_{n}}}\mathbf{r}</math> is replaced by <math>\mathbf{\mathbf{k}}

_{m_{n},k_{n}}^{p,q}(h_{n},\mathbf{M}_{n}^{\gamma },\mathbf{r})</math>.

</div>

==== Order 2 SLL-Taylor schemes ====

<div style="text-align: center;">

<small><math>\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}

Line 733 ⟶ 654:

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) }</math></small>

</div>

Line 739 ⟶ 660:

<math>\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), },</math>\qquad \qquad (177.4)</math>

</div>

where <math>\mathbf{M}_{n}</math>, <math>\mathbf{L}</math>, <math>\mathbf{r}</math> and <math>\Delta \mathbf{w}_{n}^{i}</math> are defined as in the order-1 SLL schemes, and <math>\widetilde{J}_{\alpha } </math> is order- '''2''' approximation to the multiple [[Stratonovich integral|Stratonovish integral]] <math>J_{\alpha }</math>.<ref name=":4" />

==== Order 2 SLL-RK schemes ====

[[File:FigureSSDE4.png|thumb|513x513px|'''Fig. 4, Top''': Evolution of domains in the phase plane of the harmonic oscillator (7.6), with ε=0 and ω=σ=1. Images of the initial unit circle (green) are obtained at three time moments ''T'' by the exact solution (black), and by the schemes ''SLL1'' (blue) and ''Implicit Euler'' (red) with ''h=0.05''. '''Bottom''': Expected value of the energy (solid line) along the solution of the nonlinear oscillator (7.6), with ε=1 and ω=100, and its approximation (circles) computed via [[Monte Carlo method|Monte Carlo]] with ''10000'' simulations of the ''SLL1'' scheme with ''h=1/2'' and ''p=q=6''.]]

For SDEs with a single Wiener noise ''(m=1''''')''' <ref name=":4" />

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

<div style="text-align: center;">

<math>\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) }. \qquadquad \qquad (187.5) </math>

</div>

<math>\quad \quad \quad </math>

where

: <math>\mathbf{k}_{1} =\mathbf{f}(t_{n}+\frac{h_{n}}{2},\mathbf{y}_{n}+\widetilde{

<div style="text-align: center;">

<math>\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},

</math>

</div>

: <math>\mathbf{k}_{2} =\mathbf{f}(t_{n}+\frac{h_{n}}{2},\mathbf{y}_{n}+\widetilde{

<div style="text-align: center;">

<math>\mathbf{k\phi }_}(t_{2n} =,\mathbf{fy}(t__{n}+;\frac{h_{n}}{2},)+\mathbf{y}gamma _{n-}+) -\widetildemathbf{f}

\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},</math>

f}_{t}\left( t_{n},\mathbf{y}_{n}\right) \frac{h_{n}}{2},</math>

</div>

with <math>\gamma _{\pm }=\frac{1}{h_{n}}\mathbf{g}\left( t_{n}\right) \Bigl(

with

\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}} \Bigr) </math>.

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

<div style="text-align: center;">

<math>\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\} .

</math>

</div>

==== Stability and dynamics ====

Here, <math>\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} </math> for for low dimensional SDEs, and <math>\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}) </math> for large systems of SDEs, where <math>\mathbf{M}_{n} </math>, <math>\mathbf{L} </math>, <math>\mathbf{r} </math>, <math>\Delta \mathbf{w}_{n}^{i} </math> and <math>\widetilde{J}_{\alpha } </math> are defined as in the order-2 SLL-Taylor schemes, '''p+q>1''' and <math>m_{n}>2 </math>.

By construction, the strong LL and HOLL discretizations inherit the stability and [[Random dynamical system|dynamics]] of the linear SDEs, but it is not the case of the strong LL schemes in general. LL schemes (7.2)-(7.5) with <math>p\leq q\leq p+2 </math> are ''A''-stable, including stiff and highly oscillatory linear equations.<ref name=":12" /> Moreover, for linear SDEs with [[Pullback attractor|random attractors]], these schemes also have a random attractor that [[Convergence in probability|converges in probability]] to the exact one as the stepsize decreases and preserve the [[ergodicity]] of these equations for any stepsize.<ref name=":4" /><ref name=":12" /> 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.<ref name=":4" /><ref name=":5">de la Cruz H.; Jimenez J.C.; Zubelli J.P. (2017). "Locally Linearized methods for the simulation of stochastic oscillators driven by random forces". BIT Numer. Math. 57: 123–151. [[doi:10.1007/s10543-016-0620-2]].</ref> For nonlinear SDEs with small noise (i.e., (7.1) with <math>\mathbf{g}_{i}(t)\approx 0 </math>), the paths of these SLL schemes are basically the nonrandom paths of the LL scheme (4.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.<ref name=":4" /> For instance, Fig 4 shows the evolution of domains in the phase plane and the energy of the stochastic oscillator

==== Stability and dynamics ====

<math>\begin{array}{ll}

By construction the strong LL and HOLL discretizations inherit the stability and [[Random dynamical system|dynamics]] of the linear ODEs, but it is not the case of the strong LL schemes in general. LL schemes '''(15)-(18)''' with <math>p\leq q\leq p+2 </math> are A-stable, which includes stiff and highly oscillatory linear equations. Moreover, for linear SDEs with [[Pullback attractor|random attractors]], these schemes also have a random attractor that [[Convergece of Random variables#Convergence in probability|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 <math>\mathbf{g}_{i}(t)\approx 0 </math>), 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.

dx(t)=y(t)dt, & x_{1}(0)=0.01 \\

dy(t)=-(\omega ^{2}x(t)+\epsilon x^{4}(t))dt+\sigma dw_{t}, & x_{1}(0)=0.1,

\end{array} \qquad \qquad (7.6)</math>

and their approximations by two numerical schemes.

= Weak LL methods for SDEs =

== Weak LL methods for SDEs ==

Consider the ''d''-dimensional stochastic differential equation

<div style="text-align: center;">

<math>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),\qquad t\in \left[ t_{0},T\right] , \qquad \qquad (198.1)</math>

</div>

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

=== Local Linear discretization ===

For a time discretization <math>\left( t\right) _{h}</math>, the order-<math>\mathbb{\beta }</math> <math>(=1,2)</math> ''Weak Local Linear discretization'' of the solution of the SDE '''(198.1)''' is defined by the recursive relation <ref name=":21"/>

<div style="text-align: center;">

Line 825 ⟶ 743:

<div style="text-align: center;">

<math>\mathbf{b}^{\mathbb{\beta }}(t_{n},\mathbf{z}_{n})=

\begin{cases}{ll}

\mathbf{f}_{t}(t_{n},\mathbf{z}_{n}) & \text{for }\mathbb{\beta }=1 \\

\mathbf{f}_{t}(t_{n},\mathbf{z}_{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{z}_{n})\mathbf{g}_{j}\left( t_{n}\right) & \text{for }\mathbb{\beta }

=2,

\end{cases}

.

</math>

</div>

and <math>\mathbf{\eta }(t_{n},\mathbf{z}_{n};\delta )</math> is a zero mean stochastic process with variance matrix

<div style="text-align: center;">

Line 846 ⟶ 763:

</div>

Here, <math>\mathbf{f}_{\mathbf{x}}</math>, <math>\mathbf{f}_{t}</math> denote the partial derivatives of <math>\mathbf{f}</math> with respect to the variables <math>\mathbf{x}</math> and '''t''', respectively, <math>\mathbf{f}_{\mathbf{xx}}</math> the Hessian matrix of <math>\mathbf{f}</math> with respect to <math>\mathbf{x}</math>, and <math>\mathbf{G}(t)=[\mathbf{g}_{1}(t),...\ldots, \mathbf{g}_{m}(t)]</math>. The weak Local Linear discretization <math>\mathbf{z}_{n+1}</math> [[Convergence of random variables|converges]] with order <math>\mathbb{\beta }</math> (=1,2) to the solution of '''(198.1)'''.

=== Local Linearization schemes ===

Depending on the way of computing <math>\mathbf{\phi }_{\mathbb{\beta }}</math> and <math>\mathbf{\Sigma }</math> different numerical schemes couldcan be obtained. Every numerical implementation <math>\mathbf{y}_{n}</math> of athe Weak Local Linear discretization <math>\mathbf{z}_{n}</math> is generically called ''Weak Local Linearization (WLL) scheme''.

==== Order 1 WLL scheme ====

Line 857 ⟶ 774:

<math>\mathbf{y}_{n+1}=\mathbf{y}_{n}+\mathbf{B}_{14}+(\mathbf{B}_{12}\mathbf{B}

_{11}^{\intercal })^{1/2}\mathbf{\xi }_{n} </math>

<ref name=":11">Jimenez J.C.; Carbonell F. (2015). "Convergence rate of weak Local Linearization schemes for stochastic differential equations with additive noise". J. Comput. Appl. Math. 279: 106–122. [[doi:10.1016/j.cam.2014.10.021]].</ref><ref name=":0">Carbonell F.; Jimenez J.C.; Biscay R.J. (2006). "Weak local linear discretizations for stochastic differential equations: convergence and numerical schemes". J. Comput. Appl. Math. 197: 578–596. [[doi:10.1016/j.cam.2005.11.032]].</ref>

</div>

where, for SDEs with autonomous diffusion coefficients, <math>\mathbf{B}_{11}</math>, <math>\mathbf{B}_{12}</math> and <math>\mathbf{B}_{14}</math> are the submatrices defined by the [[Block matrix|partitioned matrix]] <math>\mathbf{B}=\mathbf{P}_{p,q}(2^{-k_{n}}\mathcal{M}_{n}h_{n}))^{2^{k_{n}}}</math>, with

<div style="text-align: center;">

Line 869 ⟶ 787:

\mathbf{0} & \mathbf{0} \\

\mathbf{0} & \mathbf{0} & 0 & 1 \\

\mathbf{0} & \mathbf{0} & 0 & 0%

\end{array}%

\right] \in \mathbb{R}^{(2d+2)\times (2d+2)}, </math>

</div>

and <math>\{\mathbf{\xi }_{n}\}</math> is a sequence of ''d''-dimensional independent [[Bernoulli distribution|two-points distributed random vectors]] satisfying <math>P(\xi _{n}^{k}=\pm 1)=\frac{1}{2} </math>.

<div style="text-align: center;">

<math>P(\xi _{n}^{k}=\pm 1)=\frac{1}{2}. </math>

</div>

==== Order 2 WLL scheme ====

Line 884 ⟶ 798:

<div style="text-align: center;">

<math>\mathbf{y}_{n+1}=\mathbf{y}_{n}+\mathbf{B}_{16}+(\mathbf{B}_{14}\mathbf{B}

_{11}^{\intercal })^{1/2}\mathbf{\xi }_{n}, </math>

<ref name=":11" /><ref name=":0" />

</div>

where <math>\mathbf{B}_{11}</math>, <math>\mathbf{B}_{14}</math> and <math>\mathbf{B}_{16}</math> are the submatrices defined by the partitioned matrix <math>\mathbf{B}=\mathbf{P} _{p,q}(2^{-k_{n}}\mathcal{M}_{n}h_{n}))^{2^{k_{n}}}</math> with

<div style="text-align: center;">

Line 924 ⟶ 839:

</div>

==== Stability and dynamics ====

[[File:Figure WSDE.png|thumb|429x429px|'''''Fig. 5''''' Approximate mean of the SDE (8.2) computed via Monte Carlo with ''100'' simulations of various schemes with ''h=1/16'' and ''p=q=6''.]] By construction, the weak LL discretizations inherit the stability and [[Random dynamical system|dynamics]] of the linear SDEs, but it is not the case of the weak LL schemes in general. WLL schemes, with <math>p\leq q\leq p+2,</math> preserve the [[Moment (mathematics)|first two moments]] of the linear SDEs, and inherits the mean-square stability or instability that such solution may have.<ref name=":11" /> This includes, for instance, the equations of coupled harmonic oscillators driven by random force, and large systems of stiff linear SDEs that result from the method of lines for linear stochastic partial differential equations. Moreover, these WLL schemes preserve the [[ergodicity]] of the linear equations, and are geometrically ergodic for some classes of nonlinear SDEs.<ref name=":6">Hansen N.R. (2003) "Geometric ergodicity of discre-time approximations to multivariate diffusion". Bernoulli. 9 : 725-743, [[doi:10.3150/bj/1066223276]].</ref> For nonlinear SDEs with small noise (i.e., (8.1) with <math>\mathbf{g}_{i}(t)\approx 0</math>), the solutions of these WLL schemes are basically the nonrandom paths of the LL scheme (4.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 mean of the WLL scheme.<ref name=":11" /> For instance, Fig. 5 shows the approximate mean of the SDE

By construction the weak LL and HOLL discretizations inherit the stability and [[Random dynamical system|dynamics]] of the linear ODEs, but it is not the case of the weak LL schemes in general. WLL schemes, with <math>p\leq q\leq p+2,</math> preserve the [[Moment (mathematics)|first two moments]] of the linear SDEs, and inherits the mean-square stability or instability that such solution may have. This includes, for instance, the equations of coupled harmonic oscillators driven by random force, and large systems of stiff linear SDEs that result from the method of lines for linear stochastic partial differential equations. Moreover these WLL schemes preserve the [[ergodicity]] of the linear equations, and are geometrically ergodic for some classes of nonlinear SDEs. For nonlinear SDEs with small noise (e.g., ('''19''') with <math>\mathbf{g}_{i}(t)\approx 0</math>), the solutions of these WLL 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 mean of the WLL scheme.

= Historical notes =

Below is a time line of the main developments of the LL method.

- D.A. Pope (1963) introduces the LL discretization for ODEs and the LL scheme based on Taylor expantion. [[doi:10.1145/366707.367592|doi :10.1145/366707.367592]]

- T. Ozaki (1985) introduces the LL method for the integration and estimation of SDEs. The term\ "Local Linearization" (LL) is used for first time. [https://doi.org/10.1016/S0169-7161(85)05004-0 doi:10.1016/S0169-7161(85)05004-0]

- R. Biscay et al. (1996) reformulate the strong LL method for SDEs. [[doi:10.1007/BF00052324]]

- I. Shoji \& T. Ozaki (1997) reformulate the weak LL method for SDEs. [[doi:10.1111/1467-9892.00064|doi: 10.1111/1467-9892.00064]]

- M. Hochbruck et al. (1998) introduce the LL scheme for ODEs based on Krylov subspace approximation. [[doi:10.1137/S1064827595295337]]

- J.C. Jimenez (2002) introduces the LL scheme for ODEs and SDEs based on rational Padé approximation. [[doi:10.1016/S0893-9659(02)00041-1]]

- F.M. Carbonell et al. (2005) introduce the LL method for RDEs. [[doi:10.1007/s10543-005-2645-9|doi: 10.1007/s10543-005-2645-9]]

- J.C. Jimenez et al. (2006) introduce the LL method for DDEs. [[doi:10.1137/040607356|doi: 10.1137/040607356]]

- de la Cruz et al. (2006, 2007) introduce the first two classes of HOLL integrators for ODEs: LLT and LLRK [[doi:10.1007/11758501\_22]], [[doi:10.1016/j.amc.2006.06.096]]

- de la Cruz et al. (2010) introduce the strong HOLL method for SDEs. [[doi:10.1007/s10543-010-0272-6]]

= References =

de la Cruz H., Biscay R.J., Carbonell F., Ozaki T., Jimenez J.C., A higher order Local Linearization method for solving ordinary differential equations, Appl. Math. Comput., 185 (2007) 197-212.

http://dx.doi.org/10.1016/j.amc.2006.06.096

de la Cruz H., Biscay R.J., Jimenez J.C. and Carbonell F. Local Linearization - Runge Kutta Methods: a class of A-stable explicit integrators for dynamical systems, Math. Comput. Modelling, 57 (2013) 720-740.

http://dx.doi.org/10.1016/j.mcm.2012.08.011

de la Cruz H., Biscay R.J., Jimenez J.C., Carbonell F., Ozaki T., High Order Local Linearization methods: an approach for constructing A-stable high order explicit schemes for stochastic differential equations with additive noise, BIT Numer. Math., 50 (2010) 509--539.

http://dx.doi.org/10.1007/s10543-010-0272-6

de la Cruz H., Jimenez J.C., Zubelli J.P., Locally Linearized methods for the simulation of stochastic oscillators driven by random forces, BIT Numer. Math., 57 (2017) 123-151.

http://dx.doi.org/10.1007/s10543-016-0620-2

M. Hochbruck, A. Ostermann, Exponential multistep methods of Adams-type, BIT Numer. Math. 51 (2011) 889--908.

http://dx.doi.org/10.1007/s10543-011-0332-6

Jimenez J.C., Local Linearization methods for the numerical integration of ordinary differential equations: An overview. ICTP Technical Report 035, 2009.

http://publications.ictp.it

Jimenez J.C., Biscay R., Mora C.M., Rodriguez L.M., Dynamic properties of the Local Linearization method for initial-value problems, Appl. Math. Comput., 126 (2002) 63-80.

http://dx.doi.org/10.1016/S0096-3003(00)00100-4

Jimenez J.C., Carbonell F., Rate of convergence of local linearization schemes for random differential equations, BIT Numer. Math., 49 (2009) 357--373.

http://dx.doi.org/10.1007/s10543-009-0225-0

Jimenez J.C., Carbonell F., Convergence rate of weak Local Linearization schemes for stochastic differential equations with additive noise, J. Comput. Appl. Math., 279 (2015) 106-122.

http://dx.doi.org/10.1016/j.cam.2014.10.021

<math>dx=-t^{2}x\text{ }dt+\frac{3}{2(t+1)}e^{-t^{3}/3}\text{ }dw_{t},\qquad \qquad x(0)=1, \qquad \quad(8.2)</math>

Jimenez J.C., de la Cruz H., Convergence rate of strong Local Linearization schemes for stochastic differential equations with additive noise, BIT Numer. Math., 52 (2012) 357-382.

computed by various schemes.

http://dx.doi.org/10.1007/s10543-011-0360-2

== Historical notes ==

Jimenez J.C., Pedroso L., Carbonell F., Hernadez V., Local linearization method for numerical integration of delay differential equations, SIAM J. Numer. Analysis, 44 (2006) 2584-2609.

Below is a time line of the main developments of the Local Linearization (LL) method.

* Pope D.A. (1963) introduces the LL discretization for ODEs and the LL scheme based on Taylor expansion.<ref name=":18">Pope, D. A. (1963). "An exponential method of numerical integration of ordinary differential equations". Comm. ACM, 6(8), 491-493. [[doi:10.1145/366707.367592]].</ref>

http://dx.doi.org/10.1137/040607356

* Ozaki T. (1985) introduces the LL method for the integration and estimation of SDEs. The term "Local Linearization" is used for first time.<ref name=":19">Ozaki, T. (1985). "Non-linear time series models and dynamical systems". Handbook of statistics, 5, 25-83. [[doi:10.1016/S0169-7161(85)05004-0]].</ref>

* Biscay R. et al. (1996) reformulate the strong LL method for SDEs.<ref name=":20">Biscay, R., Jimenez, J. C., Riera, J. J., & Valdes, P. A. (1996). "Local linearization method for the numerical solution of stochastic differential equations". Annals Inst. Statis. Math. 48(4), 631-644. [[doi:10.1007/BF00052324]].</ref>

* Shoji I. and Ozaki T. (1997) reformulate the weak LL method for SDEs.<ref name=":21">Shoji, I., & Ozaki, T. (1997). "Comparative study of estimation methods for continuous time stochastic processes". J. Time Series Anal. 18(5), 485-506. [[doi:10.1111/1467-9892.00064]].</ref>

* Hochbruck M. et al. (1998) introduce the LL scheme for ODEs based on Krylov subspace approximation.<ref name=":22">Hochbruck, M., Lubich, C., & Selhofer, H. (1998). "Exponential integrators for large systems of differential equations". SIAM J. Scient. Comput. 19(5), 1552-1574. [[doi:10.1137/S1064827595295337]].</ref>

* Jimenez J.C. (2002) introduces the LL scheme for ODEs and SDEs based on rational Padé approximation.<ref name=":23">Jimenez, J. C. (2002). "A simple algebraic expression to evaluate the local linearization schemes for stochastic differential equations". Appl. Math. Letters, 15(6), 775-780. [[doi:10.1016/S0893-9659(02)00041-1]].</ref>

* Carbonell F.M. et al. (2005) introduce the LL method for RDEs.<ref name=":24">Carbonell, F., Jimenez, J. C., Biscay, R. J., & De La Cruz, H. (2005). "The local linearization method for numerical integration of random differential equations". BIT Num. Math. 45(1), 1-14. [[doi:10.1007/S10543-005-2645-9]].</ref>

* Jimenez J.C. et al. (2006) introduce the LL method for DDEs.<ref name=":13" />

* De la Cruz H. et al. (2006, 2007) and Tokman M. (2006) introduce the two classes of HOLL integrators for ODEs: the integrator-based <ref name=":1" /> and the quadrature-based.<ref name=":17" /><ref name=":2" />

* De la Cruz H. et al. (2010) introduce strong HOLL method for SDEs.<ref name=":4" />

== References ==

Jimenez J.C., Sotolongo A. and Sanchez-Bornot J.M., Locally Linearized Runge Kutta method of Dormand and Prince, Appl. Math. Comput., 247 (2014) 589-606.

{{Reflist}}

[[Category:Numerical analysis]]

http://dx.doi.org/10.1016/j.amc.2014.09.001

[[Category:Numerical integration]]