Local linearization method

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

Background

edit

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

High-order local linearization method

edit

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

Local linearization scheme

edit

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

LL methods for ODEs

edit

Consider the d-dimensional Ordinary Differential Equation (ODE)

 

with initial condition  , where   is a differentiable function.

Let   be a time discretization of the time interval   with maximum stepsize h such that   and  . After the local linearization of the equation (4.1) at the time step   the variation of constants formula yields

 

where

 

results from the linear approximation, and

 

is the residual of the linear approximation. Here,   and   denote the partial derivatives of f with respect to the variables x and t, respectively, and  

Local linear discretization

edit

For a time discretization  , the Local Linear discretization of the ODE (4.1) at each point   is defined by the recursive expression [1][2]

 

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

High-order local linear discretizations

edit

For a time discretization   a high-order local linear (HOLL) discretization of the ODE (4.1) at each point   is defined by the recursive expression [1][4][5][6]

 

where   is an order   (> 2) approximation to the residual r   The HOLL discretization (4.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][4][5][6] 1) (quadrature-based) by approximating the integral representation (4.2) of r; and 2) (integrator-based) by using a numerical integrator for the differential representation of r defined by

 

for all  , where

 


HOLL discretizations are, for instance, the followings:

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

 

which is obtained by solving (4.5) via a s-stage explicit Runge–Kutta (RK) scheme with coefficients  .

  • Local linear Taylor discretization[5]

 

which results from the approximation of   in (4.2) by its order-p truncated Taylor expansion.

  • Multistep-type exponential propagation discretization

 

which results from the interpolation of   in (4.2) by a polynomial of degree p on  , where   denotes the j-th backward difference of  .

  • Runge Kutta type Exponential Propagation discretization [7]

 

which results from the interpolation of   in (4.2) by a polynomial of degree p on  ,

  • Linealized exponential Adams discretization[8]

 

which results from the interpolation of   in (4.2) by a Hermite polynomial of degree p on  .

Local linearization schemes

edit

All numerical implementation   of the LL (or of a HOLL) discretization   involves approximations   to integrals   of the form

 

where A is a d × d matrix. Every numerical implementation   of the LL (or of a HOLL)   of any order is generically called Local Linearization scheme.[1][9]

Computing integrals involving matrix exponential

edit

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

 

where   are d-dimensional vectors,

 

 ,    , being   the d-dimensional identity matrix.

If   denotes the (pq)-Padé approximation of   and k is the smallest natural number such that   [12][9]

 

If   denotes the (m; p; q; k) Krylov-Padé approximation of  , then [12]

 

where   is the dimension of the Krylov subspace.

Order-2 LL schemes

edit

  [13][9]  

where the matrices  , L and r are defined as

 

  and   with   . For large systems of ODEs [3]

 

Order-3 LL-Taylor schemes

edit

  [5]  

where for autonomous ODEs the matrices   and   are defined as

 

 . Here,   denotes the second derivative of f with respect to x, and p + q > 2. For large systems of ODEs

 

Order-4 LL-RK schemes

edit

  [4][6]  

where

 

and

 

with   and p + q > 3. For large systems of ODEs, the vector   in the above scheme is replaced by   with  

Locally linearized Runge–Kutta scheme of Dormand and Prince

edit

  [14][15]  

where s = 7 is the number of stages,

 

with  , and   are the Runge–Kutta coefficients of Dormand and Prince and p + q > 4. The vector   in the above scheme is computed by a Padé or Krylor–Padé approximation for small or large systems of ODE, respectively.

Stability and dynamics

edit
 
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 RK4, and the order-4 LLRK4 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  , the LL schemes (4.6)-(4.9) are A-stable.[4] With q = p + 1 or q = p + 2, the LL schemes (4.6)–(4.9) are also L-stable.[4] For linear ODEs, the LL schemes (4.6)-(4.9) converge with order p + q.[4][9] In addition, with p = q = 6 and   = d, all the above described LL schemes yield to the ″exact computation″ (up to the precision of the floating-point arithmetic) of linear ODEs on the current personal computers.[4][9] This includes stiff and highly oscillatory linear equations. Moreover, the LL schemes (4.6)-(4.9) are regular for linear ODEs and inherit the symplectic structure of Hamiltonian harmonic oscillators.[5][13] These LL schemes are also linearization preserving, and display a better reproduction of the stable and unstable manifolds around hyperbolic equilibrium points and periodic orbits that other numerical schemes with the same stepsize.[5][13] For instance, Figure 1 shows the phase portrait of the ODEs

 

with  ,   and  , and its approximation by various schemes. This system has two stable stationary points and one unstable stationary point in the region  .

LL methods for DDEs

edit

Consider the d-dimensional Delay Differential Equation (DDE)

 

with m constant delays   and initial condition   for all   where f is a differentiable function,   is the segment function defined as

 

for all   is a given function, and  

Local linear discretization

edit

For a time discretization   , the Local Linear discretization of the DDE (5.1) at each point   is defined by the recursive expression [11]

 

where

 

  is the segment function defined as

 

and   is a suitable approximation to   for all   such that   Here,

 

are constant matrices and

 

are constant vectors.   denote, respectively, the partial derivatives of f with respect to the variables t and x, and  . The Local Linear discretization (5.2) converges to the solution of (5.1) with order   if   approximates   with order   for all  .

Local linearization schemes

edit
 
Fig. 2 Approximate paths of the Marchuk et al. (1991) antiviral immune model described by a stiff system of ten-dimensional nonlinear DDEs with five time delays: top, continuous Runge–Kutta (2,3) scheme; bottom, LL scheme (5.3). Step-size h = 0.01 fixed, and p = q = 6.

Depending on the approximations   and on the algorithm to compute   different Local Linearizations schemes can be defined. Every numerical implementation   of a Local Linear discretization   is generically called local linearization scheme.

Order-2 polynomial LL schemes

edit

 [11]  

where the matrices   and   are defined as

 

  and  , and  . Here, the matrices  ,  ,   and   are defined as in (5.2), but replacing   by   and   where

 

with  , is the Local Linear Approximation to the solution of (5.1) defined through the LL scheme (5.3) for all   and by   for  . For large systems of DDEs

 

with   and  . Fig. 2 Illustrates the stability of the LL scheme (5.3) and of that of an explicit scheme of similar order in the integration of a stiff system of DDEs.

LL methods for RDEs

edit

Consider the d-dimensional Random Differential Equation (RDE)

 

with initial condition   where   is a k-dimensional separable finite continuous stochastic process, and f is a differentiable function. Suppose that a realization (path) of   is given.

Local Linear discretization

edit

For a time discretization  , the Local Linear discretization of the RDE (6.1) at each point   is defined by the recursive expression [16]

 

where

 

and   is an approximation to the process   for all   Here,   and   denote the partial derivatives of   with respect to   and  , respectively.

Local linearization schemes

edit
 
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 = 1/32, and p = q = 6.

Depending on the approximations   to the process   and of the algorithm to compute  , different Local Linearizations schemes can be defined. Every numerical implementation   of the local linear discretization   is generically called local linearization scheme.

LL schemes

edit
  [16][17]

where the matrices   are defined as

 

 ,  , and p+q>1. For large systems of RDEs,[17]

 

The convergence rate of both schemes is  , where is   the exponent of the Holder condition of  .

Figure 3 presents the phase portrait of the RDE

 

 

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

Strong LL methods for SDEs

edit

Consider the d-dimensional Stochastic Differential Equation (SDE)

 

with initial condition  , where the drift coefficient   and the diffusion coefficient   are differentiable functions, and   is an m-dimensional standard Wiener process.

Local linear discretization

edit

For a time discretization   , the order-  (=1,1.5) Strong Local Linear discretization of the solution of the SDE (7.1) is defined by the recursive relation [18][19]

 

where

 

and

 

Here,

 

  denote the partial derivatives of   with respect to the variables   and t, respectively, and   the Hessian matrix of   with respect to  . The strong Local Linear discretization   converges with order   (= 1, 1.5) to the solution of (7.1).

High-order local linear discretizations

edit

After the local linearization of the drift term of (7.1) at  , the equation for the residual   is given by

 

for all  , where

 

A high-order local linear discretization of the SDE (7.1) at each point   is then defined by the recursive expression [20]

 

where   is a strong approximation to the residual   of order   higher than 1.5. The strong HOLL discretization   converges with order   to the solution of (7.1).

Local linearization schemes

edit

Depending on the way of computing   ,   and   different numerical schemes can be obtained. Every numerical implementation   of a strong Local Linear discretization   of any order is generically called Strong Local Linearization (SLL) scheme.

Order 1 SLL schemes

edit

  [21]  

where the matrices  ,   and   are defined as in (4.6),   is an i.i.d. zero mean Gaussian random variable with variance  , and p + q > 1. For large systems of SDEs,[21] in the above scheme   is replaced by  .

Order 1.5 SLL schemes

edit

 

where the matrices  ,   and   are defined as

 

 ,   is a i.i.d. zero mean Gaussian random variable with variance   and covariance   and p+q>1 [12]. For large systems of SDEs,[12] in the above scheme   is replaced by  .

Order 2 SLL-Taylor schemes

edit

 

 

where  ,  ,   and   are defined as in the order-1 SLL schemes, and   is order 2 approximation to the multiple Stratonovish integral  .[20]

Order 2 SLL-RK schemes

edit
 
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 with 10000 simulations of the SLL1 scheme with h=1/2 and p=q=6.

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

 

 

where

 
 

with  .

Here,   for low dimensional SDEs, and   for large systems of SDEs, where  ,  ,  ,   and   are defined as in the order-2 SLL-Taylor schemes, p+q>1 and  .

Stability and dynamics

edit

By construction, the strong LL and HOLL discretizations inherit the stability and 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   are A-stable, including stiff and highly oscillatory linear equations.[12] Moreover, for linear SDEs with random attractors, these schemes also have a random attractor that converges in probability to the exact one as the stepsize decreases and preserve the ergodicity of these equations for any stepsize.[20][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.[20][22] For nonlinear SDEs with small noise (i.e., (7.1) with  ), 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.[20] For instance, Fig 4 shows the evolution of domains in the phase plane and the energy of the stochastic oscillator

 

and their approximations by two numerical schemes.

Weak LL methods for SDEs

edit

Consider the d-dimensional stochastic differential equation

 

with initial condition  , where the drift coefficient   and the diffusion coefficient   are differentiable functions, and   is an m-dimensional standard Wiener process.

Local Linear discretization

edit

For a time discretization  , the order-    Weak Local Linear discretization of the solution of the SDE (8.1) is defined by the recursive relation [23]

 

where

 

with

 

and   is a zero mean stochastic process with variance matrix

 

Here,  ,   denote the partial derivatives of   with respect to the variables   and t, respectively,   the Hessian matrix of   with respect to  , and  . The weak Local Linear discretization   converges with order   (=1,2) to the solution of (8.1).

Local Linearization schemes

edit

Depending on the way of computing   and   different numerical schemes can be obtained. Every numerical implementation   of the Weak Local Linear discretization   is generically called Weak Local Linearization (WLL) scheme.

Order 1 WLL scheme

edit

  [24][25]

where, for SDEs with autonomous diffusion coefficients,  ,   and   are the submatrices defined by the partitioned matrix  , with

 

and   is a sequence of d-dimensional independent two-points distributed random vectors satisfying  .

Order 2 WLL scheme

edit

  [24][25]

where  ,   and   are the submatrices defined by the partitioned matrix   with

 

 

and

 

Stability and dynamics

edit
 
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 dynamics of the linear SDEs, but it is not the case of the weak LL schemes in general. WLL schemes, with   preserve the first two moments of the linear SDEs, and inherits the mean-square stability or instability that such solution may have.[24] 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.[26] For nonlinear SDEs with small noise (i.e., (8.1) with  ), 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.[24] For instance, Fig. 5 shows the approximate mean of the SDE

 

computed by various schemes.

Historical notes

edit

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.[2]
  • Ozaki T. (1985) introduces the LL method for the integration and estimation of SDEs. The term "Local Linearization" is used for first time.[27]
  • Biscay R. et al. (1996) reformulate the strong LL method for SDEs.[19]
  • Shoji I. and Ozaki T. (1997) reformulate the weak LL method for SDEs.[23]
  • Hochbruck M. et al. (1998) introduce the LL scheme for ODEs based on Krylov subspace approximation.[3]
  • Jimenez J.C. (2002) introduces the LL scheme for ODEs and SDEs based on rational Padé approximation.[21]
  • Carbonell F.M. et al. (2005) introduce the LL method for RDEs.[16]
  • Jimenez J.C. et al. (2006) introduce the LL method for DDEs.[11]
  • De la Cruz H. et al. (2006, 2007) and Tokman M. (2006) introduce the two classes of HOLL integrators for ODEs: the integrator-based [6] and the quadrature-based.[7][5]
  • De la Cruz H. et al. (2010) introduce strong HOLL method for SDEs.[20]

References

edit
  1. ^ a b c d Jimenez J.C. (2009). "Local Linearization methods for the numerical integration of ordinary differential equations: An overview". ICTP Technical Report. 035: 357–373.
  2. ^ a b 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.
  3. ^ a b c 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.
  4. ^ a b c d e f g h 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.
  5. ^ a b c d e f g h 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.
  6. ^ a b c d e 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.
  7. ^ a b 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.
  8. ^ 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.
  9. ^ a b c d e 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.
  10. ^ Carbonell F.; Jimenez J.C.; Pedroso L.M. (2008). "Computing multiple integrals involving matrix exponentials". J. Comput. Appl. Math. 213: 300–305. doi:10.1016/j.cam.2007.01.007.
  11. ^ a b c d 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.
  12. ^ a b c d e f 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.
  13. ^ a b c 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.
  14. ^ 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.
  15. ^ 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.
  16. ^ a b c 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.
  17. ^ a b 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.
  18. ^ 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.
  19. ^ a b 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.
  20. ^ a b c d e f g 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.
  21. ^ a b c 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.
  22. ^ 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.
  23. ^ a b 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.
  24. ^ a b c d 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.
  25. ^ a b 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.
  26. ^ Hansen N.R. (2003) "Geometric ergodicity of discre-time approximations to multivariate diffusion". Bernoulli. 9 : 725-743, doi:10.3150/bj/1066223276.
  27. ^ 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.