%% Transient Behavior and Location of Characteristic Roots in the Complex Plane
%
%% Initial Value Problem
%
% Consider the system described by the following second-order,
% ordinary differential equation and associated initial conditions:
%
% $$\frac{d^2x}{dt^2} + 2\zeta\omega_n\frac{dx}{dt} + \omega_n^2x = 0; x(0)=x0,\frac{dx}{dt}(0)=0 $$
%
% We would like to examine the behavior of the free response of the
% system as the damping ratio changes. This will be done first by solving
% the initial value problem for various damping ratios and plotting the
% results. Then, we will explore changes in the damping ratio affect the
% characteristic roots. (Recall that the characteristic roots dictate the
% free-response behavior.)
%% Set physical and derived parameters.
%
% To enable numerical solution of the differential equation, we will set
% specific values for the mass and (linear) spring rate. We will then use
% these to derive the (circular) natural frequency.
%
param.m = 1; % Mass [kg]
param.k = 1; % Spring rate [N/m]
param.wn = sqrt(param.k/param.m); % (Circular) natural frequency [rad/s]
%% Set initial conditions.
%
% To solve our 2nd order differential equation, we need two initial
% conditions: one on position, the other on velocity. Moreover, recall
% from the accompanying PowerPoint slides that we have represented our
% 2nd order linear ODE into a system of 1st order ODE's amenable to
% solution using MATLAB. Having defined position and velocity as our
% system states, we must define an initial state vector.
x0 = 1; % Position [m]
x0_dot = 0; % Velocity [m/s]
z0 = [x0 ; x0_dot]; % Initial condition on state vector
%% Set simulation parameters.
%
% We will define the duration of our simulation to be approximately 4
% periods. To accommodate changes in physical parameters, we will
% parameterize the duration by the natural frequency. We then define the
% associated time step and vector of simulation time. Note that the time
% step dt serves to define the times at which output will be provided, not
% the solver time step.
tend = 4*round((2*pi)/param.wn); % Simulate for approximately 4 periods [s]
dt = round((2*pi)/param.wn)/1000; % Time step [s]
tsim = 0:dt:tend; % Vector of simulation times [s]
%% Undamped system.
%
% We first consider an undamped system, characterized by a damping ratio of
% zero. Physically, this represents a system with zero dissipation of
% energy. We will plot the free response to confirm the expected behavior.
%
% To solve the system of 1st order ODE's, we will use the baseline ode45
% solver. The solver requires a vector of state derivatives; in this
% case, we call upon a function to define them.
param.zeta = 0; % Damping ratio
[t,z] = ode45(@(t,z) odesys( t,z,param ),tsim,z0);
x(:,1) = z(:,1); % Save solution for future plotting.
plot(t,x(:,1),'linewidth',2);grid on;hold all
ylim([-1 1])
xlabel('Time [s]');
ylabel('Position [m]')
title('Free Response: Undamped 2nd Order System');
%% Consider several damping ratios.
%
% Now we will consider four distinct values of the damping ratio
% corresponding to four distinct types of system behavior: undamped,
% underdamped, critically damped, and overdamped. These terms refer to the
% behavior of the system as it returns to its equilibrium state. We will
% plot and compare the free responses for these cases.
zetavals = [0 1/sqrt(2) 1 1.5]; % Damping ratios to consider.
case1 = 'Undamped: zeta = 0' ;
case2 = 'Underdamped: zeta = 0.707' ;
case3 = 'Critically Damped: zeta = 1' ;
case4 = 'Overdamped: zeta = 1.5' ;
for i = 2: length(zetavals)
param.zeta = zetavals(i); % Set current damping ratio.
[t,z] = ode45(@(t,z) odesys( t,z,param ),tsim,z0); % Compute free response.
x(:,i) = z(:,1); % Save solution for future plotting.
end;
figure % Plot results
plot(t,x,'linewidth',2);grid on
xlabel('Time [s]');
ylabel('Position [m]')
title('Free Response of 2nd Order System for Various Damping Ratios');
legend(case1,case2,case3,case4,0)
%% Track characteristic roots as damping ratio increases.
%
% Recalling that the free response of a system is governed by its
% characteristic roots, we will examine the changes in these roots as the
% damping ratio changes over a more finely spaced set of values. Since the
% roots may have real and imaginary parts, we will plot them in the complex
% plane. For each value of the damping ratio, we will compute the
% characteristic roots as the eigenvalues of the associated state matrix.
zeta_range = 0:0.01:1.5; % Damping ratio: undamped to overdamped.
S = zeros(2,length(zeta_range)); % Initialize characteristic roots array.
for i = 1:length(zeta_range)
S(:,i) = eig([0 1; -param.wn^2 -2*zeta_range(i)*param.wn]); % Characteristic roots.
end;
figure % Plot characteristic roots in complex plane.
plot(S,'bx','linewidth',2,'markersize',5); grid on;
set(gca,'YAxisLocation','right')
xlabel('Real Axis','Position',[-1.5 -0.1 1])
ylabel('Imaginary Axis')
title('Loci of Characteristic Roots in Complex Plane for Various Damping Ratios')
%% Show the "odesys" function.
%
% Here we record the function used to define the state matrix
% corresponding to our system of ordinary differential equations.
type odesys.m