% ### EXduffing.m ###       09.19.14 CB (revised 2020.01.30)

% Purpose: Numerically integrate the driven damped duffing oscillator 
%   m*x''+ b*x' + k*x + c*x^3 = A*sin(wt)

% Notes
% o Default params taken from wikipedia page re "Duffing Equation"
% (presumably as they point to a chaotic regime)

clear
% -----------------------------------------------------
% User input (Note: All paramters are stored in a structure)
P.y0(1) = -1.0;   % initial position [m]
P.y0(2) = 0.0;   % initial velocity [m/s]
P.b= 0.02;  % damping coefficient [kg/s]
P.k= 1.0;   % stiffness [N m]
P.m= 1;    % mass [kg]
P.c= 5;       % nonlinear stiffness term
% --- (sinusoidal) driving term
P.A= 8.0;   % amplitude [N] (set to zero to turn off) {8}
fD= 0.5/(2*pi);  % freq. (Hz) [expressed as fraction of resonant freq.] {0.5}
% --- Integration limits
P.t0 = 0.0;   % Start value
P.tf = 100.0;   % Finish value
P.dt = 0.01;  % time step
% ----------------------------------------------------------------------
% +++ bookeeping + display some basic derived quantities
P.wr= 2*pi*fD;  % convert to angular freq.
disp(sprintf('Resonant frequency ~%g [Hz]', sqrt(P.k/P.m)/(2*pi)));
Q = (sqrt(P.k/P.m))/(P.b/P.m);  % quality factor
disp(sprintf('Q-value = %g', Q));
% +++ use built-in ode45 to solve
[t y] = ode45('DUFFfunction', [P.t0:P.dt:P.tf],P.y0,[],P);

% ------------------------------------------------------
% --- visualize timecourse
figure(1); clf; plot(t,y(:,1)); hold on; grid on;
xlabel('t [s]');    ylabel('x(t) [m]')
% --- Phase plane
figure(2); clf;
plot(y(:,1), y(:,2)); hold on; grid on;
xlabel('x [m]');    ylabel('dx/dt [m/s]')