
% Test_DuffingForced 
% Script to determine Frequency Response of a Duffing oscillator using
% numerical continuation rather than shooting.
% Author HGE Meijer, University of Twente, 23-2-2022

%% Initial simulation
% Initialize Matcont
init;
% We assume that we have an m-file setup for MatCont
fun=DuffingForced();

%Parameter values and initital condition
p=[2; 0.75; 0.1; 0.04; 3;0]; %Parameters alpha,beta,gamma,delta,omega
X0=[0; 0; 1; 0]; %Initial condition; third and fourth should not be zero
P0=num2cell(p);
%Running a first simulation for transients
[~,y] = ode45(fun{2},[0 500],[0; 0; 1; 0],[],P0{:});

% Now Take the final point and simulate over a little more than one period:
% 2pi/omega=2.1 (and add +.1=2.2) 
% To initialize the continuation it is better to have smaller steps as that
% adds to accuracy more than setting RelTol or AbsTol.
opt=odeset('MaxStep',.01);
x1 = y(end,:);
[t,y] = ode45(fun{2},[0 2.15],x1,opt,P0{:});

% Let's check the response is indeed periodic, and transients are gone.
figure(1)
plot(y(:,1),y(:,2));
%% Numerical continuation of the periodic orbit
% We provide the simulated response (y) to the initializer
% In the part below, you may need to switch Backward from 1 to 0, this is
% simply unpredictable across systems.
tolerance=1e-3;
ap=5;    % omega will be the active parameter
[x0,v0]=initOrbLC(@DuffingForced,t,y,p,ap,40,4,tolerance);
opt=contset;
opt=contset(opt,'MaxNumPoints',800);
opt=contset(opt,'Backward',1);
opt=contset(opt,'Singularities',1);
opt=contset(opt,'IgnoreSingularity',[1 2 4]); %Ignore everything except LPC
[xlc1,vlc1,slc1,hlc1,flc1]=cont(@limitcycle,x0,v0,opt);
% For the backward branch, we have to re-initialize as some global
% variables change during the continuation. Else, we might get the same
% branch. We need fewer points here, and we know there are no singularities
% on this part, so we switch off monitoring for bifurcations.
[x0,v0]=initOrbLC(@DuffingForced,t,y,p,ap,40,4,tolerance);
opt=contset(opt,'MaxNumPoints',100);
opt=contset(opt,'Backward',0);
opt=contset(opt,'Singularities',0);
[xlc2,vlc2,slc2,hlc2,flc2]=cont(@limitcycle,x0,v0,opt);

%% Visualization of Results
% Determine the amplitude
A1=max(xlc1(1:4:end-2,:));
A2=max(xlc2(1:4:end-2,:));

% A simple bifurcation diagram (amplitude as function of parameter)
figure(2);clf(2);hold on;
plot(xlc1(end,:),A1,xlc2(end,:),A2);

% A polished bifurcation diagram including stability
figure(3);clf(3);hold on;
ind1=(1:slc1(2).index); % split into separate parts
ind2=(slc1(2).index:slc1(3).index);
ind3=(slc1(3).index:slc1(4).index);
plot(xlc1(end,ind1),A1(ind1),xlc1(end,ind3),A1(ind3),xlc2(end,:),A2,'Color','blue');
plot(xlc1(end,ind2),A1(ind2),'Color','red','LineStyle','--');
xlabel('omega');ylabel('Amplitude');
xlim([0.5 4]);

