%% DDE-Biftool demo Mackey-Glass Equation % % The Mackey-Glass equation is given by % % $$x'(t)=\beta \frac{x(t-\tau)}{1+x(t-\tau)^n}-\gamma x(t)$$ % % Parameters are (in this order) |beta|, |n|, |tau| (|gamma| is not part of % parameter vector). % % % (c) DDE-BIFTOOL v. 3.1(73), 31/12/2014 % Extended by Tony Humphries 14/1/2016 % % %% load DDE-Biftool into path clear all; close all diaryname=strcat('diary',date,'.txt'); diary(diaryname) % change this and also the rmpath add end of file if % you have ddebiftool in a different folder. You might % also consider adding it permanently to your matlab path addpath('C:/Program Files/MATLAB/dde_biftool_v3.1/ddebiftool/',... 'C:/Program Files/MATLAB/dde_biftool_v3.1/ddebiftool_extra_psol/',... 'C:/Program Files/MATLAB/dde_biftool_v3.1/ddebiftool_extra_nmfm/',... 'C:/Program Files/MATLAB/dde_biftool_v3.1/ddebiftool_utilities/'); %% Enable vectorization % (disable for speed comparison) x_vectorize=true; %% Set user-defined functions % using |gamma| as constant and (|beta|,|n|,|tau|) as parameters gamma=1.0; beta_ind=1; n_ind=2; tau_ind=3; if x_vectorize f=@(x,xtau,beta,n)beta*xtau./(1+xtau.^n)-gamma*x; funcs=set_funcs(... 'sys_rhs',@(xx,p)f(xx(1,1,:),xx(1,2,:),p(1),p(2)),... 'sys_tau',@()tau_ind,... 'x_vectorized',true); else f=@(x,xtau,beta,n)beta*xtau/(1+xtau^n)-gamma*x; %#ok funcs=set_funcs(... 'sys_rhs',@(xx,p)f(xx(1,1,:),xx(1,2,:),p(1),p(2)),... 'sys_tau',@()tau_ind); end %% Initial parameters and state beta0=2; n0=4; tau0=2; x0=(beta0-1)^(1/n0); %% Initialization of branch of non-trivial equilibria contpar=n_ind; nontriv_eqs=SetupStst(funcs,'x',x0,'parameter',[beta0,n0,tau0],'step',0.1,... 'contpar',contpar,'max_step',[contpar,0.1],'max_bound',[contpar,21]); %% Compute and find stability of non-trivial equilibria disp('Trivial equilibria'); figure(1);clf nontriv_eqs=br_contn(funcs,nontriv_eqs,1000); nontriv_eqs=br_stabl(funcs,nontriv_eqs,0,1); nunst_eqs=GetStability(nontriv_eqs); ind_hopf=find(nunst_eqs<2,1,'last'); % plot char values at Hopf bif figure(100) evals=nontriv_eqs.point(ind_hopf).stability.l0; plot(real(evals),imag(evals),'*') axis([min(real(evals))-1 max(real(evals))+1 min(imag(evals))-1 max(imag(evals))+1]) print('-depsc','HopfChar.eps') saveas(gcf,'HopChar','fig') fprintf('Hopf bifurcation near point %d\n',ind_hopf); %% Continue Hopf bifurcation in two parameters [hbranch,suc]=SetupHopf(funcs,nontriv_eqs,ind_hopf,... 'contpar',[n_ind,tau_ind],'max_bound',[n_ind,21;tau_ind,4],'dir',n_ind,'step',1e-3); figure(2);clf hbranch=br_contn(funcs,hbranch,200); hbranch=br_rvers(hbranch); hbranch=br_contn(funcs,hbranch,200); print('-depsc','HopfCurvComp') figure(101); get_par=@(i)arrayfun(@(x)x.parameter(i),hbranch.point); plot(get_par(n_ind),get_par(tau_ind),'b'); hold on; axis([0 20 0 4]) print('-depsc','HopfCurv') %% Compute L1 coefficient % to find if Hopf bifurcation is supercritical (L1<0) or subcritical (L1>0) [L1,L1low]=HopfLyapunovCoefficients(funcs,hbranch); fprintf('maximal L1 coefficient along Hopf branch: %g\n',max(L1)); fprintf('max of error estimate for L1 coefficient: %g\n',norm(L1-L1low,'inf')); %% Branch off at Hopf bifurcation continue in n disp('Branch off at Hopf bifurcation'); fprintf('Initial correction of periodic orbits at Hopf:\n'); [per_orb,suc]=SetupPsol(funcs,nontriv_eqs,ind_hopf,... 'print_residual_info',1,'intervals',20,'degree',4,... 'max_bound',[contpar,21],'max_step',[contpar,0.1]); if ~suc error('MackeyGlassDemo:fail',... 'MackeyGlassDemo: initialization of periodic orbit failed'); end figure(1); hold on per_orb=br_contn(funcs,per_orb,200); per_orb=br_stabl(funcs,per_orb,0,1); [nunst_per,dom]=GetStability(per_orb,'exclude_trivial',true); %% Find period doubling bifurcations in two parameters ind_pd=find(diff(nunst_per)==1); disp(['For ',num2str(per_orb.point(ind_pd(1)).parameter(contpar)),' fprintf('max defect of Floquet multiplier at -1: %g\n',max(abs(floqpd1+1))); %% Branch off at period doubling % (Solutions at far end get inaccurate.) [per2,suc]=DoublePsol(funcs,per_orb,ind_pd(1)); if ~suc error('MackeyGlassDemo:fail',... 'MackeyGlassDemo: branching off at period doubling failed'); end figure(1); per2=br_contn(funcs,per2,200); per2=br_stabl(funcs,per2,0,1); [nunst_per2,dom,triv_defect]=GetStability(per2,'exclude_trivial',true) fprintf('max defect of Floquet multiplier at 1: %g\n',max(triv_defect)); %% Continue period doublings in two parameters for secondary PD ind_pd2=find(diff(nunst_per2)==1); disp(['For ',num2str(per2.point(ind_pd2(1)).parameter(contpar)),'0))=nan; plot(ParVal,AmpVec,'b--'); hold on; plot(ParVal,AmpVecSt,'b') [xma,yma]=df_measr(0,per2); ParVal2=br_measr(per2,xma); AmpVec2=br_measr(per2,yma); AmpVecSt2=AmpVec2; AmpVecSt2(find(nunst_per2>0))=nan; plot(ParVal2,AmpVec2,'b--'); hold on; plot(ParVal2,AmpVecSt2,'b') [xma,yma]=df_measr(0,per3); ParVal3=br_measr(per3,xma); AmpVec3=br_measr(per3,yma); AmpVecSt3=AmpVec3; AmpVecSt3(find(nunst_per3>0))=nan; plot(ParVal3,AmpVec3,'b--'); hold on; plot(ParVal3,AmpVecSt3,'b') xlim([0 20]) print('-depsc','StabAmpOrbs.eps') saveas(gcf,'StabAmpOrbs','fig') %% plot phase portriats nlist=[7 7.75 8.5 8.75 8.79 8.84 9.65 9.696 9.7056 9.7451 10 20]; for j=1:numel(nlist) n=nlist(j); figure(4+j); clf; if n>ParVal(1) ind=find(ParVal>n,1); point=per_orb.point(ind); point.parameter(contpar)=n; [newpoint,success]=p_correc(funcs,point,[],[],per_orb.method.point); tvals=linspace(0,1,1000); uvals=psol_eva(newpoint.profile,newpoint.mesh,tvals,newpoint.degree); tdelvals=mod(tvals-newpoint.parameter(tau_ind)/newpoint.period,1); udelvals=psol_eva(newpoint.profile,newpoint.mesh,tdelvals,newpoint.degree); plot(uvals,udelvals); hold on; end if n>ParVal2(1) ind=find(ParVal2>n,1); point=per2.point(ind); point.parameter(contpar)=n; [newpoint,success]=p_correc(funcs,point,[],[],per2.method.point); tvals=linspace(0,1,1000); uvals=psol_eva(newpoint.profile,newpoint.mesh,tvals,newpoint.degree); tdelvals=mod(tvals-newpoint.parameter(tau_ind)/newpoint.period,1); udelvals=psol_eva(newpoint.profile,newpoint.mesh,tdelvals,newpoint.degree); plot(uvals,udelvals,'k') end if n>ParVal3(1) ind=find(ParVal3>n,1); point=per3.point(ind); point.parameter(contpar)=n; [newpoint,success]=p_correc(funcs,point,[],[],per3.method.point); tvals=linspace(0,1,1000); uvals=psol_eva(newpoint.profile,newpoint.mesh,tvals,newpoint.degree); tdelvals=mod(tvals-newpoint.parameter(tau_ind)/newpoint.period,1); udelvals=psol_eva(newpoint.profile,newpoint.mesh,tdelvals,newpoint.degree); plot(uvals,udelvals,'r') end fil=['PhasePort',num2str(n),'.eps']; figfil=['PhasePort',num2str(n),'.fig']; print('-depsc',fil) saveas(gcf,figfil,'fig') end %% Plot of period doubling bifurcation x profiles bifsols={pd1sols,pd2sols,hbranch.point}; floqpd={floqpd1,floqpd2}; get_par=@(i,k)arrayfun(@(x)x.parameter(i),bifsols{k}); figure(3) clf; subplot(3,2,[1,2]); plot(get_par(n_ind,1),get_par(tau_ind,1),'.-',... get_par(n_ind,2),get_par(tau_ind,2),'.-',... get_par(n_ind,3),get_par(tau_ind,3),'.-'); legend('PD1','PD2','Hopf'); xlabel('\beta'); ylabel('\tau'); title(sprintf(['Hopf, 1st and 2nd period doubling in Mackey-Glass eqn, ',... ' n=%g, gamma=1'],n0)); grid on for k=1:2 subplot(3,2,2+k); hold on for i=1:length(bifsols{k}) plot(bifsols{k}(i).mesh*bifsols{k}(i).period,bifsols{k}(i).profile,'-'); end hold off box on grid on title(sprintf('PD%d: time profiles of period doubling',k)); xlabel('t'); ylabel('x'); subplot(3,2,4+k); plot(1:length(bifsols{k}),floqpd{k}+1,'.-'); grid on title(sprintf('PD%d: dist crit Floq mult from -1',k)); ylabel('error'); xlabel('point along branch'); end %% save save('MGresults.mat'); %% tidy up rmpath('C:/Program Files/MATLAB/dde_biftool_v3.1/ddebiftool/',... 'C:/Program Files/MATLAB/dde_biftool_v3.1/ddebiftool_extra_psol/',... 'C:/Program Files/MATLAB/dde_biftool_v3.1/ddebiftool_extra_psol/',... 'C:/Program Files/MATLAB/dde_biftool_v3.1/ddebiftool_utilities/'); diary('off')