Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- %% ========================================================================
- % dydt = ...; % define your continuation equation (vectorized, so don’t need to tell COCO otherwise)
- x0 = zeros(22,1);
- p0 = [0.0001];
- p_range = [p0 10]; % parameter range for continuation
- prob = coco_prob(); % define a new problem structure
- ode_fcns = {@Coupled_Model, @Coupled_Modeldx,@Coupled_Modeldp}; % can append jacobian to this cell if you have it explicitly (ie. {dydt, dfdt, dfdp})
- eq_path = {'path', 'to', 'isol2ep'}; % savdir/’run name’ (eg. /wd/data/path/to/isol2ep) by default sav dir is /wd/data/’run name’. input your path,
- % helpful for saving different runs under different folder structures.
- prob = ode_isol2ep(prob, '', ode_fcns{:}, x0, 'N(beta)', p0); % initialize ‘initial sol to equilibrium point’ toolbox. x0, p0 from ode45/ode15s/etc…
- % prob = coco_set(prob, 'cont', 'h_min', 0.1,'h0', 1,'h_max',10); % stepsize preferences, can play around with these
- %
- % prob = coco_set(prob, 'cont', 'ItMX', [0 1e3], 'NPR', 1); % max steps in run. print info every step.
- %
- % prob = coco_set(prob, 'corr', 'TOL', 1e-3,'ResTOL', 1e-3); % correction tolerance, can play around with these
- disp(' ')
- disp(' initialising eq continuation.')
- disp('-------------------------------------')
- bd_eq = coco(prob, eq_path, [], 1, 'N(beta)', p_range); % 1 is the solution dimension
- %% ========================================================================
- %///setup po continuation
- labs = coco_bd_labs(bd_eq, 'HB'); % get the labels of the hopf bifurcations
- bd_po = cell(1,length(labs));
- if ~isempty(labs)
- disp(' ')
- disp(' initialising po continuation.')
- disp('-------------------------------------')
- disp(' ')
- for i = 1:length(labs) % you can parallelise continuation of multiple bifurcations to speed up the process if needed
- prob = coco_prob(); % new continuation problem structure
- prob = ode_HB2po(prob, '', eq_path, labs(i)); % HB2po toolbox, follows PO from HB
- % prob = coco_add_slot(prob, 'get_PO', @coco_savpo, [], 'bddat'); % needed to save complete PO data (but not min/max), can help you with @coco_savpo if needed
- prob = coco_add_event(prob, 'HB', 'BP', 'atlas.test.BP', 0); % stops continuation repeating itself if it rejoins the eq branch
- prob = coco_set(prob, 'corr', 'TOL', 1e-3,'ResTOL', 1e-3, 'ItMX', 15);
- prob = coco_set(prob, 'coll', 'NTST', 80); % number of mesh points on which the collocation is performed
- prob = coco_set(prob, 'cont', 'NAdapt', 5, 'NPR', 1, 'PtMX', 1000); % mesh adaption every 5 steps
- prob = coco_set(prob,...
- 'cont', 'h_min', 0.1,'h0', 1,'h_max', 5,'ItMX', 1e3, 'al_max', 45);
- disp('done.')
- disp(' ')
- disp([' starting po continuation at ', datestr(clock)])
- disp('--------------------------------------------------------')
- bd_po{i} = coco(prob,...
- {'path', 'to', strcat(sprintf('hb2po_%d', i))}, [], 1, {'N(beta)' 'po.period'}, p_range); % output N(beta) and PO period during continuation
- end
- else
- disp('no hopf points found for po continuation')
- end
- %% ========================================================================
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement