Numerical method for differential equations
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 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 r = x − z {\displaystyle \mathbf {r} =\mathbf {x} -\mathbf {z} }
Local linearization scheme [ edit ] A Local Linearization (LL) scheme 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)
d x ( t ) d t = f ( t , x ( t ) ) , t ∈ [ t 0 , T ] , ( 4.1 ) {\displaystyle {\frac {d\mathbf {x} \left(t\right)}{dt}}=\mathbf {f} \left(t,\mathbf {x} \left(t\right)\right),\qquad t\in \left[t_{0},T\right],\qquad \qquad \qquad \qquad (4.1)}
with initial condition x ( t 0 ) = x 0 {\displaystyle \mathbf {x} (t_{0})=\mathbf {x} _{0}} f {\displaystyle \mathbf {f} }
Let ( t ) h = { t n : n = 0 , . . , N } {\displaystyle \left(t\right)_{h}=\{t_{n}:n=0,..,N\}} [ t 0 , T ] {\displaystyle [t_{0},T]} h such that t n < t n + 1 {\displaystyle t_{n}<t_{n+1}} h n = t n + 1 − t n ≤ h {\displaystyle h_{n}=t_{n+1}-t_{n}\leq h} t n {\displaystyle t_{n}} variation of constants formula yields
x ( t n + h ) = x ( t n ) + ϕ ( t n , x ( t n ) ; h ) + r ( t n , x ( t n ) ; h ) , {\displaystyle \mathbf {x} (t_{n}+h)=\mathbf {x} (t_{n})+\mathbf {\phi } (t_{n},\mathbf {x} (t_{n});h)+\mathbf {r} (t_{n},\mathbf {x} (t_{n});h),}
where
ϕ ( t n , z n ; h ) = ∫ 0 h e f x ( t n , z n ) ( h − s ) ( f ( t n , z n ) + f t ( t n , z n ) s ) d s {\displaystyle \mathbf {\phi } (t_{n},\mathbf {z} _{n};h)=\int \limits _{0}^{h}e^{\mathbf {f} _{\mathbf {x} }(t_{n},\mathbf {z} _{n})(h-s)}(\mathbf {f} (t_{n},\mathbf {z} _{n})+\mathbf {f} _{t}(t_{n},\mathbf {z} _{n})s)\,ds\qquad }
results from the linear approximation, and
r ( t n , z n ; h ) = ∫ 0 h e f x ( t n , z n ) ( h − s ) g n ( s , x ( t n + s ) ) d s , ( 4.2 ) {\displaystyle \mathbf {r} (t_{n},\mathbf {z} _{n};h)=\int \limits _{0}^{h}e^{\mathbf {f} _{\mathbf {x} }(t_{n},\mathbf {z} _{n})(h-s)}\mathbf {g} _{n}(s,\mathbf {x} (t_{n}+s))\,ds,\qquad \qquad \qquad (4.2)}
is the residual of the linear approximation. Here, f x {\displaystyle \mathbf {f} _{\mathbf {x} }} f t {\displaystyle \mathbf {f} _{t}} f with respect to the variables x and t , respectively, and g n ( s , u ) = f ( s , u ) − f x ( t n , z n ) u − f t ( t n , z n ) ( s − t n ) − f ( t n , z n ) + f x ( t n , z n ) z n . {\displaystyle \mathbf {g} _{n}(s,\mathbf {u} )=\mathbf {f} (s,\mathbf {u} )-\mathbf {f} _{\mathbf {x} }(t_{n},\mathbf {z} _{n})\mathbf {u} -\mathbf {f} _{t}(t_{n},\mathbf {z} _{n})(s-t_{n})-\mathbf {f} (t_{n},\mathbf {z} _{n})+\mathbf {f} _{\mathbf {x} }(t_{n},\mathbf {z} _{n})\mathbf {z} _{n}.}
Local linear discretization [ edit ] For a time discretization ( t ) h {\displaystyle \left(t\right)_{h}} Local Linear discretization of the ODE (4.1) at each point t n + 1 ∈ ( t ) h {\displaystyle t_{n+1}\in \left(t\right)_{h}} [1] [2]
z n + 1 = z n + ϕ ( t n , z n ; h n ) , with z 0 = x 0 . ( 4.3 ) {\displaystyle \mathbf {z} _{n+1}=\mathbf {z} _{n}+\mathbf {\phi } (t_{n},\mathbf {z} _{n};h_{n}),\qquad {\text{ with }}\quad \mathbf {z} _{0}=\mathbf {x} _{0}.\qquad \qquad \qquad \qquad (4.3)}
The Local Linear discretization (4.3) converges with order 2 [3]
High-order local linear discretizations [ edit ] For a time discretization ( t ) h , {\displaystyle (t)_{h},} high-order local linear (HOLL) discretization of the ODE (4.1) at each point t n + 1 ∈ ( t ) h {\displaystyle t_{n+1}\in (t)_{h}} [1] [4] [5] [6]
z n + 1 = z n + ϕ ( t n , z n ; h n ) + r ~ ( t n , z n ; h n ) , with z 0 = x 0 , ( 4.4 ) {\displaystyle \mathbf {z} _{n+1}=\mathbf {z} _{n}+\mathbf {\phi } (t_{n},\mathbf {z} _{n};h_{n})+{\widetilde {\mathbf {r} }}(t_{n},\mathbf {z} _{n};h_{n}),\qquad {\text{ with }}\quad \mathbf {z} _{0}=\mathbf {x} _{0},\qquad \qquad \qquad (4.4)}
where r ~ {\displaystyle {\tilde {r}}} α {\displaystyle \alpha } 2 ) approximation to the residual r ( i . e . , | r ( t n , z n ; h ) − r ~ ( t n , z n ; h ) | ∝ h α + 1 ) . {\displaystyle (i.e.,\left\vert \mathbf {r} (t_{n},\mathbf {z} _{n};h)-{\widetilde {\mathbf {r} }}(t_{n},\mathbf {z} _{n};h)\right\vert \propto h^{\alpha +1}).} converges with order α {\displaystyle \alpha }
HOLL discretizations can be derived in two ways:[1] [4] [5] [6] r ; and 2) (integrator-based) by using a numerical integrator for the differential representation of r defined by
d r ( t ) d t = q ( t n , z n ; t , r ( t ) ) , with r ( t n ) = 0 , ( 4.5 ) {\displaystyle {\frac {d\mathbf {r} (t)}{dt}}=\mathbf {q} (t_{n},\mathbf {z} _{n};t\mathbf {,\mathbf {r} } (t)\mathbf {)} ,\qquad {\text{ with }}\qquad \mathbf {r} (t_{n})=\mathbf {0,} \qquad \qquad \qquad (4.5)}
for all t ∈ [ t k , t k + 1 ] {\displaystyle t\in \lbrack t_{k},t_{k+1}]}
q ( t n , z n ; s , ξ ) = f ( s , z n + ϕ ( t n , z n ; s − t n ) + ξ ) − f x ( t n , z n ) ϕ ( t n , z n ; s − t n ) − f t ( t n , z n ) ( s − t n ) − f ( t n , z n ) . {\displaystyle \mathbf {q} (t_{n},\mathbf {z} _{n};s\mathbf {,\xi } )=\mathbf {f} (s,\mathbf {z} _{n}+\mathbf {\phi } \left(t_{n},\mathbf {z} _{n};s-t_{n}\right)+\mathbf {\xi } )-\mathbf {f} _{\mathbf {x} }(t_{n},\mathbf {z} _{n})\mathbf {\phi } (t_{n},\mathbf {z} _{n};s-t_{n})-\mathbf {f} _{t}(t_{n},\mathbf {z} _{n})(s-t_{n})-\mathbf {f} (t_{n},\mathbf {z} _{n}).}
Locally Linearized Runge Kutta discretization [6] [4] z n + 1 = z n + ϕ ( t n , z n ; h n ) + h n ∑ j = 1 s b j k j , with k i = q ( t n , z n ; t n + c i h n , h n ∑ j = 1 i − 1 a i j k j ) , {\displaystyle \qquad \mathbf {z} _{n+1}=\mathbf {z} _{n}+\mathbf {\phi } (t_{n},\mathbf {z} _{n};h_{n})+h_{n}\sum _{j=1}^{s}b_{j}\mathbf {k} _{j},\quad {\text{ with }}\quad \mathbf {k} _{i}=\mathbf {q} (t_{n},\mathbf {z} _{n};{\text{ }}t_{n}+c_{i}h_{n}\mathbf {,} \mathbf {} h_{n}\sum _{j=1}^{i-1}a_{ij}\mathbf {k} _{j}),}
which is obtained by solving (4.5) via a s-stage explicit Runge–Kutta (RK) scheme with coefficients c = [ c i ] , A = [ a i j ] a n d b = [ b j ] {\displaystyle \mathbf {c} =\left[c_{i}\right],\mathbf {A} =\left[a_{ij}\right]\quad and\quad \mathbf {b} =\left[b_{j}\right]}
Local linear Taylor discretization [5]
z n + 1 = z n + ϕ ( t n , z n ; h n ) + ∫ 0 h n e ( h n − s ) f x ( t n , z n ) ∑ j = 2 p c n , j j ! s j d s , with c n , j = ( d j + 1 x ( t ) d t j + 1 − f x ( t n , z n ) d j x ( t ) d t j ) ∣ t = z n , {\displaystyle \mathbf {z} _{n+1}=\mathbf {z} _{n}+\mathbf {\phi } (t_{n},\mathbf {z} _{n};h_{n})+\int _{0}^{h_{n}}e^{(h_{n}-s)\mathbf {f} _{\mathbf {x} }\left(t_{n},\mathbf {z} _{n}\right)}\sum _{j=2}^{p}{\frac {\mathbf {c} _{n,j}}{j!}}s^{j}\,ds,{\text{ with }}\mathbf {c} _{n,j}=\left({\frac {d^{j+1}\mathbf {x} (t)}{dt^{j+1}}}-\mathbf {f} _{\mathbf {x} }(t_{n},\mathbf {z} _{n}){\frac {d^{j}\mathbf {x} (t)}{dt^{j}}}\right)\mid _{t=\mathbf {z} _{n}},} which results from the approximation of g n {\displaystyle \mathbf {g} _{n}} p truncated Taylor expansion .
Multistep-type exponential propagation discretization z n + 1 = z n + ϕ ( t n , z n ; h ) + h ∑ j = 0 p − 1 γ j ∇ j g n ( t n , z n ) , w i t h γ j = ( − 1 ) j ∫ 0 1 e ( 1 − θ ) h f x ( t n , z n ) ( − θ j ) d θ , {\textstyle \mathbf {z} _{n+1}=\mathbf {z} _{n}+\mathbf {\phi } (t_{n},\mathbf {z} _{n};h)+h\sum _{j=0}^{p-1}\gamma _{j}\nabla ^{j}\mathbf {g} _{n}(t_{n},\mathbf {z} _{n}),\quad with\quad \gamma _{j}=(-1)^{j}\int \limits _{0}^{1}e^{(1-\theta )h\mathbf {f} _{\mathbf {x} }(t_{n},\mathbf {z} _{n})}\left({\begin{array}{c}-\theta \\j\end{array}}\right)d\theta ,}
which results from the interpolation of g n {\displaystyle \mathbf {g} _{n}} p on t n , … , t n − p + 1 {\displaystyle t_{n},\ldots ,t_{n-p+1}} ∇ j g n ( t m , z m ) {\displaystyle \nabla ^{j}\mathbf {g} _{n}(t_{m},\mathbf {z} _{m})} j -th backward difference of g n ( t m , z m ) {\displaystyle \mathbf {g} _{n}(t_{m},\mathbf {z} _{m})}
Runge Kutta type Exponential Propagation discretization [7] z n + 1 = z n + ϕ ( t n , z n ; h ) + h ∑ j = 0 p − 1 γ j , p ∇ j g n ( t n , z n ) , with γ j , p = ∫ 0 1 e ( 1 − θ ) h f x ( t n , z n ) ( θ p j ) d θ , {\textstyle \mathbf {z} _{n+1}=\mathbf {z} _{n}+\mathbf {\phi } (t_{n},\mathbf {z} _{n};h)+h\sum _{j=0}^{p-1}\gamma _{j,p}\nabla ^{j}\mathbf {g} _{n}(t_{n},\mathbf {z} _{n}),\quad {\text{ with }}\quad \gamma _{j,p}=\int \limits _{0}^{1}e^{(1-\theta )h\mathbf {f} _{\mathbf {x} }(t_{n},\mathbf {z} _{n})}\left({\begin{array}{c}\theta p\\j\end{array}}\right)d\theta ,}
which results from the interpolation of g n {\displaystyle \mathbf {g} _{n}} p on t n , … , t n + ( p − 1 ) h / p {\displaystyle t_{n},\ldots ,t_{n}+(p-1)h/p}
Linealized exponential Adams discretization [8] z n + 1 = z n + ϕ ( t n , z n ; h ) + h ∑ j = 1 p − 1 ∑ l = 1 j γ j + 1 l ∇ l g n ( t n , z n ) , with γ j + 1 = ( − 1 ) j + 1 ∫ 0 1 e ( 1 − θ ) h f x ( t n , z n ) θ ( − θ j ) d θ , {\textstyle \mathbf {z} _{n+1}=\mathbf {z} _{n}+\mathbf {\phi } (t_{n},\mathbf {z} _{n};h)+h\sum _{j=1}^{p-1}\sum _{l=1}^{j}{\frac {\gamma _{j+1}}{l}}\nabla ^{l}\mathbf {g} _{n}(t_{n},\mathbf {z} _{n}),\quad {\text{ with }}\quad \gamma _{j+1}=(-1)^{j+1}\int \limits _{0}^{1}e^{(1-\theta )h\mathbf {f} _{\mathbf {x} }\left(t_{n},\mathbf {z} _{n}\right)}\theta \left({\begin{array}{c}-\theta \\j\end{array}}\right)d\theta ,}
which results from the interpolation of g n {\displaystyle \mathbf {g} _{n}} Hermite polynomial of degree p on t n , … , t n − p + 1 {\displaystyle t_{n},\ldots ,t_{n-p+1}}
Local linearization schemes [ edit ] All numerical implementation y n {\displaystyle \mathbf {y} _{n}} z n {\displaystyle \mathbf {z} _{n}} ϕ ~ j {\displaystyle {\widetilde {\phi }}_{j}} ϕ j {\displaystyle \phi _{j}}
ϕ j ( A , h ) = ∫ 0 h e ( h − s ) A s j − 1 d s , j = 1 , 2 … , {\displaystyle \phi _{j}(\mathbf {A} ,h)=\int \limits _{0}^{h}e^{(h-s)\mathbf {A} }s^{j-1}\,ds,\qquad j=1,2\ldots ,}
where A is a d × d matrix. Every numerical implementation y n {\displaystyle \mathbf {y} _{n}} z n {\displaystyle \mathbf {z} _{n}} Local Linearization scheme .[1] [9]
Computing integrals involving matrix exponential [ edit ] Among a number of algorithms to compute the integrals ϕ j {\displaystyle \phi _{j}} [10] [5] [11]
∑ i = 1 l ϕ i ( A , h ) a i = L e h H r , {\displaystyle \sum \nolimits _{i=1}^{l}\phi _{i}(\mathbf {A} ,h)\mathbf {a} _{i}=\mathbf {L} e^{h\mathbf {H} }\mathbf {r} ,}
where a i {\displaystyle \mathbf {a} _{i}} d -dimensional vectors,
H = [ A v l v l − 1 ⋯ v 1 0 0 1 ⋯ 0 0 0 0 ⋱ 0 ⋮ ⋮ ⋮ ⋱ 1 0 0 0 ⋯ 0 ] ∈ R ( d + l ) × ( d + l ) , {\displaystyle \mathbf {H} ={\begin{bmatrix}\mathbf {A} &\mathbf {v} _{l}&\mathbf {v} _{l-1}&\cdots &\mathbf {v} _{1}\\\mathbf {0} &\mathbf {0} &1&\cdots &0\\\mathbf {0} &\mathbf {0} &0&\ddots &0\\\vdots &\vdots &\vdots &\ddots &1\\\mathbf {0} &\mathbf {0} &0&\cdots &0\end{bmatrix}}\in \mathbb {R} ^{(d+l)\times (d+l)},}
L = [ I 0 d × l ] {\displaystyle \mathbf {L} =[\mathbf {I} \quad \mathbf {0} _{d\times l}]} r = [ 0 1 × ( d + l − 1 ) 1 ] ⊺ , {\displaystyle \mathbf {r} =[\mathbf {0} _{1\times (d+l-1)}\quad 1]^{\intercal },} v i = a i ( i − 1 ) ! {\displaystyle \mathbf {v} _{i}=\mathbf {a} _{i}(i-1)!} I {\displaystyle \mathbf {I} } d -dimensional identity matrix.
If P p , q ( 2 − k H h ) {\displaystyle \mathbf {P} _{p,q}(2^{-k}\mathbf {H} h)} p ; q )-Padé approximation of e 2 − k H h {\displaystyle e^{2^{-k}\mathbf {H} h}} k is the smallest natural number such that | 2 − k H h | ≤ 1 2 , t h e n {\displaystyle |2^{-k}\mathbf {H} h|\leq {\frac {1}{2}},then} [12] [9]
| ∑ i = 1 l ϕ i ( A , h ) a i − L ( P p , q ( 2 − k H h ) ) 2 k r | ∝ h p + q + 1 . {\displaystyle \left\vert \sum \nolimits _{i=1}^{l}\phi _{i}(\mathbf {A} ,h)\mathbf {a} _{i}-\mathbf {L} \left(\mathbf {\mathbf {P} } _{p,q}(2^{-k}\mathbf {H} h)\right)^{2^{k}}\mathbf {r} \right\vert \varpropto h^{p+q+1}.}
If k m , k p , q ( h , H , r ) {\displaystyle \mathbf {\mathbf {k} } _{m,k}^{p,q}(h,\mathbf {H} ,\mathbf {r} )} (m; p; q; k) Krylov-Padé approximation of e h H r {\displaystyle e^{h\mathbf {H} }\mathbf {r} } [12]
| ∑ i = 1 l ϕ i ( A , h ) a i − L k m , k p , q ( h , H , r ) | ∝ h min ( m , p + q + 1 ) , {\displaystyle \left\vert \sum \nolimits _{i=1}^{l}\phi _{i}(\mathbf {A} ,h)\mathbf {a} _{i}-\mathbf {L\mathbf {k} } _{m,k}^{p,q}(h,\mathbf {H} ,\mathbf {r} )\right\vert \varpropto h^{\min({m,p+q+1})},}
where m ≤ d {\displaystyle m\leq d}
Order-2 LL schemes [ edit ] y n + 1 = y n + L ( P p , q ( 2 − k n M n h n ) ) 2 k n r , {\displaystyle \mathbf {y} _{n+1}=\mathbf {y} _{n}+\mathbf {L} (\mathbf {P} _{p,q}(2^{-k_{n}}\mathbf {M} _{n}h_{n}))^{2^{k_{n}}}\mathbf {r,} } [13] [9] ( 4.6 ) {\displaystyle \qquad \qquad (4.6)}
where the matrices M n {\displaystyle \mathbf {M} _{n}} L and r are defined as
M n = [ f x ( t n , y n ) f t ( t n , y n ) f ( t n , y n ) 0 0 1 0 0 0 ] ∈ R ( d + 2 ) × ( d + 2 ) , {\displaystyle \mathbf {M} _{n}={\begin{bmatrix}\mathbf {f} _{\mathbf {x} }(t_{n},\mathbf {y} _{n})&\mathbf {f} _{t}(t_{n},\mathbf {y} _{n})&\mathbf {f} (t_{n},\mathbf {y} _{n})\\0&0&1\\0&0&0\end{bmatrix}}\in \mathbb {R} ^{(d+2)\times (d+2)},}
L = [ I 0 d × 2 ] {\displaystyle \mathbf {L} =\left[{\begin{array}{ll}\mathbf {I} &\mathbf {0} _{d\times 2}\end{array}}\right]} r ⊺ = [ 0 1 × ( d + 1 ) 1 ] {\displaystyle \mathbf {r} ^{\intercal }=\left[{\begin{array}{ll}\mathbf {0} _{1\times (d+1)}&1\end{array}}\right]} p + q > 1 {\displaystyle p+q>1} [3]
y n + 1 = y n + L k m n , k n p , q ( h n , M n , r ) , with m n > 2. {\displaystyle \mathbf {y} _{n+1}=\mathbf {y} _{n}+\mathbf {L\mathbf {k} } _{m_{n},k_{n}}^{p,q}(h_{n},\mathbf {M} _{n},\mathbf {r} )\mathbf {,} \qquad {\text{ with }}\qquad m_{n}>2.}
Order-3 LL-Taylor schemes [ edit ] y n + 1 = y n + L 1 ( P p , q ( 2 − k n T n h n ) ) 2 k n r 1 , {\displaystyle \mathbf {y} _{n+1}=\mathbf {y} _{n}+\mathbf {L} _{1}(\mathbf {P} _{p,q}(2^{-k_{n}}\mathbf {T} _{n}h_{n}))^{2^{k_{n}}}\mathbf {r} _{1}\mathbf {,} } [5] ( 4.7 ) {\displaystyle \qquad \qquad (4.7)}
where for autonomous ODEs the matrices T n , L 1 {\displaystyle \mathbf {T} _{n},\mathbf {L} _{1}} r 1 {\displaystyle \mathbf {r} _{1}}
T n = [ f x ( y n ) ( I ⊗ f ⊺ ( y n ) ) f x x ( y n ) f ( y n ) 0 f ( y n ) 0 0 0 0 0 0 0 1 0 0 0 0 ] ∈ R ( d + 3 ) × ( d + 3 ) , {\displaystyle \mathbf {T} _{n}=\left[{\begin{array}{cccc}\mathbf {f} _{\mathbf {x} }(\mathbf {y} _{n})&(\mathbf {I} \otimes \mathbf {f} ^{\intercal }(\mathbf {y} _{n}))\mathbf {f} _{\mathbf {xx} }(\mathbf {y} _{n})\mathbf {f} (\mathbf {y} _{n})&\mathbf {0} &\mathbf {f} (\mathbf {y} _{n})\\0&0&0&0\\0&0&0&1\\0&0&0&0\end{array}}\right]\in \mathbb {R} ^{(d+3)\times (d+3)},}
L 1 = [ I 0 d × 3 ] a n d r 1 ⊺ = [ 0 1 × ( d + 2 ) 1 ] {\displaystyle \mathbf {L} _{1}=\left[{\begin{array}{ll}\mathbf {I} &\mathbf {0} _{d\times 3}\end{array}}\right]\quad and\quad \mathbf {r} _{1}^{\intercal }=\left[{\begin{array}{ll}\mathbf {0} _{1\times (d+2)}&1\end{array}}\right]} f x x {\displaystyle \mathbf {f} _{\mathbf {xx} }} f with respect to x , and p + q > 2 . For large systems of ODEs
y n + 1 = y n + L k m n , k n p , q ( h n , T n , r ) , with m n > 3. {\displaystyle \mathbf {y} _{n+1}=\mathbf {y} _{n}+\mathbf {L\mathbf {k} } _{m_{n},k_{n}}^{p,q}(h_{n},\mathbf {T} _{n},\mathbf {r} )\mathbf {,} \qquad {\text{ with }}\qquad m_{n}>3.}
Order-4 LL-RK schemes [ edit ] y n + 1 = y n + u 4 + h n 6 ( 2 k 2 + 2 k 3 + k 4 ) , {\displaystyle \mathbf {y} _{n+1}=\mathbf {y} _{n}+\mathbf {u} _{4}+{\frac {h_{n}}{6}}(2\mathbf {k} _{2}+2\mathbf {k} _{3}+\mathbf {k} _{4}),\quad } [4] [6] ( 4.8 ) {\displaystyle \qquad \qquad (4.8)}
where
u j = L ( P p , q ( 2 − κ j M n c j h n ) ) 2 κ j r {\displaystyle \mathbf {u} _{j}=\mathbf {L} (\mathbf {P} _{p,q}(2^{-\kappa _{j}}\mathbf {M} _{n}c_{j}h_{n}))^{2^{\kappa _{j}}}\mathbf {r} }
and
k j = f ( t n + c j h n , y n + u j + c j h n k j − 1 ) − f ( t n , y n ) − f x ( t n , y n ) u j − f t ( t n , y n ) c j h n , {\displaystyle \mathbf {k} _{j}=\mathbf {f} \left(t_{n}+c_{j}h_{n},\mathbf {y} _{n}+\mathbf {u} _{j}+c_{j}h_{n}\mathbf {k} _{j-1}\right)-\mathbf {f} \left(t_{n},\mathbf {y} _{n}\right)-\mathbf {f} _{\mathbf {x} }\left(t_{n},\mathbf {y} _{n}\right)\mathbf {u} _{j}\ -\mathbf {f} _{t}\left(t_{n},\mathbf {y} _{n}\right)c_{j}h_{n},}
with k 1 ≡ 0 , c = [ 0 1 2 1 2 1 ] , {\displaystyle \mathbf {k} _{1}\equiv \mathbf {0} ,c=\left[{\begin{array}{cccc}0&{\frac {1}{2}}&{\frac {1}{2}}&1\end{array}}\right],} p + q > 3 . For large systems of ODEs, the vector u j {\displaystyle \mathbf {u} _{j}} u j = L k m j , k j p , q ( c j h n , M n , r ) {\displaystyle \mathbf {u} _{j}=\mathbf {L\mathbf {k} } _{m_{j},k_{j}}^{p,q}(c_{j}h_{n},\mathbf {M} _{n},\mathbf {r} )} m j > 4. {\displaystyle m_{j}>4.}
Locally linearized Runge–Kutta scheme of Dormand and Prince [ edit ] y n + 1 = y n + u s + h n ∑ j = 1 s b j k j and y ^ n + 1 = y n + u s + h n ∑ j = 1 s b ^ j k j , {\displaystyle \mathbf {y} _{n+1}=\mathbf {y} _{n}+\mathbf {u} _{s}+h_{n}\sum _{j=1}^{s}b_{j}\mathbf {k} _{j}\qquad {\text{ and }}\qquad {\widehat {\mathbf {y} }}_{n+1}=\mathbf {y} _{n}+\mathbf {u} _{s}+h_{n}\sum _{j=1}^{s}{\widehat {b}}_{j}\mathbf {k} _{j},\quad } [14] [15] ( 4.9 ) {\displaystyle \qquad \qquad (4.9)}
where s = 7 is the number of stages,
k j = f ( t n + c j h n , y n + u j + h n ∑ i = 1 s − 1 a j , i k i ) − f ( t n , y n ) − f x ( t n , y n ) u j − f t ( t n , y n ) c j h n , {\displaystyle \mathbf {k} _{j}=\mathbf {f(} t_{n}+c_{j}h_{n},\mathbf {y} _{n}+\mathbf {u} _{j}+h_{n}\sum _{i=1}^{s-1}a_{j,i}\mathbf {k} _{i})-\mathbf {f} \left(t_{n},\mathbf {y} _{n}\right)-\mathbf {f} _{\mathbf {x} }\left(t_{n},\mathbf {y} _{n}\right)\mathbf {u} _{j}\ -\mathbf {f} _{t}\left(t_{n},\mathbf {y} _{n}\right)c_{j}h_{n},}
with k 1 ≡ 0 {\displaystyle \mathbf {k} _{1}\equiv \mathbf {0} } a j , i , b j , b ^ j a n d c j {\displaystyle a_{j,i},b_{j},{\widehat {b}}_{j}\quad and\quad c_{j}} Runge–Kutta coefficients of Dormand and Prince and p + q > 4. The vector u j {\displaystyle \mathbf {u} _{j}}
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 RK 4and 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 p ≤ q ≤ p + 2 {\displaystyle p\leq q\leq p+2} A -stable[4] q = p + 1 or q = p + 2, the LL schemes (4.6)–(4.9) are also L -stable[4] p + q .[4] [9] p = q = 6 and m n {\displaystyle m_{n}} d , all the above described LL schemes yield to the ″exact computation″ (up to the precision of the floating-point arithmetic ) of linear ODEs on the current personal computers.[4] [9] 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] stable and unstable manifolds around hyperbolic equilibrium points and periodic orbits that other numerical schemes with the same stepsize.[5] [13] phase portrait of the ODEs
d x 1 d t = − 2 x 1 + x 2 + 1 − μ f ( x 1 , λ ) ( 4.10 ) d x 2 d t = x 1 − 2 x 2 + 1 − μ f ( x 2 , λ ) ( 4.11 ) {\displaystyle {\begin{aligned}&{\frac {dx_{1}}{dt}}=-2x_{1}+x_{2}+1-\mu f(x_{1},\lambda )\qquad \qquad (4.10)\\[6pt]&{\frac {dx_{2}}{dt}}=x_{1}-2x_{2}+1-\mu f(x_{2},\lambda )\qquad \qquad \quad (4.11)\end{aligned}}} with f ( u , λ ) = u ( 1 + u + λ u 2 ) − 1 {\displaystyle f(u,\lambda )=u(1+u+\lambda u^{2})^{-1}} μ = 15 {\displaystyle \mu =15} λ = 57 {\displaystyle \lambda =57} stable stationary points and one unstable stationary point in the region 0 ≤ x 1 , x 2 ≤ 1 {\displaystyle 0\leq x_{1},x_{2}\leq 1}
LL methods for DDEs [ edit ] Consider the d -dimensional Delay Differential Equation (DDE)
d x ( t ) d t = f ( t , x ( t ) , x t ( − τ 1 ) , … , x t ( − τ m ) ) , t ∈ [ t 0 , T ] , ( 5.1 ) {\displaystyle {\frac {d\mathbf {x} (t)}{dt}}=\mathbf {f} (t,\mathbf {x} (t),\mathbf {x} _{t}(-\tau _{1}),\ldots ,\mathbf {x} _{t}(-\tau _{m})),\qquad t\in [t_{0},T],\qquad \qquad (5.1)}
with m constant delays τ i > 0 {\displaystyle \tau _{i}>0} x t 0 ( s ) = φ ( s ) {\displaystyle \mathbf {x} _{t_{0}}(s)=\mathbf {\varphi } (s)} s ∈ [ − τ , 0 ] , {\displaystyle s\in [-\tau ,0],} f is a differentiable function, x t : [ − τ , 0 ] ⟶ R d {\displaystyle \mathbf {x} _{t}:[-\tau ,0]\longrightarrow \mathbb {R} ^{d}}
x t ( s ) := x ( t + s ) , s ∈ [ − τ , 0 ] , {\displaystyle \mathbf {x} _{t}(s):=\mathbf {x} (t+s),{\text{ }}s\in [-\tau ,0],}
for all t ∈ [ t 0 , T ] , φ : [ − τ , 0 ] ⟶ R d {\displaystyle t\in [t_{0},T],\mathbf {\varphi } :[-\tau ,0]\longrightarrow \mathbb {R} ^{d}} τ = max { τ 1 , … , τ m } . {\displaystyle \tau =\max \left\{\tau _{1},\ldots ,\tau _{m}\right\}.}
Local linear discretization [ edit ] For a time discretization ( t ) h {\displaystyle (t)_{h}} Local Linear discretization of the DDE (5.1) at each point t n + 1 ∈ ( t ) h {\displaystyle t_{n+1}\in (t)_{h}} [11]
z n + 1 = z n + Φ ( t n , z n , h n ; z ~ t n 1 , … , z ~ t n m ) , ( 5.2 ) {\displaystyle \mathbf {z} _{n+1}=\mathbf {z} _{n}+\Phi (t_{n},\mathbf {z} _{n},h_{n};{\widetilde {\mathbf {z} }}_{t_{n}}^{1},\ldots ,{\widetilde {\mathbf {z} }}_{t_{n}}^{m}),\qquad \qquad (5.2)}
where
Φ ( t n , z n , h n ; z ~ t n 1 , … , z ~ t n m ) = ∫ 0 h n e A n ( h n − u ) [ ∑ i = 1 m B n i ( z ~ t n i ( u − τ i ) − z ~ t n i ( − τ i ) ) + d n ] d u + ∫ 0 h n ∫ 0 u e A n ( h n − u ) c n d r d u {\displaystyle \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}
z ~ t n i : [ − τ i , 0 ] ⟶ R d {\displaystyle {\widetilde {\mathbf {z} }}_{t_{n}}^{i}:\left[-\tau _{i},0\right]\longrightarrow \mathbb {R} ^{d}}
z ~ t n i ( s ) := z ~ i ( t n + s ) , s ∈ [ − τ i , 0 ] , {\displaystyle {\widetilde {\mathbf {z} }}_{t_{n}}^{i}(s):={\widetilde {\mathbf {z} }}^{i}(t_{n}+s),{\text{ }}s\in [-\tau _{i},0],}
and z ~ i : [ t n − τ i , t n ] ⟶ R d {\displaystyle {\widetilde {\mathbf {z} }}^{i}:\left[t_{n}-\tau _{i},t_{n}\right]\longrightarrow \mathbb {R} ^{d}} x ( t ) {\displaystyle \mathbf {x} (t)} t ∈ [ t n − τ i , t n ] {\displaystyle t\in \lbrack t_{n}-\tau _{i},t_{n}]} z ~ i ( t n ) = z n . {\displaystyle {\widetilde {\mathbf {z} }}^{i}(t_{n})=\mathbf {z} _{n}.}
A n = f x ( t n , z n , z ~ t n 1 ( − τ 1 ) , … , z ~ t n m ( − τ d ) ) , B n i = f x t ( − τ i ) ( t n , z n , z ~ t n 1 ( − τ 1 ) , … , z ~ t n m ( − τ d ) ) {\displaystyle \mathbf {A} _{n}=\mathbf {f} _{x}(t_{n},\mathbf {z} _{n},{\widetilde {\mathbf {z} }}_{t_{n}}^{1}(-\tau _{1}),\ldots ,{\widetilde {\mathbf {z} }}_{t_{n}}^{m}(-\tau _{d})),{\text{ }}\mathbf {B} _{n}^{i}=\mathbf {f} _{x_{t}(-\tau _{i})}(t_{n},\mathbf {z} _{n},{\widetilde {\mathbf {z} }}_{t_{n}}^{1}(-\tau _{1}),\ldots ,{\widetilde {\mathbf {z} }}_{t_{n}}^{m}(-\tau _{d}))}
are constant matrices and
c n = f t ( t n , z n , z ~ t n 1 ( − τ 1 ) , … , z ~ t n m ( − τ d ) ) and d n = f ( t n , z n , z ~ t n 1 ( − τ 1 ) , … , z ~ t n m ( − τ d ) ) {\displaystyle \mathbf {c} _{n}=\mathbf {f} _{t}(t_{n},\mathbf {z} _{n},{\widetilde {\mathbf {z} }}_{t_{n}}^{1}(-\tau _{1}),\ldots ,{\widetilde {\mathbf {z} }}_{t_{n}}^{m}(-\tau _{d})){\text{ and }}\mathbf {d} _{n}=\mathbf {f(} t_{n},\mathbf {z} _{n},{\widetilde {\mathbf {z} }}_{t_{n}}^{1}(-\tau _{1}),\ldots ,{\widetilde {\mathbf {z} }}_{t_{n}}^{m}(-\tau _{d}))}
are constant vectors. f t , f x a n d f x t ( − τ i ) {\displaystyle \mathbf {f} _{t},\mathbf {f} _{x}\quad and\quad \mathbf {f} _{x_{t}(-\tau _{i})}} f with respect to the variables t and x x t ( − τ i ) {\displaystyle \mathbf {x} _{t}(-\tau _{i})} α = min { 2 , r } , {\displaystyle \alpha =\min\{2,r\},} z ~ t n i {\displaystyle {\widetilde {\mathbf {z} }}_{t_{n}}^{i}} z t n i {\displaystyle \mathbf {z} _{t_{n}}^{i}} r ( i . e . , | z t n i ( u − τ i ) − z ~ t n i ( u − τ i ) | ∝ h n r {\displaystyle 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}} u ∈ [ 0 , h n ] ) {\displaystyle u\in \lbrack 0,h_{n}])}
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; botom, LL scheme (5.3). Step-size h = 0.01 fixed, and p = q = 6.Depending on the approximations z ~ t n i {\displaystyle {\widetilde {\mathbf {z} }}_{t_{n}}^{i}} ϕ {\displaystyle \mathbf {\phi } } y n {\displaystyle \mathbf {y} _{n}} z n {\displaystyle \mathbf {z} _{n}} local linearization scheme .
Order-2 polynomial LL schemes [ edit ] y n + 1 = y n + L ( P p , q ( 2 − k n M n h n ) ) 2 k n r , {\displaystyle \mathbf {y} _{n+1}=\mathbf {y} _{n}+\mathbf {L} (\mathbf {P} _{p,q}(2^{-k_{n}}\mathbf {M} _{n}h_{n}))^{2^{k_{n}}}\mathbf {r,} \quad } [11] ( 5.3 ) {\displaystyle \qquad (5.3)}
where the matrices M n , L {\displaystyle \mathbf {M} _{n},\mathbf {L} } r {\displaystyle \mathbf {r} }
M n = [ A n c n + ∑ i = 1 m B n i α n i d n 0 0 1 0 0 0 ] ∈ R ( d + 2 ) × ( d + 2 ) , {\displaystyle \mathbf {M} _{n}={\begin{bmatrix}\mathbf {A} _{n}&\mathbf {c} _{n}+\sum \limits _{i=1}^{m}\mathbf {B} _{n}^{i}\mathbf {\alpha } _{n}^{i}&\mathbf {d} _{n}\\0&0&1\\0&0&0\end{bmatrix}}\in \mathbb {R} ^{(d+2)\times (d+2)},}
L = [ I 0 d × 2 ] {\displaystyle \mathbf {L} =\left[{\begin{array}{ll}\mathbf {I} &\mathbf {0} _{d\times 2}\end{array}}\right]} r ⊺ = [ 0 1 × ( d + 1 ) 1 ] , h n ≤ τ {\displaystyle \mathbf {r} ^{\intercal }=\left[{\begin{array}{ll}\mathbf {0} _{1\times (d+1)}&1\end{array}}\right],h_{n}\leq \tau } p + q > 1 {\displaystyle p+q>1} A n {\displaystyle \mathbf {A} _{n}} B n i {\displaystyle \mathbf {B} _{n}^{i}} c n {\displaystyle \mathbf {c} _{n}} d n {\displaystyle \mathbf {d} _{n}} z {\displaystyle \mathbf {z} } y {\displaystyle \mathbf {y} } α n i = ( y ( t n + 1 − τ i ) − y ( t n − τ i ) ) / h n , {\displaystyle \mathbf {\alpha } _{n}^{i}=(\mathbf {y} (t_{n+1}-\tau _{i})-\mathbf {y} (t_{n}-\tau _{i}))/h_{n},}
y ( t ) = y n t + L ( P p , q ( 2 − k n M n t ( t − t n t ) ) ) 2 k n r , {\displaystyle \mathbf {y} \left(t\right)=\mathbf {y} _{n_{t}}+\mathbf {L} (\mathbf {P} _{p,q}(2^{-k_{n}}\mathbf {M} _{n_{t}}(t-t_{n_{t}})))^{2^{k_{n}}}\mathbf {r} ,}
with n t = max { n = 0 , 1 , 2 , . . . , : t n ≤ t and t n ∈ ( t ) h } {\displaystyle n_{t}=\max\{n=0,1,2,...,:t_{n}\leq t{\text{ and }}t_{n}\in \left(t\right)_{h}\}} Local Linear Approximation to the solution of (5.1) defined through the LL scheme (5.3) for all t ∈ [ t 0 , t n ] {\displaystyle t\in \lbrack t_{0},t_{n}]} y ( t ) = φ ( t ) {\displaystyle \mathbf {y} \left(t\right)=\mathbf {\varphi } \left(t\right)} t ∈ [ t 0 − τ , t 0 ] {\displaystyle t\in \left[t_{0}-\tau ,t_{0}\right]}
y n + 1 = y n + L k m n , k n p , q ( h n , M n , r ) a n d y ( t ) = y n t + L k m n t , k n t p , q ( t − t n t , M n t , r ) , {\displaystyle \mathbf {y} _{n+1}=\mathbf {y} _{n}+\mathbf {L\mathbf {k} } _{m_{n},k_{n}}^{p,q}(h_{n},\mathbf {M} _{n},\mathbf {r} )\quad and\quad \mathbf {y} \left(t\right)=\mathbf {y} _{n_{t}}+\mathbf {L\mathbf {k} } _{m_{n_{t}},k_{n_{t}}}^{p,q}(t-t_{n_{t}},\mathbf {M} _{n_{t}},\mathbf {r} ),}
with p + q > 1 {\displaystyle p+q>1} m n > 2 {\displaystyle m_{n}>2}
LL methods for RDEs [ edit ] Consider the d- dimensional Random Differential Equation (RDE)
d x ( t ) d t = f ( x ( t ) , ξ ( t ) ) , t ∈ [ t 0 , T ] , ( 6.1 ) {\displaystyle {\frac {d\mathbf {x} \left(t\right)}{dt}}=\mathbf {f} (\mathbf {x} (t),\mathbf {\xi } (t)),\quad t\in \left[t_{0},T\right],\qquad \qquad \qquad (6.1)} with initial condition x ( t 0 ) = x 0 , {\displaystyle \mathbf {x} (t_{0})=\mathbf {x} _{0},} ξ {\displaystyle \mathbf {\xi } } k -dimensional separable finite continuous stochastic process , and f is a differentiable function. Suppose that a realization (path) of ξ {\displaystyle \mathbf {\xi } }
Local Linear discretization [ edit ] For a time discretization ( t ) h {\displaystyle \left(t\right)_{h}} Local Linear discretization of the RDE (6.1) at each point t n + 1 ∈ ( t ) h {\displaystyle t_{n+1}\in \left(t\right)_{h}} [16]
z n + 1 = z n + ϕ ( t n , z n ; h n ) , with z 0 = x 0 , {\displaystyle \mathbf {z} _{n+1}=\mathbf {z} _{n}+\mathbf {\phi } (t_{n},\mathbf {z} _{n};h_{n}),\qquad {\text{ with }}\qquad \mathbf {z} _{0}=\mathbf {x} _{0},}
where
ϕ ( t n , z n ; h n ) = ∫ 0 h n e f x ( z n , ξ ( t n ) ) ( h n − u ) ( f ( z n , ξ ( t n ) ) + f ξ ( z n , ξ ( t n ) ) ( ξ ~ ( t n + u ) − ξ ~ ( t n ) ) ) d u {\displaystyle \mathbf {\phi } (t_{n},\mathbf {z} _{n};h_{n})=\int \limits _{0}^{h_{n}}e^{\mathbf {f} _{\mathbf {x} }(\mathbf {z} _{n},\mathbf {\xi } (t_{n}))(h_{n}-u)}(\mathbf {f(z} _{n},\mathbf {\xi } (t_{n}))+\mathbf {f} _{\mathbf {\xi } }(\mathbf {z} _{n},\mathbf {\xi } (t_{n}))({\widetilde {\mathbf {\xi } }}(t_{n}+u)-{\widetilde {\mathbf {\xi } }}(t_{n})))\,du}
and ξ ~ {\displaystyle {\widetilde {\mathbf {\xi } }}} ξ {\displaystyle \mathbf {\xi } } t ∈ [ t 0 , T ] . {\displaystyle t\in \left[t_{0},T\right].} f x {\displaystyle \mathbf {f} _{x}} f ξ {\displaystyle \mathbf {f} _{\xi }} f {\displaystyle \mathbf {f} } x {\displaystyle \mathbf {x} } ξ {\displaystyle \xi }
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 ξ ~ {\displaystyle {\widetilde {\mathbf {\xi } }}} ξ {\displaystyle \mathbf {\xi } } ϕ {\displaystyle \mathbf {\phi } } y n {\displaystyle \mathbf {y} _{n}} z n {\displaystyle \mathbf {z} _{n}} local linearization scheme.
LL schemes [ edit ] y n + 1 = y n + L ( P p , q ( 2 − k n M n h n ) ) 2 k n r , {\displaystyle \mathbf {y} _{n+1}=\mathbf {y} _{n}+\mathbf {L} (\mathbf {P} _{p,q}(2^{-k_{n}}\mathbf {M} _{n}h_{n}))^{2^{k_{n}}}\mathbf {r,} \quad } [16] [17] where the matrices M n , L a n d r {\displaystyle \mathbf {M} _{n},\quad \mathbf {L} \quad and\quad \mathbf {r} }
M n = [ f x ( y n , ξ ( t n ) ) f ξ ( y n , ξ ( t n ) ( ξ ( t n + 1 ) − ξ ( t n ) ) / h n f ( y n , ξ ( t n ) ) 0 0 1 0 0 0 ] {\displaystyle \mathbf {M} _{n}=\left[{\begin{array}{ccc}\mathbf {f} _{\mathbf {x} }\left(\mathbf {y} _{n},\mathbf {\xi } (t_{n})\right)&\mathbf {f} _{\mathbf {\xi } }(\mathbf {y} _{n},\mathbf {\xi } (t_{n})(\mathbf {\xi } (t_{n+1})-\mathbf {\xi } (t_{n}))/h_{n}&\mathbf {f} \left(\mathbf {y} _{n},\mathbf {\xi } (t_{n})\right)\\0&0&1\\0&0&0\end{array}}\right]}
L = [ I 0 d × 2 ] {\displaystyle \mathbf {L} =\left[{\begin{array}{ll}\mathbf {I} &\mathbf {0} _{d\times 2}\end{array}}\right]} r ⊺ = [ 0 1 × ( d + 1 ) 1 ] {\displaystyle \mathbf {r} ^{\intercal }=\left[{\begin{array}{ll}\mathbf {0} _{1\times (d+1)}&1\end{array}}\right]} p+q>1 . For large systems of RDEs,[17]
y n + 1 = y n + L k m n , k n p , q ( h n , M n , r ) , p + q > 1 a n d m n > 2. {\displaystyle \mathbf {y} _{n+1}=\mathbf {y} _{n}+\mathbf {L\mathbf {k} } _{m_{n},k_{n}}^{p,q}(h_{n},\mathbf {M} _{n},\mathbf {r} ),\quad p+q>1\quad and\quad m_{n}>2.}
The convergence rate of both schemes is m i n { 2 , 2 γ } {\displaystyle min\{2,2\gamma \}} γ {\displaystyle \gamma } ξ {\displaystyle \mathbf {\xi } }
Figure 3 presents the phase portrait of the RDE
d x 1 d t = − x 2 + ( 1 − x 1 2 − x 2 2 ) x 1 sin ( w H ( t ) ) 2 , x 1 ( 0 ) = 0.8 ( 6.2 ) {\displaystyle {\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)}
d x 2 d t = x 1 + ( 1 − x 1 2 − x 2 2 ) x 2 sin ( w H ( t ) ) 2 , x 2 ( 0 ) = 0.1 , ( 6.3 ) {\displaystyle {\frac {dx_{2}}{dt}}=x_{1}+(1-x_{1}^{2}-x_{2}^{2})x_{2}\sin(w^{H}(t))^{2},\qquad \qquad x_{2}(0)=0.1,\qquad (6.3)}
and its approximation by two numerical schemes, where w H {\displaystyle w^{H}} fractional Brownian process with Hurst exponent H=0.45 .
Strong LL methods for SDEs [ edit ] Consider the d -dimensional Stochastic Differential Equation (SDE)
d x ( t ) = f ( t , x ( t ) ) d t + ∑ i = 1 m g i ( t ) d w i ( t ) , t ∈ [ t 0 , T ] , ( 7.1 ) {\displaystyle d\mathbf {x} (t)=\mathbf {f} (t,\mathbf {x} (t))dt+\sum \limits _{i=1}^{m}\mathbf {g} _{i}(t)d\mathbf {w} ^{i}(t),\quad t\in \left[t_{0},T\right],\qquad \qquad \qquad (7.1)}
with initial condition x ( t 0 ) = x 0 {\displaystyle \mathbf {x} (t_{0})=\mathbf {x} _{0}} f {\displaystyle \mathbf {f} } g i {\displaystyle \mathbf {g} _{i}} w = ( w 1 , … , w m ) {\displaystyle \mathbf {w=(\mathbf {w} } ^{1},\ldots ,\mathbf {w} ^{m}\mathbf {)} } m -dimensional standard Wiener process .
Local linear discretization [ edit ] For a time discretization ( t ) h {\displaystyle \left(t\right)_{h}} γ {\displaystyle \mathbb {\gamma } } Strong Local Linear discretization of the solution of the SDE (7.1) is defined by the recursive relation [18] [19]
z n + 1 = z n + ϕ γ ( t n , z n ; h n ) + ξ ( t n , z n ; h n ) , w i t h z 0 = x 0 , {\displaystyle \mathbf {z} _{n+1}=\mathbf {z} _{n}+\mathbf {\phi } _{\mathbb {\gamma } }(t_{n},\mathbf {z} _{n};h_{n})+\mathbf {\xi } (t_{n},\mathbf {z} _{n};h_{n}),\quad with\quad \mathbf {z} _{0}=\mathbf {x} _{0},}
where
ϕ γ ( t n , z n ; δ ) = ∫ 0 δ e f x ( t n , y n ) ( δ − u ) ( f ( t n , z n ) + a γ ( t n , z n ) u ) d u {\displaystyle \mathbf {\phi } _{\mathbb {\gamma } }(t_{n},\mathbf {z} _{n};\delta )=\int _{0}^{\delta }e^{\mathbf {f} _{\mathbf {x} }(t_{n},\mathbf {y} _{n})(\delta -u)}(\mathbf {f(} t_{n},\mathbf {z} _{n})+\mathbf {a} ^{\mathbb {\gamma } }(t_{n},\mathbf {z} _{n})u)du}
and
ξ ( t n , z n ; δ ) = ∑ i = 1 m ∫ t n t n + δ e f x ( t n , z n ) ( t n + δ − u ) g i ( u ) d w i ( u ) . {\displaystyle \mathbf {\xi } \left(t_{n},\mathbf {z} _{n};\delta \right)=\sum \limits _{i=1}^{m}\int \nolimits _{t_{n}}^{t_{n}+\delta }e^{\mathbf {f} _{\mathbf {x} }(t_{n},\mathbf {z} _{n})(t_{n}+\delta -u)}\mathbf {g} _{i}(u)\,d\mathbf {w} ^{i}(u).}
Here,
a γ ( t n , z n ) = { f t ( t n , z n ) for γ = 1 f t ( t n , z n ) + 1 2 ∑ j = 1 m ( I ⊗ g j ⊺ ( t n ) ) f x x ( t n , z n ) g j ( t n ) for γ = 1.5 , {\displaystyle \mathbf {a} ^{\mathbb {\gamma } }(t_{n},\mathbf {z} _{n})=\left\{{\begin{array}{cl}\mathbf {f} _{t}(t_{n},\mathbf {z} _{n})&{\text{for }}\qquad \mathbb {\gamma } =1\\\mathbf {f} _{t}(t_{n},\mathbf {z} _{n})+{\frac {1}{2}}\sum \limits _{j=1}^{m}(\mathbf {I} \otimes \mathbf {g} _{j}^{\intercal }(t_{n}))\mathbf {f} _{\mathbf {xx} }(t_{n},\mathbf {z} _{n})\mathbf {g} _{j}(t_{n})&{\text{for }}\quad \mathbb {\gamma } =1.5,\end{array}}\right.}
f x , f t {\displaystyle \mathbf {f} _{\mathbf {x} },\mathbf {f} _{t}} f {\displaystyle \mathbf {f} } x {\displaystyle \mathbf {x} } t f x x {\displaystyle \mathbf {f} _{\mathbf {xx} }} f {\displaystyle \mathbf {f} } x {\displaystyle \mathbf {x} } z n + 1 {\displaystyle \mathbf {z} _{n+1}} converges with order γ {\displaystyle \mathbb {\gamma } }
High-order local linear discretizations [ edit ] After the local linearization of the drift term of (7.1) at ( t n , z n ) {\displaystyle (t_{n},\mathbf {z} _{n})} r {\displaystyle \mathbf {r} }
d r ( t ) = q γ ( t n , z n ; t , r ( t ) ) d t + ∑ i = 1 m g i ( t ) d w i ( t ) , r ( t n ) = 0 {\displaystyle d\mathbf {r} (t)=\mathbf {q} _{\gamma }(t_{n},\mathbf {z} _{n};t\mathbf {,\mathbf {r} } (t))\,dt+\sum \limits _{i=1}^{m}\mathbf {g} _{i}(t)\,d\mathbf {w} ^{i}(t)\mathbf {,} \qquad \mathbf {r} (t_{n})=\mathbf {0} }
for all t ∈ [ t n , t n + 1 ] {\displaystyle t\in \lbrack t_{n},t_{n+1}]}
q γ ( t n , z n ; s , ξ ) = f ( s , z n + ϕ γ ( t n , z n ; s − t n ) + ξ ) − f x ( t n , z n ) ϕ γ ( t n , z n ; s − t n ) − a γ ( t n , z n ) ( s − t n ) − f ( t n , z n ) . {\displaystyle \mathbf {q} _{\gamma }(t_{n},\mathbf {z} _{n};s\mathbf {,\xi } )=\mathbf {f} (s,\mathbf {z} _{n}+\mathbf {\phi } _{\gamma }(t_{n},\mathbf {z} _{n};s-t_{n})+\mathbf {\xi } )-\mathbf {f} _{\mathbf {x} }(t_{n},\mathbf {z} _{n})\mathbf {\phi } _{\gamma }(t_{n},\mathbf {z} _{n};s-t_{n})-\mathbf {a} ^{\gamma }(t_{n},\mathbf {z} _{n})(s-t_{n})-\mathbf {f} (t_{n},\mathbf {z} _{n}).}
A high-order local linear discretization of the SDE (7.1) at each point t n + 1 ∈ ( t ) h {\displaystyle t_{n+1}\in (t)_{h}} [20]
z n + 1 = z n + ϕ γ ( t n , z n ; h n ) + r ~ ( t n , z n ; h n ) , with z 0 = x 0 , {\displaystyle \mathbf {z} _{n+1}=\mathbf {z} _{n}+\mathbf {\phi } _{\gamma }(t_{n},\mathbf {z} _{n};h_{n})+{\widetilde {\mathbf {r} }}(t_{n},\mathbf {z} _{n};h_{n}),\qquad {\text{ with }}\qquad \mathbf {z} _{0}=\mathbf {x} _{0},}
where r ~ {\displaystyle {\widetilde {\mathbf {r} }}} r {\displaystyle \mathbf {r} } α {\displaystyle \alpha } 1.5 . The strong HOLL discretization z n + 1 {\displaystyle \mathbf {z} _{n+1}} α {\displaystyle \alpha }
Local linearization schemes [ edit ] Depending on the way of computing ϕ γ {\displaystyle \mathbf {\phi } _{\mathbb {\gamma } }} ξ {\displaystyle \mathbf {\xi } } r ~ {\displaystyle {\widetilde {\mathbf {r} }}} y n {\displaystyle \mathbf {y} _{n}} z n {\displaystyle \mathbf {z} _{n}} Strong Local Linearization (SLL) scheme .
Order 1 SLL schemes [ edit ] y n + 1 = y n + L ( P p , q ( 2 − k n M n h n ) ) 2 k n r + ∑ i = 1 m g i ( t n ) Δ w n i , {\displaystyle \mathbf {y} _{n+1}=\mathbf {y} _{n}+\mathbf {L} (\mathbf {P} _{p,q}(2^{-k_{n}}\mathbf {M} _{n}h_{n}))^{2^{k_{n}}}\mathbf {r+} \sum \limits _{i=1}^{m}\mathbf {g} _{i}(t_{n})\Delta \mathbf {w} _{n}^{i},\quad } [21] ( 7.2 ) {\displaystyle \qquad \qquad (7.2)}
where the matrices M n {\displaystyle \mathbf {M} _{n}} L {\displaystyle \mathbf {L} } r {\displaystyle \mathbf {r} } Δ w n i {\displaystyle \Delta \mathbf {w} _{n}^{i}} i.i.d. zero mean Gaussian random variable with variance h n {\displaystyle h_{n}} p + q > 1. For large systems of SDEs,[21] ( P p , q ( 2 − k n M n h n ) ) 2 k n r {\displaystyle (\mathbf {P} _{p,q}(2^{-k_{n}}\mathbf {M} _{n}h_{n}))^{2^{k_{n}}}\mathbf {r} } k m n , k n p , q ( h n , M n , r ) {\displaystyle \mathbf {\mathbf {k} } _{m_{n},k_{n}}^{p,q}(h_{n},\mathbf {M} _{n},\mathbf {r} )}
Order 1.5 SLL schemes [ edit ] y n + 1 = y n + L ( P p , q ( 2 − k n M n h n ) ) 2 k n r + ∑ i = 1 m ( g i ( t n ) Δ w n i f x ( t n , y ~ n ) g i ( t n ) Δ z n i + d g i ( t n ) d t ( Δ w n i h n − Δ z n i ) ) , ( 7.3 ) {\displaystyle \mathbf {y} _{n+1}=\mathbf {y} _{n}+\mathbf {L} (\mathbf {P} _{p,q}(2^{-k_{n}}\mathbf {M} _{n}h_{n}))^{2^{k_{n}}}\mathbf {r} +\sum \limits _{i=1}^{m}\left(\mathbf {g} _{i}(t_{n})\Delta \mathbf {w} _{n}^{i}\mathbf {f} _{\mathbf {x} }(t_{n},{\widetilde {\mathbf {y} }}_{n})\mathbf {g} _{i}(t_{n})\Delta \mathbf {z} _{n}^{i}+{\frac {d\mathbf {g} _{i}(t_{n})}{dt}}(\Delta \mathbf {w} _{n}^{i}h_{n}-\Delta \mathbf {z} _{n}^{i})\right),\qquad \qquad (7.3)}
From Wikipedia the free encyclopedia
![{\displaystyle {\frac {d\mathbf {x} \left(t\right)}{dt}}=\mathbf {f} \left(t,\mathbf {x} \left(t\right)\right),\qquad t\in \left[t_{0},T\right],\qquad \qquad \qquad \qquad (4.1)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/b689077ae401751f25375eb3338ed66cc57f756a)
![{\displaystyle [t_{0},T]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/986ba7ea2bc36ce31beb5c5f4faffbfb6405f69b)
![{\displaystyle t\in \lbrack t_{k},t_{k+1}]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/dacd0fb31482f9ac95b8c3261d9dc00f1055dca4)
![{\displaystyle \mathbf {c} =\left[c_{i}\right],\mathbf {A} =\left[a_{ij}\right]\quad and\quad \mathbf {b} =\left[b_{j}\right]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/6fbc317004aaf306e62ecc1c6d08e4f49d9d0543)
![{\displaystyle \mathbf {L} =[\mathbf {I} \quad \mathbf {0} _{d\times l}]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/75b8ac2636b6ebd86e0bba8f06886336cfc6ec82)
![{\displaystyle \mathbf {r} =[\mathbf {0} _{1\times (d+l-1)}\quad 1]^{\intercal },}](https://wikimedia.org/api/rest_v1/media/math/render/svg/3710227203d25399ba90fb8e8f94a9d34ff29aa9)
![{\displaystyle \mathbf {L} =\left[{\begin{array}{ll}\mathbf {I} &\mathbf {0} _{d\times 2}\end{array}}\right]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/05d54f18477ecdb4b717abf25a75b96d17729bcd)
![{\displaystyle \mathbf {r} ^{\intercal }=\left[{\begin{array}{ll}\mathbf {0} _{1\times (d+1)}&1\end{array}}\right]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/f82230fb66e66705351377690e9f14c8d4a41224)
![{\displaystyle \mathbf {T} _{n}=\left[{\begin{array}{cccc}\mathbf {f} _{\mathbf {x} }(\mathbf {y} _{n})&(\mathbf {I} \otimes \mathbf {f} ^{\intercal }(\mathbf {y} _{n}))\mathbf {f} _{\mathbf {xx} }(\mathbf {y} _{n})\mathbf {f} (\mathbf {y} _{n})&\mathbf {0} &\mathbf {f} (\mathbf {y} _{n})\\0&0&0&0\\0&0&0&1\\0&0&0&0\end{array}}\right]\in \mathbb {R} ^{(d+3)\times (d+3)},}](https://wikimedia.org/api/rest_v1/media/math/render/svg/04b4535597fcece8be4c40e7b4ae75229d8f9831)
![{\displaystyle \mathbf {L} _{1}=\left[{\begin{array}{ll}\mathbf {I} &\mathbf {0} _{d\times 3}\end{array}}\right]\quad and\quad \mathbf {r} _{1}^{\intercal }=\left[{\begin{array}{ll}\mathbf {0} _{1\times (d+2)}&1\end{array}}\right]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/f495cb0102ab17316c499936bd236cdb35f4ebbd)
![{\displaystyle \mathbf {k} _{1}\equiv \mathbf {0} ,c=\left[{\begin{array}{cccc}0&{\frac {1}{2}}&{\frac {1}{2}}&1\end{array}}\right],}](https://wikimedia.org/api/rest_v1/media/math/render/svg/5c7569810b131c3d8e9ad24b6201f9454d90fbbc)
![{\displaystyle {\begin{aligned}&{\frac {dx_{1}}{dt}}=-2x_{1}+x_{2}+1-\mu f(x_{1},\lambda )\qquad \qquad (4.10)\\[6pt]&{\frac {dx_{2}}{dt}}=x_{1}-2x_{2}+1-\mu f(x_{2},\lambda )\qquad \qquad \quad (4.11)\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/855ec9f888ae9f2a03221800414cbf76ef9edd46)
![{\displaystyle {\frac {d\mathbf {x} (t)}{dt}}=\mathbf {f} (t,\mathbf {x} (t),\mathbf {x} _{t}(-\tau _{1}),\ldots ,\mathbf {x} _{t}(-\tau _{m})),\qquad t\in [t_{0},T],\qquad \qquad (5.1)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/682b3725e57d7e2ee5324a0a139b7e2d1063767b)
![{\displaystyle s\in [-\tau ,0],}](https://wikimedia.org/api/rest_v1/media/math/render/svg/9e9d0232ea2938ec04dafc8665c88403f10670c3)
![{\displaystyle \mathbf {x} _{t}:[-\tau ,0]\longrightarrow \mathbb {R} ^{d}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/026aa256631bc0e92437d8cea3a76cd0c6c2aded)
![{\displaystyle \mathbf {x} _{t}(s):=\mathbf {x} (t+s),{\text{ }}s\in [-\tau ,0],}](https://wikimedia.org/api/rest_v1/media/math/render/svg/857fd392ec867a13d23dcb2ea2d705d132a706d6)
![{\displaystyle t\in [t_{0},T],\mathbf {\varphi } :[-\tau ,0]\longrightarrow \mathbb {R} ^{d}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/332fbe38ea2b9aa1ca97db545170bba5e188883f)
![{\displaystyle \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}](https://wikimedia.org/api/rest_v1/media/math/render/svg/d522057eb0209091646bcc9221b24d9e13d0e371)
![{\displaystyle {\widetilde {\mathbf {z} }}_{t_{n}}^{i}:\left[-\tau _{i},0\right]\longrightarrow \mathbb {R} ^{d}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/d3f4bc0cc5b10a10d0844d2265aababd5028d6e4)
![{\displaystyle {\widetilde {\mathbf {z} }}_{t_{n}}^{i}(s):={\widetilde {\mathbf {z} }}^{i}(t_{n}+s),{\text{ }}s\in [-\tau _{i},0],}](https://wikimedia.org/api/rest_v1/media/math/render/svg/d4c1d308208faf6bc44261b2ce1b128a9f394195)
![{\displaystyle {\widetilde {\mathbf {z} }}^{i}:\left[t_{n}-\tau _{i},t_{n}\right]\longrightarrow \mathbb {R} ^{d}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/922ac10d7369b3ff0ffe9487ebb8a509a1305574)
![{\displaystyle t\in \lbrack t_{n}-\tau _{i},t_{n}]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/1f287009ce5d61c9d7ed2a843c93dfdef6226a50)
![{\displaystyle u\in \lbrack 0,h_{n}])}](https://wikimedia.org/api/rest_v1/media/math/render/svg/956d5d79b8e61251940a26652b1ecead7daa8b6c)
![{\displaystyle \mathbf {r} ^{\intercal }=\left[{\begin{array}{ll}\mathbf {0} _{1\times (d+1)}&1\end{array}}\right],h_{n}\leq \tau }](https://wikimedia.org/api/rest_v1/media/math/render/svg/05fdb2050667c7b12025db96c59b015211204a57)
![{\displaystyle t\in \lbrack t_{0},t_{n}]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/8ed17e1a454d79e9e85aa3fdb316727d17eb069b)
![{\displaystyle t\in \left[t_{0}-\tau ,t_{0}\right]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/99fafe4165d02d310ce7c3f55c6fbfc1742d2d65)
![{\displaystyle {\frac {d\mathbf {x} \left(t\right)}{dt}}=\mathbf {f} (\mathbf {x} (t),\mathbf {\xi } (t)),\quad t\in \left[t_{0},T\right],\qquad \qquad \qquad (6.1)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/06e7f4c5cbe530e829f8edf84e2e49d0829a2ef0)
![{\displaystyle t\in \left[t_{0},T\right].}](https://wikimedia.org/api/rest_v1/media/math/render/svg/aedde214e225a352a24c1c922043d56a7d90588b)
![{\displaystyle \mathbf {M} _{n}=\left[{\begin{array}{ccc}\mathbf {f} _{\mathbf {x} }\left(\mathbf {y} _{n},\mathbf {\xi } (t_{n})\right)&\mathbf {f} _{\mathbf {\xi } }(\mathbf {y} _{n},\mathbf {\xi } (t_{n})(\mathbf {\xi } (t_{n+1})-\mathbf {\xi } (t_{n}))/h_{n}&\mathbf {f} \left(\mathbf {y} _{n},\mathbf {\xi } (t_{n})\right)\\0&0&1\\0&0&0\end{array}}\right]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/de2e549c00febeadc80e1df68d44a6a66abd159a)
![{\displaystyle d\mathbf {x} (t)=\mathbf {f} (t,\mathbf {x} (t))dt+\sum \limits _{i=1}^{m}\mathbf {g} _{i}(t)d\mathbf {w} ^{i}(t),\quad t\in \left[t_{0},T\right],\qquad \qquad \qquad (7.1)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/464be59af5611f27dd6d946030a7202ce481c427)
![{\displaystyle t\in \lbrack t_{n},t_{n+1}]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/b15f801777044ac19f6f4e6551c3e78da54f33c5)
