Continuation of folds of periodic orbits
(c) DDE-BIFTOOL v. 3.1.1(75), 31/12/2014
The extension ddebiftool_extra_psol is able to continue local bifurcations of periodic orbits in two parameters. This demo shows how one can continue folds (saddle-nodes) of periodic orbits i nthe neuron problem. (requires running of demo1_psol.html first).
The extension is installed in a separate folder. In addition, a folder with utilities for user convenience is loaded.
addpath('../../ddebiftool_extra_psol/'); addpath('../../ddebiftool_utilities/'); %#ok<*ASGLU,*NOPTS,*NASGU>
The functions neuron_sys_rhs and neuron_sys_deri are not vectorized. In order to speed up computations we re-define neuron_sys_rhs, and replace neuron_sys_seri with the default finite-difference approximation
neuron_sys_rhs=@(xx,par)[... -par(1)*xx(1,1,:)+par(2)*tanh(xx(1,4,:))+par(3)*tanh(xx(2,3,:));.... -par(1)*xx(2,1,:)+par(2)*tanh(xx(2,4,:))+par(4)*tanh(xx(1,2,:))]; vfuncs=set_funcs(... 'sys_rhs',neuron_sys_rhs,... 'sys_tau',@()[5,6,7],... 'x_vectorized',true);
For the branch computed in demo1_psol.html one fold occured at the maximal parameter value. So we extract its index indmax.
Then we call SetupPOfold to initialize the branch and to create the functions for the extended system. For the core DDE-Biftool routines the fold of periodic orbits is of the same type as a standard periodic orbit. However, the user-provided right-hand side has been extended (eg, foldfuncs.sys_rhs is different from funcs.sys_rhs). SetupPOfold has three mandatory parameters:
- funcs containing the user-problem functions (in the demo: vfuncs),
- branch, the branch on which the fold occurred (in the demo: branch5), and
- ind, the index in branch.point near which the fold occurs (in the demo indmax).
Output (to be fed into br_contn):
- foldfuncs: functions for the extended system (derived from funcs)
- FPObranch: branch of folds of periodic orbits with two points
The other parameters are optional name-value pairs. Important are
- 'contpar': typically two integers ,the two continuation parameters,
- 'dir': in which direction the first step is taken), and
- 'step': the length of the first step. Note that name-value pairs are also passed on to fields of the branch structure: the output FPObranch inherits all fields from the input branch. Additional name-value argument pairs canbe used to change selected fields.
[foldfuncs,branch7]=SetupPOfold(vfuncs,branch5,indmax,'contpar',[ind_a21,ind_taus],... 'dir',ind_taus,'print_residual_info',1,'step',0.01,'plot_measure',,... 'min_bound',[ind_a21,0],'max_bound',[ind_a21,4; ind_taus,10],... 'max_step',[ind_a21,0.2; ind_taus,0.5]);
it=1, res=0.000574744 it=2, res=0.0083463 it=3, res=0.000139207 it=4, res=4.6582e-08 it=5, res=2.22616e-10 it=1, res=0.000267241 it=2, res=2.96116e-10 it=1, res=0.0699482 it=2, res=0.014741 it=3, res=6.29929e-06 it=4, res=2.28572e-10 it=1, res=0.000175272 it=2, res=8.75278e-10
The output of SetupPOfold can be fed directly into br_contn to perform a branch continuation. Note that the computation is typically slower than a standard continuation of periodic orbits because the system dimension is twice as large and additional delays have been introduced.
figure(6); branch7.method.point.print_residual_info=0; branch7=br_contn(foldfuncs,branch7,100); xlabel('a21');ylabel('tau_s'); title('Continuation of fold of periodic orbits');
BR_CONTN warning: boundary hit.
Since the fold continuation solves an extended system (with additional components) it does not make sense to compute the linear stability of the resulting orbits directly. However, foldfuncs has an additional field 'get_comp'. This is a function that extacts the components of the extended system that correspond to the solution of the original system. The output of get_comp is an array of periodic orbits that are located precisely on the fold.
We compute Floquet multipliers for these orbits using the utility function GetStability(psol_array,...). Its optional inputs instruct it to ignore the two eigenvalues closest to 1 when computing stability. It also returns triv_defect which measures the distance of the supposedly trivial eigenvalues (as computed) to their correct values. GetStability needs the optional input 'funcs' if its first argument does not yet have a stability structure.
[nunst_pf,dom,triv_defect,pf_orbits]=GetStability(pf_orbits,'funcs',vfuncs,... 'exclude_trivial',true,'locate_trivial',@(x)[1,1]); fprintf('max number of unstable Floquet multipliers: %d\n',max(nunst_pf)); n_orbits=length(pf_orbits); figure(16); a21_pfold=arrayfun(@(x)x.parameter(4),branch7.point); taus_pfold=arrayfun(@(x)x.parameter(7),branch7.point); plot(1:n_orbits,triv_defect); xlabel('point number'); ylabel('defect of trivial Floquet multiplier');
max number of unstable Floquet multipliers: 0
Save results (end of tutorial demo, but try also demo1_hcli.html)