-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathhopfstochcomb.m
141 lines (119 loc) · 4.63 KB
/
hopfstochcomb.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
133
134
135
136
137
138
139
140
141
function [Xdet, Xsto, Fext, pxxd, pxxs, fxx,pxxdcpk,pxxscpk,pxxsin1,cpk] = hopfstochcomb(mu,fosc,NoiseSTD,tvec,F1,F2,f1,f2)
%
% This function simulates the normal form of the supercritical Hopf
% bifurcation, given by two planar equations:
%
% x_dot = mu*x - omega*y - x*(x^2 + y^2) + F
% y_dot = omega*x + mu*y - y*(x^2 + y^2) + F
%
% where mu is the control parameter. For mu>0, the system will oscillate at
% an amplitude that grows with sqrt(mu). Alternatively, one may express the
% above equations in polar coordinates, making the amplitude relationship
% with respect to mu more apparent:
%
% rho_dot = rho*(mu + i*omega - rho^2)
%
% Here we simulate both the deterministic and stochastic cases for the
% supercritical Hopf bifurcation
%
% [Xdet Xsto] = hopfstoch(mu,fosc,NoiseSTD,tvec,Fextmax,fr)
%
% Xdet : deterministic result
% Xsto : stochastic result
%
%
% tvec : time vector
% mu : control parameter
% xNoiseSTD : standard deviation of stochastic noise in x
% yNoiseSTD : standard deviation of stochastic noise in y
% fosc : frequency of oscillation on the unstable side of the bifurcation
% Fextmax : amplitude of periodic forcing
% fr : frequency of driving
%
% By modifying the code, you can also add a step function or external
% forcing.
%
%
% Initial condition
xzero = 1;
yzero = -1;
% Add external forcing if desired
sinusoidalstim = 1; pulsestim = 0; % pulse or sinusoid?
%Fextmax = 1e8; % amplitude of sinusoidal stim OR pulse
%fr = 1.01; % frequency of stimulation
pulsestart = 1; % start of pulse
pulseend = 2; % end of pulse
% Decrease time step size by factor of Dtfac to ensure convergence
Dtfac = 10^2;
Dt = (tvec(2)-tvec(1))/Dtfac;
N = tvec(end)/Dt;
%Set the default random number stream
RandStream.setGlobalStream(RandStream('mt19937ar','seed',1))
xdW = sqrt(Dt)*randn(1,N); % White noise increments
ydW = sqrt(Dt)*randn(1,N); % White noise increments
xNoiseSTD = NoiseSTD;
yNoiseSTD = NoiseSTD;
xdet = zeros(1,N); xdet(1) = xzero;
xsto = zeros(1,N); xsto(1) = xzero;
ydet = zeros(1,N); ydet(1) = yzero;
ysto = zeros(1,N); ysto(1) = yzero;
% External forcing
if sinusoidalstim == 1
Ftime = linspace(tvec(1),tvec(end),N);
Fext = F1*sin(2*pi*f1*Ftime) + F2*sin(2*pi*f2*Ftime);
Fext = hilbert(Fext);
elseif pulsestim == 1
Ftime = linspace(tvec(1),tvec(end),N);
Fext = ((Ftime<pulseend)-(Ftime<pulsestart))*Fextmax;
else
Fext = zeros(1,N);
end
% Euler-Murayama Method with Ito Integration
for j = 2:N
%Deterministic integral
xdet(j) = xdet(j-1) + Dt*(mu*xdet(j-1) - 2*pi*fosc*ydet(j-1) - xdet(j-1)*(xdet(j-1)^2 + ydet(j-1)^2) + real(Fext(j)));
ydet(j) = ydet(j-1) + Dt*(2*pi*fosc*xdet(j-1) + mu*ydet(j-1) - ydet(j-1)*(xdet(j-1)^2 + ydet(j-1)^2) + imag(Fext(j)));
%Stochastic integral
xsto(j) = xsto(j-1) + Dt*(mu*xsto(j-1) - 2*pi*fosc*ysto(j-1) - xsto(j-1)*(xsto(j-1)^2 + ysto(j-1)^2) + real(Fext(j))) + xNoiseSTD*xdW(j);
ysto(j) = ysto(j-1) + Dt*(2*pi*fosc*xsto(j-1) + mu*ysto(j-1) - ysto(j-1)*(xsto(j-1)^2 + ysto(j-1)^2) + imag(Fext(j))) + yNoiseSTD*ydW(j);
end
Xdet = zeros(2,length(tvec)-1);
Xsto = zeros(2,length(tvec)-1);
%Return vectors at times specified by Time.
Xdet(1,:) = xdet(1:Dtfac:N);
Xdet(2,:) = ydet(1:Dtfac:N);
Xsto(1,:) = xsto(1:Dtfac:N);
Xsto(2,:) = ysto(1:Dtfac:N);
Fext = Fext(1:Dtfac:N);
L=length(Xdet);
NFFT = (2^3)*2^nextpow2(L);
Fs = 1/(tvec(2)-tvec(1));
[pxxd,fxx]=pwelch(Xdet(1,:),[],[],NFFT,Fs);
[pxxs,fxx]=pwelch(Xsto(1,:),[],[],NFFT,Fs);
[pxxsine,fxx2]=pwelch(1*sin(6*pi*tvec),[],[],NFFT,Fs);
qw=findnearest(fxx2,3);
pxxsin1 = pxxsine(qw-50:qw+50);
winfunc = rectwin(length(tvec));
pxxsin1 = sqrt(max(pxxsin1).*(2.*Fs.*length(tvec).*(sum(abs(winfunc).^2)./length(tvec))))./length(tvec);
F1=min([f1 f2]);
F2=max([f1 f2]);
cpk=find(fxx==(2*f1-f2));
if isempty(cpk)==1
cpk=findnearest(fxx,(2*f1-f2));
end
pxxdcpk=max(pxxd(cpk-50:cpk+50));
pxxscpk=max(pxxs(cpk-50:cpk+50));
pxxdcpk=(sqrt(pxxdcpk.*(2.*Fs.*length(tvec).*(sum(abs(winfunc).^2)./length(tvec))))./length(tvec))./pxxsin1;
pxxscpk=(sqrt(pxxscpk.*(2.*Fs.*length(tvec).*(sum(abs(winfunc).^2)./length(tvec))))./length(tvec))./pxxsin1;
% Make a plot of the data?
plotyn=0;
if plotyn==1
figure;
subplot(1,2,1);hold on;plot(tvec(2:end),Xsto(1,:),'r');plot(tvec(2:end),Xdet(1,:),'k');title('Black=deterministic; Red=stochastic; real part only');
subplot(1,2,2);hold on;plot(Xsto(1,:),Xsto(2,:),'r');plot(Xdet(1,:),Xdet(2,:),'k');title('Black=deterministic; Red=stochastic');
figure;
subplot(1,2,1);plot(fxx,pxxd,'k');axis([0 5 0 max(max(pxxd))]);
subplot(1,2,2);plot(fxx,pxxs,'r');axis([0 5 0 max(max(pxxs))]);
end
end