Continuation and stability of steady states (equilibria)

(c) DDE-BIFTOOL v. 3.1.1(75), 31/12/2014

Once the user-defined functions are prepared, DDE-Biftool can compute and continue equilibria of the DDE, and compute their linearized stability, thus detecting local bifurcations. This demo requires demo1_funcs.html to have run beforehand.

Contents

%#ok<*ASGLU,*NOPTS,*NASGU>
%

Initial guess for steady state

It is clear that the neuron DDE has a steady state solution $(x_1^*,x_2^*)=(0,0)$ for all values of the parameters. We define a first steady state solution using the parameter values $\kappa=0.5$, $\beta=-1$, $a_{12}=1$, $a_{21}=2.34$, $\tau_1=\tau_2=0.2$ and $\tau_s=1.5$. Remember that we chose par $=[\kappa, \beta, a_{12}, a_{21},\tau_1,\tau_2, \tau_s]$.

stst.kind='stst';
stst.parameter=[1/2, -1, 1, 2.34, 0.2, 0.2, 1.5];
stst.x=[0;0]

stst = 
         kind: 'stst'
    parameter: [0.5 -1 1 2.34 0.2 0.2 1.5]
            x: [2x1 double]

Linear stability of initial equilibrium

We get default point method parameters and correct the point, which, being already a correct solution, remains unchanged. Computing and plotting stability of the corrected point reveals it has one unstable real mode, see figure.

flag_newhheur=1; % flag_newhheur=1 is the default choice if this argument is omitted
method=df_mthod(funcs,'stst',flag_newhheur);
method.stability.minimal_real_part=-1
[stst,success]=p_correc(funcs,stst,[],[],method.point)
% compute its stability:
stst.stability=p_stabil(funcs,stst,method.stability)
figure(1); clf;
p_splot(stst); % plot its stability:

method = 
    continuation: [1x1 struct]
     bifurcation: [1x1 struct]
           point: [1x1 struct]
       stability: [1x1 struct]
stst = 
         kind: 'stst'
    parameter: [0.5 -1 1 2.34 0.2 0.2 1.5]
            x: [2x1 double]
success =
     1
stst = 
         kind: 'stst'
    parameter: [0.5 -1 1 2.34 0.2 0.2 1.5]
            x: [2x1 double]
    stability: [1x1 struct]

Ask for roots with more negative real part

In both figures, approximations $(\times)$ and corrections $(*)$ are nearly indistinguishable.

method.stability.minimal_real_part=-2;
stst.stability=p_stabil(funcs,stst,method.stability); % recompute stability:
figure(2); clf;
p_splot(stst); % replot stability

Figures: Spectrum of equilibrium

Approximated $(\times)$ and corrected $(*)$ roots of the characteristic equation of neuron system at its steady state solution $(x_1^*,x_2^*)=(0,0)$. Real parts computed up to $\Re(\lambda)\geq -\frac{1}{\tau}$ (top), $\Re(\lambda)\geq -2$ (bottom).

Initialize branch of trivial equilibria

We will use this point as a first point to compute a branch of steady state solutions. First, we obtain an empty branch with free parameter $a_{21}$, limited by $a_{21}\in[0,5]$ and $\Delta a_{21}\leq 0.2$ between points.

% get an empty branch with ind_a21 as a free parameter:
branch1=df_brnch(funcs,ind_a21,'stst')
branch1.parameter
branch1.parameter.min_bound
% set bounds for continuation parameter
branch1.parameter.min_bound(1,:)=[ind_a21 0];
branch1.parameter.max_bound(1,:)=[ind_a21 5];
branch1.parameter.max_step(1,:)=[ind_a21 0.2];
% use stst as a first branch point:
branch1.point=stst;

branch1 = 
       method: [1x1 struct]
    parameter: [1x1 struct]
        point: []
ans = 
         free: 4
    min_bound: [3x2 double]
    max_bound: []
     max_step: []
ans =
     5     0
     6     0
     7     0

Extend and continue branch of trivial equilibria

To obtain a second starting point we change parameter value $a_{21}$ slightly and correct again.Because we know how the branch of steady state solutions continued in $a_{21}$ looks like (it is constant at $(x_1^*,x_2^*)=(0,0)$) we disable plotting during continuation by setting the corresponding continuation method parameter to zero.

stst.parameter(ind_a21)=stst.parameter(ind_a21)+0.1;
[stst,success]=p_correc(funcs,stst,[],[],method.point)
% use as a second branch point:
branch1.point(2)=stst;
branch1.method.continuation.plot=0;

% continue in one direction:
[branch1,s,f,r]=br_contn(funcs,branch1,100)
% turn the branch around:
branch1=br_rvers(branch1);
% continue in the other direction:
[branch1,s,f,r]=br_contn(funcs,branch1,100)

stst = 
         kind: 'stst'
    parameter: [0.5 -1 1 2.44 0.2 0.2 1.5]
            x: [2x1 double]
    stability: []
success =
     1
BR_CONTN warning: boundary hit.
branch1 = 
       method: [1x1 struct]
    parameter: [1x1 struct]
        point: [1x16 struct]
s =
    15
f =
     0
r =
     0
BR_CONTN warning: boundary hit.
branch1 = 
...

Stability of branch of equilibria

During continuation, sixteen points were successfully computed ($s=16$) before the right boundary $a_{21}=5$ was hit (signalled by a warning). No corrections failed ($f=0$) and no computed points were later rejected ($r=0$). Reversing the order of the branch points allows to continue to the left.

After obtaining suitable measure structures we plot the real part of the approximated and corrected roots of the characteristic equation along the branch, (see figure). Notice the strange behaviour (coinciding of several complex pairs of roots) at $a_{21}=0$. At this parameter value one of the couplings between the neurons is broken. In fact, for $a_{21}=0$, the evolution of the second component is independent of the evolution of the first.

branch1.method.stability.minimal_real_part=-2;
branch1=br_stabl(funcs,branch1,0,0);

% obtain suitable scalar measures to plot stability along branch:
[xm,ym]=df_measr(1,branch1)
figure(3); clf;
br_plot(branch1,xm,ym,'b'); % plot stability along branch:
ym.subfield='l0';
br_plot(branch1,xm,ym,'c');
plot([0 5],[0 0],'-.');
axis([0 5 -2 1.5]);
xlabel('a21');ylabel('\Re\lambda');
% plot stability versus point number:
figure(4); clf;
br_plot(branch1,[],ym,'b');
br_plot(branch1,[],ym,'b.');
plot([0 30],[0 0],'-.');
xlabel('point number along branch');ylabel('\Re(\lambda)');

xm = 
       field: 'parameter'
    subfield: ''
         row: 1
         col: 4
        func: ''
ym = 
       field: 'stability'
    subfield: 'l1'
         row: 'all'
         col: 1
        func: 'real'

Figures: Stability of equilibria

Real parts of the approximated (top) and corrected (top,bottom) roots of the characteristic equation versus $a_{21}$ (top) respectively the point number along the branch (bottom).

Save and continue with Hopf bifurcations: demo1_hopf.html

save('demo1_stst_results.mat');


Published with MATLAB® 7.14