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 form:
by transforming the underlying differential equation into the stochastic mapping:
where
is the discrete time
,
is the deterministic mapping matrix constructed using
,
and
,
are the stochastic mapping matrices constructed using
,
and
,
and
are the deterministic and stochastic additive vectors, constructed using
and
, respectively.
The vector
is the discretized state space vector:
The first moment dynamics is described by the expected value of the stochastic mapping, leading to a deterministic mapping:
while the second moment dynamics is described by the expected value of the outer product of the mapping:
where coefficient matrice
,
and additive vector
are constructed using the mapping matrices
,
and vectors
,
, 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
,
,
and vectors
,
are also constant.
The integer
is chosen in a way, that
(the discretized "history function" contains all possible delayed values) and
is the dimension of the state space
.
With the use of the discrete mappings, the moment stability of the original system can be investigated (approximately), by the spectral radius
of the coefficient matrices
over a period:
Furthermore, the stationary 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,
address = {New York, NY},
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}},
url = {http://link.springer.com/10.1007/978-1-4614-0335-7},
volume = {178},
year = {2011}
}
Usage with examples
Installation
julia> ] add StochasticSemiDiscretizationMethod
Stochastic Hayes equations
Here
(First example in paper [1])
using StochasticSemiDiscretizationMethod
function createHayesProblem(a,β)
AMx = ProportionalMX(a*ones(1,1));
τ1=1.
BMx1 = DelayMX(τ1,zeros(1,1));
cVec = Additive(1)
noiseID = 1
αMx1 = stCoeffMX(noiseID,ProportionalMX(zeros(1,1)))
βMx11 = stCoeffMX(noiseID,DelayMX(τ1,β*ones(1,1)))
σ = stAdditive(1,Additive(ones(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]
# spectralRadiusOfMapping(mapping) = 0.3835887415448961
# 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.
foo(a,b) = log(spectralRadiusOfMapping(DiscreteMapping_M2(createHayesProblem(a,b),method,τmax,
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
Here
(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.]);
cVec = Additive(2)
noiseID = 1
αMx1 = stCoeffMX(noiseID,ProportionalMX(@SMatrix [0. 0.; α 0.]))
βMx11 = stCoeffMX(noiseID,DelayMX(τ1,@SMatrix [0. 0.; β 0.]))
σVec = stAdditive(1,Additive(@SVector [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;
# spectralRadiusOfMapping(mapping) = 0.5059591707964441
# 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
foo(A,B) = log(spectralRadiusOfMapping(DiscreteMapping_M2(createSLDOProblem(A,B,0.05,0.3*A,0.3*B,0.),method,τmax,idxs,
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)