0

0

0

1

# SemiDiscretizationMethod.jl

Julia package to investigate the behaviour of linear delay differential equations based on the book 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:

$\dot{\mathbf{x}}(t)&space;=&space;\mathbf{A}(t)&space;\mathbf{x}(t)&space;+&space;\sum_{j=1}^g&space;\mathbf{B}_j(t)&space;\mathbf{x}(t-\tau_j(t))+\mathbf{c}(t)$

by transforming the underlying differential equation into the mapping:

$\mathbf{y}_{n+1}&space;=&space;\mathbf{F}_n\mathbf{y}_n+\mathbf{f}_n,$

where $n$ is the discrete time $(t_n&space;=&space;n&space;\Delta&space;t)$ , $\mathbf{F}_n$ is the mapping matrix constructed using $\mathbf{A}(t)$ , $\mathbf{B}(t)$ and $\tau_j(t)$ , while the vector $\mathbf{y}_n$ is the discretized state space vector:

$\mathbf{y}_n&space;=&space;\left(\mathbf{x}(t_n)^\top,&space;\mathbf{x}(t_{n-1})^\top,\ldots,\mathbf{x}(t_{n-r})\right)^\top\!.$

Each coefficient matrices of delay differential equations are periodic, with a principle period of $P$ , namely:

$A(t)=A(t+P),\;&space;B_j(t)=B_j(t+P),\;&space;\tau_j(t)=\tau_j(t+P))&space;and&space;c(t)=c(t+P).$ Furthermore, the integer $r$ is chosen in a way, that $\inline&space;r\Delta&space;t\geq&space;\max_{t&space;\in&space;\left[0,P\right],j=1\ldots&space;g}\tau_j(t)$ (the discretized "history function" contains all possible delayed values).

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

$\rho\left(\mathbf{F}_{n+p-1}\cdot\mathbf{F}_{n+p-2}\cdot\ldots\cdot\mathbf{F}_{n}\right)\left\{&space;\begin{matrix}&space;<1&space;&&space;\Rightarrow&space;&&space;\text{the&space;mapping&space;is&space;stable}\\&space;>1&space;&&space;\Rightarrow&space;&&space;\text{the&space;mapping&space;is&space;unstable}&space;\end{matrix}&space;\right.$

Furthermore, the stationary solution can be determined by the periodic fix point (stationary orbit) of the mapping.

# Citing

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

@book{Insperger2011,
author = {Insperger, Tam{julia-observer-html-cut-paste-0__work#39;{a}}s and St{julia-observer-html-cut-paste-0__work#39;{e}}p{julia-observer-html-cut-paste-0__work#39;{a}}n, G{julia-observer-html-cut-paste-0__work#39;{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}
}


If you are interested in the behaviour of your linear delay model in the presence of Gaussian white noise, please consider the StochasticSemiDiscretizationMethod.jl package.

# Usage with examples

## Installation

julia> ] add SemiDiscretizationMethod


## Hayes equations

$\dot{x}(t)&space;=&space;a&space;\,x(t)&space;+&space;b&space;\,x(t-1)&space;+&space;1.$

Here

using SemiDiscretizationMethod

function createHayesProblem(a,b)
AMx =  ProportionalMX(a*ones(1,1));
τ1=1.
BMx1 = DelayMX(τ1,b*ones(1,1));
LDDEProblem(AMx,[BMx1],cVec)
end

hayes_lddep=createHayesProblem(-1.,-1.); # LDDE problem for Hayes equation
method=SemiDiscretization(1,0.1) # 1st order semi discretization with Δt=0.1
τmax=1. # the largest τ of the system
mapping=DiscreteMapping(hayes_lddep,method,τmax,n_steps=1,calculate_additive=true); #The discrete mapping of the system

@show spectralRadiusOfMapping(mapping); # spectral radius ρ of the mapping matrix (ρ>1 unstable, ρ<1 stable)
@show fixPointOfMapping(mapping); # stationary solution of the hayes equation (equilibrium position)

# fixPointOfMapping(mapping) = [0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5]


### Stability borders of the Hayes Equation

using MDBM

using Plots
gr();
using LaTeXStrings

method=SemiDiscretization(4,0.1);
τmax=1.

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

axis=[Axis(-15.0:15.,:a),
Axis(-15.0:15.,:b)]

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

scatter(stab_border_points...,xlim=(-15.,15.),ylim=(-15.,15.),
label="",title="Stability border of the Hayes equation",xlabel=L"a",ylabel=L"b",
guidefontsize=14,tickfont = font(10),markersize=2,markerstrokewidth=0)
 html
![](./assets/HayesStability.png)
## Delay Mathieu equation
<!-- $$\ddot{x}(t) + a_1 \,\dot{x}(t)+(\delta + \varepsilon \cos(t))x(t)= b_0 \,x(t-2\pi) + \sin(2t)$$ -->
<a href="https://www.codecogs.com/eqnedit.php?latex=\ddot{x}(t)&space;&plus;&space;a_1&space;\,\dot{x}(t)&plus;(\delta&space;&plus;&space;\varepsilon&space;\cos(t))x(t)=&space;b_0&space;\,x(t-2\pi)&space;&plus;&space;\sin(2t)" target="_blank"><img src="https://latex.codecogs.com/gif.latex?\ddot{x}(t)&space;&plus;&space;a_1&space;\,\dot{x}(t)&plus;(\delta&space;&plus;&space;\varepsilon&space;\cos(t))x(t)=&space;b_0&space;\,x(t-2\pi)&space;&plus;&space;\sin(2t)." title="\ddot{x}(t) + a_1 \,\dot{x}(t)+(\delta + \varepsilon \cos(t))x(t)= b_0 \,x(t-2\pi) + \sin(2t)." /></a>

Here
<!-- $$\mathbf{x}(t) = \begin{bmatrix}x(t) \ \dot{x}(t)\end{bmatrix}, \quad \mathbf{A}(t) = \begin{bmatrix} 0 & 1 \ -\delta - \varepsilon \cos(t) & -a_1 \end{bmatrix}, \quad \mathbf{B}_1(t) = \begin{bmatrix}0 & 0 \ b_0 & 0\end{bmatrix}, \quad \tau_1(t) \equiv 2\pi, \quad \text{and} \quad \mathbf{c}(t) = \begin{bmatrix} 0 \ \sin(2t) \end{bmatrix}.$$ -->
<img src="https://latex.codecogs.com/gif.latex?\mathbf{x}(t)&space;=&space;\begin{bmatrix}x(t)&space;\&space;\dot{x}(t)\end{bmatrix},&space;\quad&space;\mathbf{A}(t)&space;=&space;\begin{bmatrix}&space;0&space;&&space;1&space;\&space;-\delta&space;-&space;\varepsilon&space;\cos(t)&space;&&space;-a_1&space;\end{bmatrix},&space;\quad&space;\mathbf{B}_1(t)&space;=&space;\begin{bmatrix}0&space;&&space;0&space;\&space;b_0&space;&&space;0\end{bmatrix}," title="\mathbf{x}(t) = \begin{bmatrix}x(t) \ \dot{x}(t)\end{bmatrix}, \quad \mathbf{A}(t) = \begin{bmatrix} 0 & 1 \ -\delta - \varepsilon \cos(t) & -a_1 \end{bmatrix}, \quad \mathbf{B}_1(t) = \begin{bmatrix}0 & 0 \ b_0 & 0\end{bmatrix}," />
<br>

(Page 77 of the book)


function createMathieuProblem(δ,ε,b0,a1;T=2π) AMx = ProportionalMX(t->@SMatrix [0. 1.; -δ-ε*cos(2π/T*t) -a1]); τ1=2π # if function is needed, the use τ1 = t->foo(t) BMx1 = DelayMX(τ1,t->@SMatrix [0. 0.; b0 0.]); cVec = Additive(t->@SVector [0.,sin(4π/T*t)]) LDDEProblem(AMx,[BMx1],cVec) end;


julia
τmax=2π # the largest τ of the system
P=2π #Principle period of the system (sin(t)=sin(t+P))
mathieu_lddep=createMathieuProblem(3.,2.,-0.15,0.1,T=P); # LDDE problem for Hayes equation
method=SemiDiscretization(1,0.01) # 1st order semi discretization with Δt=0.01
# if P = τmax, then n_steps is automatically calculated
mapping=DiscreteMapping(mathieu_lddep,method,τmax,
n_steps=Int((P+100eps(P))÷method.Δt),calculate_additive=true); #The discrete mapping of the system

@show spectralRadiusOfMapping(mapping); # spectral radius ρ of the mapping matrix (ρ>1 unstable, ρ<1 stable)
fp=fixPointOfMapping(mapping); # stationary solution of the hayes equation (equilibrium position)

# fixPointOfMapping plotted
plot(0.0:method.Δt:P,fp[1:2:end],
xlabel=L"-s",title=L"t \in [nP,(n+1)P],\quad n \to \infty",guidefontsize=14,linewidth=3,
label=L"x(t-s)",legendfontsize=11,tickfont = font(10))
plot!(0.0:method.Δt:P,fp[2:2:end],
xlabel=L"-s",linewidth=3,
label=L"\dot{x}(t-s)")
plot!(0.0:method.Δt:P,sin.(2*(0.0:method.Δt:P)),linewidth=3,label=L"\sin(2t)")


### Stability Chart of the delayed Mathieu equation

a1=0.1;
ε=1;
τmax=2π;
T=1π;
method=SemiDiscretization(2,T/40);

n_steps=Int((T+100eps(T))÷method.Δt)))); # No additive term calculated

axis=[Axis(-1:0.2:5.,:δ),
Axis(-2:0.2:1.5,:b0)]

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

scatter(stab_border_points...,xlim=(-1.,5),ylim=(-2.,1.5),
label="",title="Stability border of the delay Mathieu equation",xlabel=L"\delta",ylabel=L"b_0",
guidefontsize=14,tickfont = font(10),markersize=2,markerstrokewidth=0)
 html
![](./assets/MathieuStability.png)
`

04/02/2019