dummy-link

SemiDiscretizationMethod

Semi Discretization method for Periodic Linear Delay Differential Equations

Readme

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:

by transforming the underlying differential equation into the mapping:

where is the discrete time , is the mapping matrix constructed using , and , while the vector is the discretized state space vector:

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

Furthermore, the integer is chosen in a way, that (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 of the product of the mapping matrices over a period:

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,
address = {New York, NY},
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}},
url = {http://link.springer.com/10.1007/978-1-4614-0335-7},
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

Here


(Page 13 of the book)

using SemiDiscretizationMethod
function createHayesProblem(a,b)
    AMx =  ProportionalMX(a*ones(1,1));
    τ1=1. 
    BMx1 = DelayMX(τ1,b*ones(1,1));
    cVec = Additive(ones(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)

# spectralRadiusOfMapping(mapping) = 0.941189374166563
# 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.

foo(a,b) = log(spectralRadiusOfMapping(DiscreteMapping(createHayesProblem(a,b),method,τmax,
    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>
<img src="https://latex.codecogs.com/gif.latex?\tau_1(t)&space;\equiv&space;2\pi,&space;\quad&space;\text{and}&space;\quad&space;\mathbf{c}(t)&space;=&space;\begin{bmatrix}&space;0&space;\&space;\sin(2t)&space;\end{bmatrix}." title="\tau_1(t) \equiv 2\pi, \quad \text{and} \quad \mathbf{c}(t) = \begin{bmatrix} 0 \ \sin(2t) \end{bmatrix}." />

(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)
# spectralRadiusOfMapping(mapping) = 0.5131596340374617
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);

foo(δ,b0) = log(spectralRadiusOfMapping(DiscreteMapping(createMathieuProblem(δ,ε,b0,a1,T=T),method,τmax,
    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)

First Commit

04/02/2019

Last Touched

11 months ago

Commits

13 commits

Used By: