-
Notifications
You must be signed in to change notification settings - Fork 0
/
main_pE_all.m
132 lines (111 loc) · 3.72 KB
/
main_pE_all.m
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
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
main_init;
addpath('C:\Users\bai_f\Documents\research\Axon_growth\DelayedFeedback\mine_dde23');
%% Initial values of bifurcation parameters
% tau = N/nu/(1-rho_bulk)
% N==1, nu==1, rho_bulk==1/2, then we have the following initial values.
rhobulkK = 0.5;
rhobulkD = 0.5;
pE = 1;
wE = 5;
pI = 6;
wI = 5;
vK = 1;
vD = 1;
N = 10;
contpar=ind_pE;
namepar='pE';
file_id=4;
%% Setup the first fixed point
stst.kind = 'stst';
stst.parameter = [rhobulkK, rhobulkD, pE, wE, pI, wI, vK, vD, N];
% This values are approximations obtained from dde23.
E1_st=0.444444435008581;
E2_st=0.555555548432940;
I1_st=0.0197344203680364;
I2_st=0.0157996522811142;
stst.x = [E1_st; E2_st; I1_st; I2_st];
% Correct the approximations above.
method = df_mthod(funcs, 'stst');
% If there is an error message caused by running the following line,
% you should check set_funcs to see whether the parameter positions are
% correct or not.
[stst, success] = p_correc(funcs, stst, [], [], method.point);
E1_st = stst.x(1);
E2_st = stst.x(2);
I1_st = stst.x(3);
I2_st = stst.x(4);
%% Initialize the branch
stst_branch0 = SetupStst(funcs,'x',[E1_st; E2_st; I1_st; I2_st],'parameter',stst.parameter,...
'contpar',contpar,'max_step',[contpar,0.1],'min_bound',...
[contpar 1],'max_bound',[contpar 50],...
'newheuristics_tests',0);
%% Continuation of the branch
figure('Name',append('Generating branch for ',namepar),'NumberTitle','off');
clf;
ax1=gca;
title(ax1,sprintf(append('Generating branch for ',namepar)));
[stst_branch0] = br_contn(funcs,stst_branch0,500);
[stst_branch_wbifs,stst_testfuncs]=LocateSpecialPoints(funcs,stst_branch0);
nunst_stst=GetStability(stst_branch_wbifs);
%% Plot bifurcation diagram
par_stst=getpar(stst_branch_wbifs,contpar);
E1_stst=getx(stst_branch_wbifs,1);
fig_bf = figure('Name',append('Bifurcation diagram for ',namepar),'NumberTitle','off');
clf;
ax1=gca;
hold on;
plot(ax1,par_stst(nunst_stst==0),E1_stst(nunst_stst==0),'g.','DisplayName','stable');
plot(ax1,par_stst(nunst_stst>0),E1_stst(nunst_stst>0),'b.','DisplayName','unstable');
plot(ax1,bgetpar(stst_branch_wbifs,contpar,'hopf'),bgetx(stst_branch_wbifs,1,'hopf'),'ks','DisplayName','hopf');
legend(ax1,'location','north');
xlabel(ax1,namepar);
ylabel(ax1,'E1');
title(ax1,sprintf(append('Bifurcation diagram for ',namepar)));
%% Extract frequency and amplitude of oscillation
fprintf('----- Extract frequencies -----\n');
JK = vK*J(rhobulkK);
JD = vD*J(rhobulkD);
tauK = delay(rhobulkK, N, vK);
tauD = delay(rhobulkD, N, vD);
delays = [tauK, tauD];
tspan = [0, 400];
t_start = 100;
t_end = 400;
dt = 0.1;
length_par = length(par_stst);
freq_list = [];
max_list = [];
min_list = [];
amp_list = [];
par_list = [];
for i=1:length_par
if nunst_stst(i)>0
pE = par_stst(i)
sol = dde23(@rhs_dde23, delays, @history, tspan);
[freq, amp, T, amp_max, amp_min] = extract_freq_amp(sol, 1, t_start, t_end, dt);
if amp<1e-2
freq = 0;
end
par_list(end+1) = pE;
freq_list(end+1) = freq;
max_list(end+1) = amp_max;
min_list(end+1) = amp_min;
amp_list(end+1) = amp;
end
end
figure('Name','freq','NumberTitle','off')
plot(par_list, freq_list, 'b.');
xlabel('pE');
ylabel('Frequency');
figure('Name','amp','NumberTitle','off')
plot(par_list, amp_list, 'b.');
xlabel('pE');
ylabel('Amplitude');
figure(fig_bf);
ax1 = gca;
hold on;
plot(ax1,par_list,max_list,'k.', 'DisplayName','max');
plot(ax1,par_list,min_list,'r.', 'DisplayName','min');
legend(ax1,'location','north');
%% output
save_data_1par(namepar, file_id, N, par_stst, E1_stst, nunst_stst, par_list, freq_list, amp_list, max_list, min_list);