IFAC2020/PhaseODE_NonCycling.m at master · giansimone/IFAC2020 · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
%% PhaseODE_NonCycling.m
%%% MARCH 4, 2020

function dxdt = PhaseODE_NonCycling(~, x, N, L, D, kappa, h)

%% Set model parameters
omega = 2 * pi / 80;

Theta_c = .25 * 2 * pi;


alpha_a = 20;

d_e = 1; % Degradation rate of external molecule

d_i = 1; % Degradation rate of internal molecule

d_a = 1;

alpha_i = .05;

eta = 2; % Coupling constant


%% Preallocate array dxdt, Theta, A, I, E
dxdt = nan(4*N,1);

Theta = zeros(N,1);

A = zeros(N,1);

I = zeros(N,1);

E = zeros(N,1);


%% Retrieve state
Theta(:,1) = mod(x(1:4:4*N), 2*pi);

A(:,1) = x(2:4:4*N);

I(:,1) = x(3:4:4*N);

E(:,1) = x(4:4:4*N);


%% Compute switching function
sigma = Theta < Theta_c;

u = zeros(N,1);

if h == +Inf

    u = double(I >= kappa);

else

    u = (I .^ h) ./ (I .^ h + kappa .^ h);

end


%% ODEs
dxdt((4*1-3):4:(4*N-3),1) = omega .* (1 + sigma .* (u - 1));

dxdt((4*1-2):4:(4*N-2),1) = - d_a .* A + alpha_a .* sigma;

dxdt((4*1-1):4:(4*N-1),1) = - d_i .* I + alpha_i .* A + eta .* (E - I);

dxdt(4*1:4:4*N,1) = - d_e .* E + eta .* (I - E) - D .* L * E;