## StochasticSemiDiscretizationMethod.jl

Julia package to investigate the behaviour of the first and second moments of stochastic linear delay differential equations
Author HTSykora
Popularity
1 Star
Updated Last
2 Years Ago
Started In
April 2019

# StochasticSemiDiscretizationMethod.jl

Julia package to investigate the behaviour of the first and second moments of stochastic linear delay differential equations based on the paper [1] Stochastic semi‐discretization for linear stochastic delay differential equations and the book [2] Semi-Discretization for Time-Delay Systems (by Insperger and Stepan).

This package provides a tool to approximate the stability properties and stationary behaviour of linear periodic delay systems of the forms:

$$\mathrm{d} \mathbf{x}(t) = \left(\mathbf{A}(t) \mathbf{x}(t) + \sum_{j=1}^g \mathbf{B}(t) \mathbf{x}(t-\tau_j(t))+\mathbf{c}(t)\right)\mathrm{d}t + \sum_{k=1}^w\left(\boldsymbol{\alpha}^k(t) + \sum_{j=1}^g \boldsymbol{\beta}^k_j(t) \mathbf{x}(t-\tau_j(t)) + \boldsymbol{\sigma}^k(t) \right)\mathrm{d}W^k(t)$$

by transforming the underlying differential equation into the stochastic mapping:

$$\mathbf{y}_{n+1} = \left(\mathbf{F}n+\sum{k=1}^w\mathbf{G}^k_n\right)\mathbf{y}_n + \left(\mathbf{f}n + \sum{k=1}^w\mathbf{g}^k_n\right),$$

where $n$ refers to the discrete time $t_n = n \Delta t$, $\mathbf{F}_n$ is the deterministic mapping matrix constructed using $\mathbf{A}$, $\mathbf{B}$ and $\tau_j$, $\mathbf{G}^k_n$ are the stochastic mapping matrices constructed using $\boldsymbol{\alpha}(t)$, $\boldsymbol{\beta}(t)$ and $\tau_j$, $\mathbf{f}_n$ and $\mathbf{g}^k_n$ are the deterministic and stochastic additive vectors, constructed using $\mathbf{c}$ and $\mathbf{\sigma}^k$, respectively. The vector $\mathbf{y}_n$ is the discretized state space vector:

$$\mathbf y_{n} = \left(\mathbf{x}(t_n)^\top, \mathbf{x}(t_{n-1})^\top,\ldots,\mathbf{x}(t_{n-r})\right)^\top!.$$

The first moment dynamics is described by the expected value of the stochastic mapping, leading to a deterministic mapping:

$$\mathbb{E}\left(\mathbf{y}_{n+1}\right) = \mathbf{F}_n\left(\mathbf{y}_n\right)+\mathbf{f}_n,$$ while the second moment dynamics is described by the expected value of the outer product of the mapping:

$$\mathbf M_n = \mathbb{E}\left(\mathbf y_n \mathbf y_n^\top\right), \quad \Rightarrow \quad \mathbf m_n = \mathrm{vec}\left(\mathbf M_n\right) := \left[ M_{n,11}, M_{n,22},\dots , M_{n,12}, M_{n,23},\dots, M_{n,1,\left(r+1\right)d} \right]^\top,$$

$$\mathbf m_{n+1} =\mathbf H_n \, \mathbf m_n + \mathbf h_{1,n}\mathbb{E}\left(\mathbf y_n\right)+\mathbf h_n,$$

where coefficient matrices $\mathbf{H}n$, $\mathbf{h}{1,n}$ and additive vector $\mathbf{h}_n$ are constructed using the mapping matrices $\mathbf{F}_n$, $\mathbf{G}^k_n$ and vectors $\mathbf{f}_n$, $\mathbf{g}^k_n$, respectively. Note, that since the coefficient matrices of the original system are constant, the statistical properties of the coefficient matrices and additive vectors, hence the matrices $\mathbf{F}n$, $\mathbf H_n$, $\mathbf{h}{1,n}$ and vectors $\mathbf{f}_n$, $\mathbf{h}n$ are also constant. The integer $r$ is chosen in a way, that $r \Delta t \geq \max{t \in \left[ 0, P \right], j = 1 \ldots g}\tau_j(t)$ (the discretized "history function" contains all possible delayed values) and $d$ is the dimension of the state space $\left(\mathbf{x}(t) \in \mathbb{R}^d\right)$..

With the use of the discrete mappings, the moment stability of the original system can be investigated (approximately), by the spectral radius $\rho$ of the coefficient matrices $\mathbf{F}_n$ over a period:

$$\rho\left(\prod_{i=0}^{p-1}{F}_{n+i}\right): \quad \begin{matrix} <1 & \Rightarrow & \text{the mapping is stable}\\ >1 & \Rightarrow & \text{the mapping is unstable} \end{matrix}$$

Furthermore, the steady-state first and second moments can be determined as the fix points of corresponding moment mappings.

# Citing

If you use this package as part of your research, teaching, or other activities, we would be grateful if you could cite the paper [1] and the book it is based on (BibTeX entries):

@article{Sykora2019,
author = {Sykora, Henrik T and Bachrathy, Daniel and Stepan, Gabor},
doi = {10.1002/nme.6076},
journal = {International Journal for Numerical Methods in Engineering},
keywords = { stability, stochastic problems, time delay,differential equations},
number = {ja},
title = {{Stochastic semi-discretization for linear stochastic delay differential equations}},
url = {https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.6076},
volume = {0}
}

@book{Insperger2011,
author = {Insperger, Tam{\'{a}}s and St{\'{e}}p{\'{a}}n, G{\'{a}}bor},
doi = {10.1007/978-1-4614-0335-7},
isbn = {978-1-4614-0334-0},
publisher = {Springer New York},
series = {Applied Mathematical Sciences},
title = {{Semi-Discretization for Time-Delay Systems}},
volume = {178},
year = {2011}
}


# Usage with examples

## Installation

julia> ] add StochasticSemiDiscretizationMethod

## Stochastic Hayes equations

$$\mathrm{d}{x}(t) = a ,x(t)\mathrm{d}t + \left(\beta ,x(t-1) + 1\right) \mathrm{d}W(t),$$ Here

$$\mathbf{A}(t) \equiv \begin{bmatrix} a \end{bmatrix}, \quad \mathbf{B}_1(t) \equiv \begin{bmatrix}0\end{bmatrix}, \quad \mathbf{c}(t) \equiv \begin{bmatrix} 0 \end{bmatrix},$$ $$\boldsymbol{\alpha}^1(t) \equiv \begin{bmatrix} 0 \end{bmatrix}, \quad \boldsymbol{\beta}^1_1(t) \equiv \begin{bmatrix}\beta\end{bmatrix}, \quad \tau^1_1(t) \equiv 1 \quad \mathrm{and} \quad \boldsymbol{\sigma}^1(t) \equiv \begin{bmatrix} 1 \end{bmatrix}.$$

(First example in paper [1])

using StochasticSemiDiscretizationMethod
function createHayesProblem(a,β)
AMx =  ProportionalMX(a*ones(1,1));
τ1=1.
BMx1 = DelayMX(τ1,zeros(1,1));
noiseID = 1
αMx1 = stCoeffMX(noiseID,ProportionalMX(zeros(1,1)))
βMx11 = stCoeffMX(noiseID,DelayMX(τ1,β*ones(1,1)))
LDDEProblem(AMx,[BMx1],[αMx1],[βMx11],cVec,[σ])
end
hayes_lddep=createHayesProblem(-6.,2.); # LDDE problem for Hayes equation
method=SemiDiscretization(0,0.1) # 0th order semi discretization with Δt=0.1
τmax=1. # the largest τ of the system
# Second Moment mapping
mapping=DiscreteMapping_M2(hayes_lddep,method,τmax,n_steps=10,calculate_additive=true); #The discrete mapping of the system
@show spectralRadiusOfMapping(mapping); # spectral radius ρ of the mapping matrix (ρ>1 unstable, ρ<1 stable)
statM2=VecToCovMx(fixPointOfMapping(mapping), length(mapping.M1_Vs[1])); # stationary second moment matrix of the hayes equation (equilibrium position)
@show statM2[1,1]

# statM2[1,1] = 0.12505625506304247

### Stability borders of the Hayes Equation

using MDBM
using Plots
gr();
using LaTeXStrings
method=SemiDiscretization(0,0.1);
τmax=1.

n_steps=10))); # No additive term calculated

axis=[Axis(-9.0:1.0,:a),
Axis(-5.0:5.0,:β)]

iteration=4;
stab_border_points=getinterpolatedsolution(solve!(MDBM_Problem(foo,axis),iteration));

scatter(stab_border_points...,xlim=(-9.,1.),ylim=(-5.,5.),
label="",title="2nd moment stability border of the Hayes equation",xlabel=L"A",ylabel=L"$\beta$",
guidefontsize=14,tickfont = font(10),markersize=2,markerstrokewidth=0)

## Stochastic Linear Delay Oscillator

$$\dot{x}(t) = v(t)$$

$$\dot{v}(t) + 2\zeta v(t) + A x(t) = B x(t-2\pi) + \left(\alpha x(t) + \beta x(t-2\pi) + \sigma\right)\Gamma(t)$$

(Second example in paper [1])

function createSLDOProblem(A,B,ζ,α,β,σ)
AMx =  ProportionalMX(@SMatrix [0. 1.;-A -2ζ]);
τ1=2π
BMx1 = DelayMX(τ1,@SMatrix [0. 0.; B 0.]);
noiseID = 1
αMx1 = stCoeffMX(noiseID,ProportionalMX(@SMatrix [0. 0.; α 0.]))
βMx11 = stCoeffMX(noiseID,DelayMX(τ1,@SMatrix [0. 0.; β 0.]))
LDDEProblem(AMx,[BMx1],[αMx1],[βMx11],cVec,[σVec])
end
SLDOP_lddep=createSLDOProblem(1.,0.1,0.1,0.1,0.1,0.5); # LDDE problem for Hayes equation
method=SemiDiscretization(5,(2π+100eps())/10) # 5th order semi discretization with Δt=2π/10
τmax=2π # the largest τ of the system
# Second Moment mapping
mapping=DiscreteMapping_M2(SLDOP_lddep,method,τmax,n_steps=10,calculate_additive=true); #The discrete mapping of the system

@show spectralRadiusOfMapping(mapping); # spectral radius ρ of the mapping matrix (ρ>1 unstable, ρ<1 stable)
statM2=VecToCovMx(fixPointOfMapping(mapping), length(mapping.M1_Vs[1])); # stationary second moment matrix of the hayes equation (equilibrium position)
@show statM2[1,1] |> sqrt;

# statM2[1, 1] |> sqrt = 0.9042071549857947

### Stability borders of the Delay Oscillator

method=SemiDiscretization(5,2π/30);
τmax=2π+100eps()
idxs=[1,2,3:2:(StochasticSemiDiscretizationMethod.rOfDelay(τmax,method)+1)*2...]

# ζ=0.05, α=0.3*A, β=0.3*B
n_steps=30),nev=8)); # No additive term calculated

axis=[Axis(-1.0:0.6:5.0,:A),
Axis(LinRange(-1.5,1.5,12),:B)]

iteration=4;
stab_border_points=getinterpolatedsolution(solve!(MDBM_Problem(foo,axis),iteration));

scatter(stab_border_points...,xlim=(-1.,5.),ylim=(-1.5,1.4),
label="",title="2nd moment stability border of the Delay Oscillator",xlabel=L"A",ylabel=L"$B$",
guidefontsize=14,tickfont = font(10),markersize=2,markerstrokewidth=0)

### Used By Packages

No packages found.